中南大学MATLAB程序设计实践学长有爱奉献，下载填上信息即可上交，没有下载券的自行百度。

所需m文件照本文档做即可，即新建（FILE）→脚本（NEW-Sscript）→复制本文档代码→运行（会跳出保存界面，文件名默认不要修改，保存）→结果。

第一题需要把数据文本文档和m文件放在一起。

全部测试无误，放心使用。

本文档针对做牛顿法求非线性函数题目的同学，当然第一题都一样，所有人都可以用。

←记得删掉这段话班级：学号：姓名：一、《MATLAB程序设计实践》Matlab基础表示多晶体材料织构的三维取向分布函数（f＝f（φ1，φ，φ2））是一个非常复杂的函数，难以精确的用解析函数表达，通常采用离散空间函数值来表示取向分布函数，是三维取向分布函数的一个实例。

由于数据量非常大，不便于分析，需要借助图形来分析。

请你编写一个matlab程序画出如下的几种图形来分析其取向分布特征：（1）用Slice函数给出其整体分布特征；"~（2）用pcolor或contour函数分别给出(φ2＝0, 5, 10, 15, 20, 25, 30, 35 … 90)切面上f分布情况(需要用到subplot函数)；(3) 用plot函数给出沿α取向线(φ1=0~90，φ＝45，φ2＝0)的f分布情况。

