matlab程序设计实践-牛顿法解非线性方程

1、精选优质文档-倾情为你奉上精选优质文档-倾情为你奉上专心-专注-专业专心-专注-专业精选优质文档-倾情为你奉上专心-专注-专业中南大学MATLAB程序设计实践学长有爱奉献，下载填上信息即可上交，没有下载券的自行百度。所需m文件照本文档做即可，即新建（FILE）脚本（NEW-Sscript）复制本文档代码运行（会跳出保存界面，文件名默认不要修改，保存）结果。第一题需要把数据文本文档和m文件放在一起。全部测试无误，放心使用。本文档针对做牛顿法求非线性函数题目的同学，当然第一题都一样，所有人都可以用。 记得删掉这段话班级： 学号：姓名： 一、MATLAB程序设计实践Matlab基础表示多晶体材料织构

2、的三维取向分布函数（ff（1，2）是一个非常复杂的函数，难以精确的用解析函数表达，通常采用离散空间函数值来表示取向分布函数，Data.txt是三维取向分布函数的一个实例。由于数据量非常大，不便于分析，需要借助图形来分析。请你编写一个matlab程序画出如下的几种图形来分析其取向分布特征：（1）用Slice函数给出其整体分布特征；（2）用pcolor或contour函数分别给出(20, 5, 10, 15, 20, 25, 30, 35 90)切面上f分布情况(需要用到subplot函数)；(3) 用plot函数给出沿取向线(1=090，45，20)的f分布情况。备注：表示以下数据为20的数据，

3、即f（1，0）此方向表示f随着从0,5,10,15, 20 到90的变化而变化此方向表示f随着1从0,5,10,15, 20 到90的变化而变化数据说明部分，与作图无关data.txt数据格式说明解：将文件Data.txt内的数据按照要求读取到矩阵f(phi1,phi,phi2)中，代码如下：fid=fopen(data.txt);for i=1:18 tline=fgetl(fid);endphi1=1;phi=1;phi2=1;line=0;f=zeros(19,19,19);while feof(fid) tline=fgetl(fid); data=str2num(tline); li

4、ne=line+1; if mod(line,20)=1 phi2=(data/5)+1; phi=1; else for phi1=1:19 f(phi1,phi,phi2)=data(phi1); end phi=phi+1; endendfclose(fid);将以上代码保存为readdata.m在MATLAB中运行，运行结果如下图所示：将以下代码保存为code1.m文件：fopen(readdata.m);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)运行结果如下图所示：将以下代

5、码保存为code2.m文件：fopen(readdata.m);readdata;for i=1:19 subplot(5,4,i) pcolor(f(:,:,i)nd运行结果如下图所示：将以下代码保存为code3.m文件：fopen(readdata.m);readdata;for i=1:19 subplot(5,4,i) contour(f(:,:,i)end运行结果如下图所示：1=090，45，20所对应的f(1,2)即为f(:,10,1)。将以下代码保存为code4.m文件：fopen(readdata.m);readdata;plot(0:5:90,f(:,10,1),-bo)te

6、xt(60,6,phi=45 phi2=0)运行结果如下图所示：二 MATLAB程序设计实践科学计算（24）班级：学号：姓名： 、编程实现以下科学计算算法，并举一例应用之。（参考书籍精通科学计算，王正林等著，电子工业出版社，年）“牛顿法非线性方程求解”解：弦截法本质是一种割线法，它从两端向中间逐渐逼近方程的根；牛顿法本质上是一种切线法，它从一端向一个方向逼近方程的根，其递推公式为：初始值可以取和的较大者，这样可以加快收敛速度。和牛顿法有关的还有简化牛顿法和牛顿下山法。在MATLAB中编程实现的牛顿法的函数为：NewtonRoot。功能：用牛顿法求函数在某个区间上的一个零点。调用格式：root=

7、NewtonRoot其中，为函数名； 为区间左端点； 为区间右端点 eps为根的精度； root为求出的函数零点。，NY输入参数值nargin=3f1=0eps=1.0e-4结束f2=0f1*f20两端点函数值乘积大于0 ！YNdfadfbYNNNYY开始tolepsY输出结果N流程图：root=aroot=broot=a-fa/dfaroot=b-fb/dfbroot=r1-fx/dfx牛顿法的matlab程序代码如下：function root=NewtonRoot(f,a,b,eps) %牛顿法求函数f在区间a,b上的一个零点%函数名:f%区间左端点:a%区间右端点:b%根的精度:eps

8、%求出的函数零点:rootif(nargin=3)eps=1.0e-4;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;end if (f1*f20) disp(两端点函数值乘积大于0 ！); return;else tol=1; fun=diff(sym(f); %求导数 fa=subs(sym(f),findsym(sym(f),a); fb=subs(sym(f),findsym(sym(f),b); dfa=subs(sym(f

9、un),findsym(sym(fun),a); dfb=subs(sym(fun),findsym(sym(fun),b); if(dfadfb) %初始值取两端点导数较大者 root=a-fa/dfa; else root=b-fb/dfb; end while(toleps) 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 例：求方程3x2exp（x）0

10、的一根解：在MATLAB命令窗口输入： r=NewtonRoot(3*x2-exp(x),3,4)输出结果：X=3.7331流程图：开始 调用M文件输入方程程序运行输出结果结束、编程解决以下科学计算问题。1)解：这个方程可用下列步骤来解 （1）用eig函数求出矩阵K- M的特征值L和特征向量U，U和L满足 （2）在原始方程Mx+Kx=0两端各乘以及在中间乘以对角矩阵*U，得 *M*U*+*K*x=0 作变量置换，得 这是一个对角矩阵方程，即可把它分两个方程： 这意味着两种振动模态可以解耦，令，是第i个模式的固有频率（i=1,2）。 (3)由上述的解耦模态中，给出初始条件，化为，即可求出其分量，

11、。 设位置和速度的初始条件分别为，则这三个步 骤得到的最后公式为输入参数值结束开始流程图：构成参数矩阵设定初始条件 i=1:round(tf/dt)+1循环计算矩阵指数绘制图象系统解耦的振动模态的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=0.1; %步长%构成二阶参数矩阵M=m1,0

12、;0,m2;K=k1+k2,-k2;-k2,k2;C=c1+c2,-c2;-c2,c2;%构成四阶参数矩阵A=zeros(2,2),eye(2);-MK,-MC;%四元变量的初始条件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),ylabe

13、l(y);运行M文件，依下图所示在MATLAB命令窗口中输入数据：即可得出该振动的两种模态2) 解：第一步，在MATLAB命令窗口输入命令pdetool打开工具箱，调整x坐标范围为0 4,y坐标范围为0 3.通过options选项的Axes Linits设定如下图所示。第二步，设定矩形区域。点击工具箱栏中的按钮“”，拖动鼠标画出一矩形，并双击该矩形，设定矩形大小，如下图所示。第三步，设边界条件。点击工具栏中的按钮“”，并双击矩形区域的相应的边线在弹出的对话框中设定边界条件。如下图所示，分别为各边框的边界条件。第四步，设定方程。单击工具栏中的按钮“”，在PDE模式下选择方程类型，如下图所示，并在其中设定参数。