收藏 分享(赏)

中南大学材料院matlab程序设计实践.doc

上传人:精品资料 文档编号:8231283 上传时间:2019-06-15 格式:DOC 页数:33 大小:832.31KB
下载 相关 举报
中南大学材料院matlab程序设计实践.doc_第1页
第1页 / 共33页
中南大学材料院matlab程序设计实践.doc_第2页
第2页 / 共33页
中南大学材料院matlab程序设计实践.doc_第3页
第3页 / 共33页
中南大学材料院matlab程序设计实践.doc_第4页
第4页 / 共33页
中南大学材料院matlab程序设计实践.doc_第5页
第5页 / 共33页
点击查看更多>>
资源描述

1、中南大学MATLAB 程序设计实践班级:学号:姓名:一、 MATLAB 程序设计实践Matlab 基础表示多晶体材料织构的三维取向分布函数(f f( 1, 2) )是一个非常复杂的函数,难以精确的用解析函数表达,通常采用离散空间函数值来表示取向分布函数,Data.txt 是三维取向分布函数的一个实例。由于数据量非常大,不便于分析,需要借助图形来分析。请你编写一个 matlab 程序画出如下的几种图形来分析其取向分布特征:(1)用 Slice 函数给出其整体分布特征;(2)用 pcolor 或 contour 函数分别给出(20, 5, 10, 15, 20, 25, 30, 35 90)切面上

2、 f 分布情况( 需要用到 subplot 函数) ;(3) 用plot 函数给出沿 取向线(1=090, 45,20)的 f 分布情况。备注:data.txt 数据格式说明数据说明部分,与作图无关 此方向表示 f 随着1 从 0,5,10,15, 20 到 90 的变化而变化此方向表示 f 随着 从 0,5,10,15, 20 到90 的变化而变化表示以下数据为20 的数据,即f( 1, ,0)解:(1)将文件 Data.txt 内的数据按照要求读取到矩阵 f(phi1,phi,phi2)中,代码如下:fid=fopen(data.txt);for i=1:18tline=fgetl(fid

3、);endphi1=1;phi=1;phi2=1;line=0;f=zeros(19,19,19);while feof(fid)tline=fgetl(fid);data=str2num(tline);line=line+1;if mod(line,20)=1phi2=(data/5)+1;phi=1;elsefor phi1=1:19f(phi1,phi,phi2)=data(phi1);endphi=phi+1;endendfclose(fid);将以上代码保存为 readtext.m 文件并在 MATLAB 中运行,运行结果如下图所示:将以下代码保存为 code1_1.m 文件:fop