(备注：数据格式说明解：（1）（（2）将文件内的数据按照要求读取到矩阵f(phi1,phi,phi2)中，代码如下：fid=fopen('');for i=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;数据说明部分，与作图无关此方向表示f随着φ1从0,5,10,15,20 …到90的变化而变化此方向表示f随着φ从0,5,10,15, 20 …到90的变化而变化表示以下数据为φ2＝0的数据，即f（φ1，φ，0）if mod(line,20)==1phi2=(data/5)+1;phi=1;else~for phi1=1:19f(phi1,phi,phi2)=data(phi1);endphi=phi+1;endendfclose(fid);。

将以上代码保存为在MATLAB中运行，运行结果如下图所示：！将以下代码保存为文件：fopen('');,readdata;[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））（3）将以下代码保存为文件：fopen('');readdata;for i=1:19subplot(5,4,i)pcolor(f(:,:,i))nd)运行结果如下图所示：|将以下代码保存为文件：fopen('');readdata;for i=1:19…subplot(5,4,i)contour(f(:,:,i))end运行结果如下图所示：（3）φ1=0~90，φ＝45，φ2＝0所对应的f(φ1,φ,φ2)即为f(:,10,1)。

将以下代码保存为文件：fopen('');…readdata;plot([0:5:90],f(:,10,1),'-bo')text(60,6,'\phi=45 \phi2=0')运行结果如下图所示：》#；二 《MATLAB 程序设计实践》科学计算（24）班级： 学号： 姓名： >１、编程实现以下科学计算算法，并举一例应用之。

（参考书籍《精通ＭＡＬＡＢ科学计算》，王正林等著，电子工业出版社，２００９年）“牛顿法非线性方程求解”解：弦截法本质是一种割线法，它从两端向中间逐渐逼近方程的根；牛顿法本质上是一种切线法，它从一端向一个方向逼近方程的根，其递推公式为：-=+n n x x 1)()('n n x f x f 初始值可以取)('a f 和)('b f 的较大者，这样可以加快收敛速度。

和牛顿法有关的还有简化牛顿法和牛顿下山法。

在MATLAB 中编程实现的牛顿法的函数为：NewtonRoot 。

功能：用牛顿法求函数在某个区间上的一个零点。

/调用格式：root=NewtonRoot )(```eps b a f其中，f为函数名；a为区间左端点；b为区间右端点eps为根的精度；root为求出的函数零点。

，—牛顿法的matlab程序代码如下：|function root=NewtonRoot(f,a,b,eps)%牛顿法求函数f在区间[a,b]上的一个零点%函数名:f%区间左端点:a%区间右端点:b%根的精度:eps》%求出的函数零点:rootif(nargin==3)eps=;endf1=subs(sym(f),findsym(sym(f)),a);f2=subs(sym(f),findsym(sym(f)),b);if (f1==0)root=a;【endif (f2==0)root=b;endif (f1*f2>0)disp('两端点函数值乘积大于0 ！');return;"elsetol=1;fun=diff(sym(f)); %求导数fa=subs(sym(f),findsym(sym(f)),a);fb=subs(sym(f),findsym(sym(f)),b);dfa=subs(sym(fun),findsym(sym(fun)),a);dfb=subs(sym(fun),findsym(sym(fun)),b);if(dfa>dfb) %初始值取两端点导数较大者、root=a-fa/dfa;elseroot=b-fb/dfb;endwhile(tol>eps)r1=root;fx=subs(sym(f),findsym(sym(f)),r1);dfx=subs(sym(fun),findsym(sym(fun)),r1); %求该点的导数值、root=r1-fx/dfx; %迭代的核心公式tol=abs(root-r1);endend例：求方程3x^2－exp（x）＝0的一根解：在MATLAB命令窗口输入：>> r=NewtonRoot('3*x^2-exp(x)',3,4)输出结果：X=流程图：…２、编程解决以下科学计算问题。

1)解：这个方程可用下列步骤来解（1）用eig 函数求出矩阵K-λ M 的特征值L 和特征向量U ，U 和L 满足⎥⎦⎤⎢⎣⎡===21''00****λλU K U L IU M U（2）在原始方程Mx+Kx=0两端各乘以'U 及在中间乘以对角矩阵'U *U ，得<'U *M*'U *U*'x +'U *K*'U *x=0作变量置换x U x x U z z z **'21'21=⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡=，得0*''=+z L z这是一个对角矩阵方程，即可把它分两个方程：022''211''1=+=+z z z z λλ这意味着两种振动模态可以解耦，令i i λω=2，i ω是第i 个模式的固有频率（i=1,2）。

(3)由上述的解耦模态中，给出初始条件0x ，0d x ，化为0z ，0d z ，即可求出其分量1z ，z 。

设位置和速度的初始条件分别为[]Tx x x 02010=，[]Td d d x x x 02010=，则这三个步骤得到的最后公式为[][][])sin 1cos .()(0021t Mx u t Mx u u t x i d T i ii Ti i i ωωω+=∑=【系统解耦的振动模态的MATLAB 代码如下： function erziyoudu()流程图：%输入各原始参数m1=input('m1=');m2=input('m2='); %质量k1=input('k1=');k2=input('k2='); %刚度%输入阻尼系数%c1=input('c1=');c2=input('c2=');%给出初始条件及时间向量x0=[1;0];xd0=[0;-1];.tf=50; %步数dt=; %步长%构成二阶参数矩阵M=[m1,0;0,m2];(K=[k1+k2,-k2;-k2,k2];C=[c1+c2,-c2;-c2,c2];%构成四阶参数矩阵A=[zeros(2,2),eye(2);-M\K,-M\C];】%四元变量的初始条件y0=[x0;xd0];%设定计算点，作循环计算for i=1:round(tf/dt)+1t(i)=dt*(i-1);y(:,i)=expm(A*t(i))*y0;%循环计算矩阵指数end%按两个分图绘制x1、x2曲线！subplot(2,1,1),plot(t,y(1,:)),gridxlabel('t'),ylabel('y');subplot(2,1,2),plot(t,y(2,:)),gridxlabel('t'),ylabel('y');\运行M文件，依下图所示在MATLAB命令窗口中输入数据：~即可得出该振动的两种模态$*2)解：第一步，在MATLAB命令窗口输入命令pdetool打开工具箱，调整x坐标范围为[0 4],y 坐标范围为[0 3].通过options选项的Axes Linits设定如下图所示。

第二步，设定矩形区域。

点击工具箱栏中的按钮“”，拖动鼠标画出一矩形，并双击该矩形，设定矩形大小，如下图所示。

第三步，设边界条件。

点击工具栏中的按钮“”，并双击矩形区域的相应的边线在弹出的对话框中设定边界条件。

如下图所示，分别为各边框的边界条件。

第四步，设定方程。

单击工具栏中的按钮“”，在PDE模式下选择方程类型，如下图所示，并在其中设定参数。

第五步，单击工具栏中的按钮“”，拆分区域为若干子区域，如下图所示。

第六步，单击工具栏中的按钮“”，将子区域细化，从而保证结果更精确，如下图所示。

第七步，单击工具栏中的按钮“”，设置所画曲线的特性，如下图所示，并作出其解的三维图。

第八步，单击图2-8在标出的“plot”按钮，或单击工具栏中的按钮“”，可作出解的三维图，如下图所示。