中南大学材料专业MATLAB实践-4号

1一、《MATLAB程序设计实践》Matlab基础班级：学号：姓名：表示多晶体材料织构的三维取向分布函数（f＝f（φ1，φ，φ2））是一个非常复杂的函数，难以精确的用解析函数表达，通常采用离散空间函数值来表示取向分布函数，Data.txt是三维取向分布函数的一个实例。由于数据量非常大，不便于分析，需要借助图形来分析。请你编写一个matlab程序画出如下的几种图形来分析其取向分布特征：（1）用Slice函数给出其整体分布特征；2（2）用pcolor或contour函数分别给出(φ2＝0,5,10,15,20,25,30,35…90)切面上f分布情况(需要用到subplot函数)；(3)用plot函数给出沿α取向线(φ1=0~90，φ＝45，φ2＝0)的f分布情况。3备注：data.txt数据格式说明将文件Data.txt内的数据按照要求读取到矩阵f(phi1,phi,phi2)中，代码如下：fid=fopen('data.txt');%读取数据文件Data.txtfori=1:18tline=fgetl(fid);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;ifmod(line,20)==1phi2=(data/5)+1;phi=1;elseforphi1=1:19数据说明部分，与作图无关此方向表示f随着φ1从0,5,10,15,20…到90的变化而变化此方向表示f随着φ从0,5,10,15,20…到90的变化而变化表示以下数据为φ2＝0的数据，即f（φ1，φ，0）4f(phi1,phi,phi2)=data(phi1);endphi=phi+1;endendfclose(fid);建立新的脚本，将以上代码保存为m文件，命名为readtext.m，在MATLAB运行后。其变量f结果为：(1).用Slice函数给出其整体分布特征程序代码：fopen('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])%运用slice函数绘制图形将程序代码保存为code1_1.m文件，其运行结果：(2).○1.用subplot函数和pcolor函数给出(φ2＝0,5,10,15,20,25,30,35…90)切面上f分布情况；5程序代码：fopen('readtext.m');readtext;fori=1:19subplot(5,4,i)pcolor(f(:,:,i))%运用pcolor函数绘制图形end将上述代码保存为code1_2_1.m文件，其运行结果：○2.使用subplot函数和contour函数给出(φ2＝0,5,10,15,20,25,30,35…90)切面上f分布情况；程序代码：fopen('readtext.m');%运用contour函数绘制图形readtext;fori=1:19subplot(5,4,i)contour(f(:,:,i))end将上述代码保存为code1_2_2.m文件，运行结果：6（3）.用plot函数给出沿α取向线(φ1=0~90，φ＝45，φ2＝0)的f分布情况。程序代码：fopen('readtext.m');readtext;plot([0:5:90],f(:,10,1),'-bo')%运用plot函数绘制图形text(60,6,'\phi=45\phi2=0')将上述代码保存为code1_3.m文件，运行结果：7二《MATLAB程序设计实践》科学计算（04）班级１、编程实现以下科学计算算法，并举一例应用之。（参考书籍《精通ＭＡＬＡＢ科学计算》，王正林等著，电子工业出版社，２００９年）“切比雪夫逼近”算法说明：当一个连续函数定义在区间[-1,1]上时，它可以展开成切比雪夫级数。即：0()()nnnfxfTx其中()nTx为n次切比雪夫多项式，具体表达可通过递推得出：0()1Tx，1()Txx11()2nnnTxxTxTx它们之间满足如下的正交关系：1210,()(),021,0nmnmTxTxdxnmxnm在实际应用中，可根据所需的精度来截取有限的项数，切比雪夫级数中的系数由下式决定：10211()1fxfdxx81212()()1nTxfxfdxx流程图：程序源代码（m文件）：functionf=Chebyshev(y,k,x0)%用切比雪夫多项式逼近已知函数开始定义sym型变量t以及切比雪夫多项式矩阵T(n)，并规定前两项为1，t根据前两项多项式计算其系数c(1)，c(2)由第一项和第二项多项式及其系数计算一级切比雪夫逼近值i从3渐渐增加到k+1由T(i-1)和T(i-2)计算T(i)由T(i)计算其系数c(i)由T(i)和c(i)计算i级切比雪夫逼近值。ik+1？否是输入三项?是求出逼近x0处六位精度下的值是输出k项切比雪夫多项式否结束9%已知函数：y%逼近已知函数所需项数：k%逼近点的x坐标：x0%求得的切比雪夫逼近多项式或在x0处的逼近值：fsymst;%定义sym型变量tT(1:k+1)=t;T(1)=sym('1');T(2)=t;%定义切比雪夫多项式矩阵T(n)，并规定前两项为1，tc(1:k+1)=sym('0');c(1)=int(subs(y,findsym(sym(y)),sym('t'))*T(1)/sqrt(1-t^2),t,-1,1)/pi;c(2)=2*int(subs(y,findsym(sym(y)),sym('t'))*T(2)/sqrt(1-t^2),t,-1,1)/pi;%根据前两项多项式计算其系数c(1)，c(2)f=c(1)+c(2)*t;%由第一项和第二项多项式及其系数计算切比雪夫逼近值fori=3:k+1%从第3项到第k+1项进行循环T(i)=2*t*T(i-1)-T(i-2);%计算第i项切比雪夫多项式矩阵T(i)c(i)=2*int(subs(y,findsym(sym(y)),sym('t'))*T(i)/sqrt(1-t^2),t,-1,1)/2;%计算第i项系数f=f+c(i)*T(i);%计算切比雪夫逼近值f=vpa(f,6);if(i==k+1)if(nargin==3)f=subs(f,'t',x0);%输入三项，则输出x0处六级精度的切比雪夫逼近值endf=vpa(f,6);endend根据代码可知，在MATLAB中编程实现的切比雪夫逼近法函数为：Chebyshev。功能：用切比雪夫多项式逼近已知函数。调用格式：Chebyshev(y,k,x0)f其中，y为已知函数；k为逼近已知函数所需项数；f是求得的切比雪夫逼近多项式或在x0处的逼近值。应用实例：切比雪夫应用实例。用切比雪夫公式（取6项）逼近函数12x，并求当x=1时的函数值。解：利用程序求解方程，在MATLAB命令窗口中输入：Chebyshev('1/(2-x)',5)%调用创建的函数，输出切比雪夫多项式的510个项结果：ans=0.277013*t+0.0647768*t*(2.0*t^2-1.0)-0.00501051*t*(2.0*t*(t-2.0*t*(2.0*t^2-1.0))+2.0*t^2-1.0)-0.0186995*t*(t-2.0*t*(2.0*t^2-1.0))+0.24175*t^2+0.456475再在MATLAB命令窗口中输入：Chebyshev('1/(2-x)',5,1)%调用创建的函数，输出当x=1时的函数值逼近值结果：ans=1.06372在MATLAB命令窗口计算符号表达式x=1的极限：limit(sym('1/(2-x)'),'x',1)%计算其真值，并与其进行比较结果：ans=1输出结果：绝对误差：1.06372-1=0.06372；相对误差：0.06372/1*100%=6.372%２、编程解决以下科学计算问题。111)问题分析：1.相关参数字母设定：如图，由在P的作用，点A移动到点A’，设在水平方向移动距离为dx，在竖直方向移动距离为dy，P与水平的夹角为α。设各杆的轴力为Ni(i=1,2,3,…n)；各杆与水平面夹角分别αi(i=1,2,3,…n)；各杆的变化长度分别为dLi(i=1,2,3,…n)；2.解题思路：（1）取一根杆分析，由题意可知每根杆的伸长量为dx和dy在杆上的分量；故可得方程○1：(2)将各轴力和外力P水平（x）竖直（y）正交分解，由受力平衡可得：x方向有方程○2：12y方向有方程○3：（3）联立上述○1○2○3可得n+2等式，可解得n个轴力，以及A点位移dx和dy，联立的线性方程组如下：...（4）建立[Pcosα,Psinα，0，0，0…0]’的常数矩阵，以及如下所示的系数矩阵：（5）再用求逆法求解此线性方程组，即用常数矩阵除以系数矩阵，得出结果。程序源代码：13Ei=input('请输入各杆的刚度:(注意用[]括起来)');%输入刚度矩阵EiLi=input('请输入各杆的长度:(注意用[]括起来)');%输入杆的长度矩阵LiAi=input('请输入各杆的横截面积:(注意用[]括起来)');%输入杆的横截面积矩阵Aiai=input('请输入各杆与水平面的夹角:(注意用[]括起来)');%输入杆与水平面的夹角矩阵aiP=input('请输入外力P:');%输入外加力Pa=input('请输入P与水平面的夹角:');%输入外加力P与x的夹角n1=length(Ei);n2=length(Li);n3=length(Ai);n4=length(ai);if(n1~=n2|n1~=n3|n1~=n4)disp('输入数据错误')elsen=n1;end%判断赋值维度是否一致Ki=Li./(Ei.*Ai);C=zeros(n+2,1);C(1,1)=P*cos(a);C(2,1)=P*sin(a);C(3:n+2,1)=zeros(n,1);%建立方程组等号右边常数的矩阵D=zeros(n+2,n+2);D(1,:)=[cos(ai),0,0];D(2,:)=[sin(ai),0,0];for(i=1:n)D(i+2,i)=Ki(i);endD(3:n+2,n+1)=(-cos(ai));D(3:n+2,n+2)=(-sin(ai));%建立方程组系数矩阵x=D\C;x=x';%求解该线性方程组，得出位移以及每根杆的轴力disp('节点在x、y方向上的位移分别：')x(n+1:n+2)disp('各杆的轴力分别为：')x(1:n)%输出结果流程图：开始结束14应用实例：设三根杆组成的支架如下图所示，挂一重物P=3000N。设L=3m，各杆的横截面积分别为：A1=15010-6m2，A2=20010-6m2，A3=30010-6m2，材料的弹性模量均为E=200109N/m2，求各杆所受力的大小以及C点位移解：将前面的程序源代码保存为main.m文件运行，根据提示依次输入题中数据：[200e9,200e9,200e9];[3/sin(pi/3),3/sin(pi/2),3/sin(pi/4)];[150e-6,200e-6,300e-6];[pi/3,pi/2,3*pi/4];3000;0;输出结果：节点在x、y方向上的位移分别：输入所需数据：各杆的刚度Ei、横截面积Ai、长度Li，以及夹角ai，外加力P以及其与x方向夹角输入Ei，Ai，Li，αi维度是否一致建立系数矩阵D，以及常数矩阵C用求逆法解线性微分方程输出A点位移，以及各杆轴力输出：输入数据错误是否15ans=1.0e-03*0.33990.0420各杆的轴力分别为：ans=1.0e+03*1.78650.5595-2.97