4、en(readtext.m);readtext;x,y,z=meshgrid(0:5:90,0:5:90,0:5:90);slice(x,y,z,f,45,90,45,90,0,45)运行结果如右图所示:(2)将以下代码保存为 code1_2_1.m 文件:fopen(readtext.m);readtext;for i=1:19subplot(5,4,i)pcolor(f(:,:,i)End运行结果如右图所示:将以下代码保存为 code1_2_2.m 文件:fopen(readtext.m);readtext;for i=1:19subplot(5,4,i)contour(f(:,:,i)e

5、nd运行结果如右图所示:(3)1=090, 45,20 所对应的 f(1,2)即为 f(:,10,1)。将以下代码保存为 code1_3.m 文件:fopen(readtext.m);readtext;plot(0:5:90,f(:,10,1),-bo)text(60,6,phi=45 phi2=0)运行结果如下图所示:二 MATLAB 程序设计实践科学计算(07)、编程实现以下科学计算算法,并举一例应用之。 (参考书籍精通科学计算 ,王正林等著,电子工业出版社,年)“傅里叶逼近”解:(1)用法说明:对于连续周期函数,只要计算出其傅里叶展开级数即可,在 Matlab 中编程实现的连续函数的傅里

6、叶逼近法函数为:FZZ。功能:用傅里叶级数逼近已知的连续周期函数。调用格式:A0,A,B=FZZ(func,T,n).其中,func 为已知函数;T 为已知函数的周期;N 为展开级数的项数;A0 为展开后的常数项;A 为展开后的余弦项系数;B 为展开后的正弦项系数。(2)源程序代码:function A0,A,B=FZZ(func,T,n)syms t;func=subs(sym(func),sym(x),sym(t);A0=int(sym(func),t,-T/2,T/2)/T;for(k=1:n)A(k)=int(func*cos(2*pi*k*t/T),t,-T/2,T/2)*2/T;A

7、(k)=vpa(A(k),4);B(k)= int(func*sin(2*pi*k*t/T),t,-T/2,T/2)*2/T;B(k)= vpa(B(k),4);end(3)输出结果:傅里叶逼近应用实例。用傅里叶级数(取 5 项)逼近函数 x,输出系数值。 A0,A,B=FZZ(x,2*pi,5)A0 =0A =0.,0.,0.,0.,0.B =2.,-1.,.6667, -.5000,.4000(4)流程图:开始调用傅里叶级数A0,A,B=FZZ(func,T,n)输入已知函数 func根据 func 计算并输出傅里叶级数常数项 A0计算并输出傅里叶级数的各余弦项系数A(k) ,精度为 4

8、位小数求出并输出傅里叶级数的各正弦项系数B(k)精度为 4 位小数结束流程图:开始调用傅里叶级数A0,A,B=FZZ(func,T,n)输入已知函数func=sinx()根据 func=sinx 计算并输出傅里叶级数常数项 A0计算并输出傅里叶级数的各余弦项系数A(1)到 A(5)、编程解决以下科学计算问题。1)解:(1)算法说明:系统的受力平衡表达式:f-c*v-K*x=M*a; (1)求出并输出傅里叶级数的各正弦项系数 B(1)到B(5)结束其中 x 为位移,v 为速度,a 为加速度用 dx/dt 代替 v,得:f-4*dx/dt-100*x=1*d2x/dt2; (2)运用函数分段求解位

9、移的表达式,然后用 plot 画图(2)源程序代码: dsolve(D2x=t/0.015-4*Dx-100*x,x(0)=0,Dx(0)=0)ans =-23/900*exp(-2*t)*sin(4*6(1/2)*t)*6(1/2)+2/75*exp(-2*t)*cos(4*6(1/2)*t)-2/75+2/3*t t1=0:0.01:0.15; x1=-23/900*exp(-2.*t1).*sin(4*6(1/2).*t1).*6(1/2)+2/75*exp(-2.*t1).*cos(4*6(1/2).*t1)-2/75+2/3.*t1; subs(x1,t1,0.15)ans =Col

10、umns 1 through 11 0 0.0000 0.0001 0.0003 0.0007 0.0013 0.0022 0.0035 0.0051 0.0071 0.0096Columns 12 through 16 0.0125 0.0160 0.0199 0.0243 0.0292 subs(diff(x1),t1,0.15)ans =Columns 1 through 11 0.0000 0.0001 0.0002 0.0004 0.0006 0.0009 0.0013 0.0016 0.0020 0.0025 0.0029Columns 12 through 15 0.0034 0

11、.0039 0.0044 0.0049 dsolve(D2x=10-4*Dx-100*x,x(0.15)=0.0292,Dx(0.15)=0.0049)ans =-1/240000*exp(-2*t)*sin(4*6(1/2)*t)*exp(3/10)*(1367*cos(3/5*6(1/2)*6(1/2)+16992*sin(3/5*6(1/2)+exp(-2*t)*cos(4*6(1/2)*t)*(1367/240000*exp(3/10)*6(1/2)*sin(3/5*6(1/2)-177/2500*exp(3/10)*cos(3/5*6(1/2)+1/10 t2=0.15:0.1:1.

12、2; x2=-1/240000*exp(-2.*t2).*sin(4*6(1/2).*t2)*exp(3/10)*(1367*cos(3/5*6(1/2)*6(1/2)+16992*sin(3/5*6(1/2)+exp(-2.*t2).*cos(4*6(1/2).*t2)*(1367/240000*exp(3/10)*6(1/2)*sin(3/5*6(1/2)-177/2500*exp(3/10)*cos(3/5*6(1/2)+1/10; plot(t1,x1,t2,x2)(3)运行结果:(4)流程图: 开始根据受力情况写出方程(1)将方程(1)化成方程(2)用 dsolve 函数求解前0.1

13、5s 的微分方程符号解,解出位移方程 x1利用 x(0.15) 、v(0.15)和表达式(2)用 dsolve 求解位移方程 x2根据方程 x1 求 0.15s 时的位移x(0.15)和速度 v(0.15)2)解:=(1-x)*u/(1+x)2 + 1/(1+x)2 (2)22)1(xud用 plot 画出 x1 和 x2的图形,即 1.2s 以内的波形结束由表达式(2)知,这是一个椭圆型偏微分方程-div(c*grad(u)+a*u=f所以得到系数为 c=-1;a=(x-1)./(1+x).2; f=1./(1+x).2操作步骤如下:1) 先打开 pdetool,在 command 窗口输入

14、 pdetool2) 设定矩形区域如下3)设定边界条件,假设 y 与 x 的边界条件一致4)设定矩形区域边界 如下:在 y 方向默认 u=05)设定椭圆偏微分方程系数c=-1;a=(x-1)./(1+x).2;f=1./(1+x).26)画三角形并细化7)点击工具栏中的三维画图标准,画图(2)整个过程产生的代码:function pdemodelpde_fig,ax=pdeinit;pdetool(appl_cb,1);set(ax,DataAspectRatio,1 1.5 1);set(ax,PlotBoxAspectRatio,1 0.66666666666666663 1);set(a

15、x,XLim,-0.5 1.5);set(ax,YLim,-0.5 1.5);set(ax,XTickMode,auto);set(ax,YTickMode,auto);% Geometry description:pderect(0 1 1 0,R1);set(findobj(get(pde_fig,Children),Tag,PDEEval),String,R1)% Boundary conditions:pdetool(changemode,0)pdesetbd(4,.dir,.1,.1,.1)pdesetbd(3,.dir,.1,.1,.0)pdesetbd(2,.dir,.1,.1,

16、.0.5)pdesetbd(1,.dir,.1,.1,.0)% Mesh generation:setappdata(pde_fig,Hgrad,1.3);setappdata(pde_fig,refinemethod,regular);pdetool(initmesh)pdetool(refine)% PDE coefficients:pdeseteq(1,.-1.0,.(x-1)./(1+x).2,.1./(1+x).2,.1.0,.0:10,.0.0,.0.0,.0 100)setappdata(pde_fig,currparam,.-1.0 ;.(x-1)./(1+x).2;.1./(

17、1+x).2 ;.1.0 )% Solve parameters:setappdata(pde_fig,solveparam,.str2mat(0,1872,10,pdeadworst,.0.5,longest,0,1E-4,fixed,Inf)% Plotflags and user data strings:setappdata(pde_fig,plotflags,1 1 1 1 1 1 1 1 1 1 0 1 1 0 1 0 0 1);setappdata(pde_fig,colstring,);setappdata(pde_fig,arrowstring,);setappdata(pde_fig,deformstring,);

展开阅读全文
相关资源
猜你喜欢
相关搜索

当前位置:首页 > 企业管理 > 管理学资料

本站链接:文库   一言   我酷   合作


客服QQ:2549714901微博号:道客多多官方知乎号:道客多多

经营许可证编号: 粤ICP备2021046453号世界地图

道客多多©版权所有2020-2025营业执照举报