牛顿法非线性方程求解

《MATLAB程序设计实践》课程考核---第39-40页题１：编程实现以下科学计算算法，并举一例应用之。（参考书籍《精通MATLAB科学计算》，王正林等著，电子工业出版社，2009年）“牛顿法非线性方程求解”弦截法本质是一种割线法，它从两端向中间逐渐逼近方程的根；牛顿法本质上是一种切线法，它从一端向一个方向逼近方程的根，其递推公式为：nnxx1)()('nnxfxf初始值可以取)('af和)('bf的较大者，这样可以加快收敛速度。和牛顿法有关的还有简化牛顿法和牛顿下山法。在MATLAB中编程实现的牛顿法的函数为：NewtonRoot。功能：用牛顿法求函数在某个区间上的一个零点。调用格式：root=NewtonRoot)(```epsbaf其中，f为函数名；a为区间左端点；b为区间右端点eps为根的精度；root为求出的函数零点。，牛顿法的matlab程序代码如下：functionroot=NewtonRoot(f,a,b,eps)%牛顿法求函数f在区间[a,b]上的一个零点%函数名:f%区间左端点:aNY输入参数值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%区间右端点:b%根的精度:eps%求出的函数零点: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;endif(f1*f20)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(dfadfb)%初始值取两端点导数较大者root=a-fa/dfa;elseroot=b-fb/dfb;endwhile(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例：求方程3x^2－exp（x）＝0的一根解：在MATLAB命令窗口输入：r=NewtonRoot('3*x^2-exp(x)',3,4)输出结果：X=3.73312、编程解决以下科学计算问题1）二自由度可解耦系统的振动模态分析，二自由度振动系统如图所示，其一般方程为：0)12(2)12(2220221)21(221)21(11........xxKxxcxmxKxKKxcxccxm可写成矩阵形式：0...KxxCxM设C=0，即无阻尼情况，则系统可解耦为两种独立的振动模态。输入参数值结束开始流程图：构成参数矩阵设定初始条件i=1:round(tf/dt)+1循环计算矩阵指数绘制图象系统解耦的振动模态的MATLAB代码如下：functionerziyoudu()%输入各原始参数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;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];%设定计算点，作循环计算fori=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)用GUI方式解下列PDE：;0,4sin0),3(30,40,030402222yyxxuxuuyyuyxyuxu解：第一步，在MATLAB命令窗口输入命令pdetool打开工具箱，调整x坐标范围为[05],y坐标范围为[05].通过options选项的AxesLinits设定如下图所示。第二步，设定矩形区域。点击工具箱栏中的按钮“”，拖动鼠标画出一矩形，并双击该矩形，设定矩形大小，如下图所示。第三步，设边界条件。点击工具栏中的按钮“”，并双击矩形区域的相应的边线在弹出的对话框中设定边界条件。如下图所示，分别为各边框的边界条件。第四步，设定方程。单击工具栏中的按钮“”，在PDE模式下选择方程类型，如下图所示，并在其中设定参数。第五步，单击工具栏中的按钮“”，拆分区域为若干子区域，如下图所示。第六步，单击工具栏中的按钮“”，将子区域细化，从而保证结果更精确，如下图所示。第七步，单击工具栏中的按钮“”，设置所画曲线的特性，如下图所示，并作出其解的三维图。第八步，单击图2-8在标出的“plot”按钮，或单击工具栏中的按钮“”，可作出解的三维图，如下图所示。简要流程图：开始绘制要求区域图设置边界条件设置方程参数剖分网格作出解的三维图结束