牛顿法非线性方程求解.docVIP

《MATLAB 程序设计实践》课程考核 ---第39-40页 题１ ： 编程实现以下科学计算算法，并举一例应用之。（参考书籍《精通MATLAB科学计算》，王正林等著，电子工业出版社，2009 年） “牛顿法非线性方程求解” 弦截法本质是一种割线法，它从两端向中间逐渐逼近方程的根；牛顿法本质上是一种切线法，它从一端向一个方向逼近方程的根，其递推公式为： 初始值可以取和的较大者，这样可以加快收敛速度。 和牛顿法有关的还有简化牛顿法和牛顿下山法。 在MATLAB中编程实现的牛顿法的函数为：NewtonRoot。 功能：用牛顿法求函数在某个区间上的一个零点。 调用格式：root=NewtonRoot 其中，为函数名； 为区间左端点； 为区间右端点 eps为根的精度； root为求出的函数零点。 ， 牛顿法的matlab程序代码如下： function root=NewtonRoot(f,a,b,eps) %牛顿法求函数f在区间[a,b]上的一个零点 %函数名:f %区间左端点:a %区间右端点:b %根的精度:eps %求出的函数零点:root if(nargin==3) eps=1.0e-4; end f1=subs(sym(f),findsym(sym(f)),a); f2=subs(sym(f),findsym(sym(f)),b); if (f1==0) root=a; end if (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(fun),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); end end 例：求方程3x^2－exp（x）＝0的一根 解：在MATLAB命令窗口输入： r=NewtonRoot(3*x^2-exp(x),3,4) 输出结果： X=3.7331 2、编程解决以下科学计算问题 二自由度可解耦系统的振动模态分析，二自由度振动系统如图所示，其一般方程为： 可写成矩阵形式： 设C=0，即无阻尼情况，则系统可解耦为两种独立的振动模态。 系统解耦的振动模态的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;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)+1 t(i)=dt*(i-1); y(:,i)=expm(A*t(i))*y0;%循环计算矩阵指数 end %按两个分图绘制x1、x2曲线 subplot(2,1,1),plot(t,y(1,:)),grid xlabel(t),ylabel(y); subplot(2,1,2),plot(t,y(2,:)),grid xlabel(t),ylabel(y); 运行M