牛顿法非线性方程求解

《MATLAB 程序设计实践》课程考核

---第37-38页

题１ ： 编程实现以下科学计算算法，并举一例应用之。（参考书籍《精

通ＭＡT ＬＡＢ科学计算》，王正林等著，电子工业出版社，２００９ 年）

“牛顿法非线性方程求解” 弦截法本质是一种割线法，它从两端向中间逐渐逼近方程的根；牛顿法本质上是一种切线法，它从一端向一个方向逼近方程的根，其递推公式为：

-

=+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

%求出的函数零点: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*f2>0)

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(dfa>dfb) %初始值取两端点导数较大者

root=a-fa/dfa;

else

root=b-fb/dfb;

end

while(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);

end

end

例：求方程3x^2－exp（x）＝0的一根

解：在MATLAB命令窗口输入：

>> r=NewtonRoot('3*x^2-exp(x)',3,4)

输出结果：

X=3.7331

2、编程解决以下科学计算问题

1） 二自由度可解耦系统的振动模态分析，二自由度振动系统如图所示，其一般方程为：

)12(2)12(2220

221)21(221)21(11.

.

..

.

.

..

=-+-+=-++-++x x K x x c x m x K x K K x c x c c x m

可写成矩阵形式：0.

..

=++Kx x C x M

设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文件，依下图所示在MATLAB命令窗口中输入数据：

即可得出该振动的两种模态