牛顿法解非线性方程组实验报告

实验名称： 牛顿法解非线性方程组

1 引言

我们已经知道，线性方程组我们可以采取Jacobi迭代法，G-S迭代法以及SOR迭代方法求解。而在科学技术领域里常常提出求解非线性方程组的问题，例如，用非线性函数拟合实验数据问题、非线性网络问题，用差分法求解非线性微分方程问题等。

我们在解非线性方程组时，也考虑用迭代法求解，其思路和解非线性方程式一样，首先要将F(x)=0转化为等价的方程组

xigi(x1,x2,,xn),(i1,2,n)

或者简记为x=g（x），其中gi:RnR,g:RnRn

g1(x)x1g(x),xx2Rn2g(x) gn(x)xn迭代法：首先从某个初始向量x(0)开始，按下述逐次代入方法构造一向量序列{x(k)}：

(k1)(k)k)xigi(x1,,x(),(i1,2,,n) n(k)k)其中，x(k)(x1,x(,2k)T,x()。 n或写成向量形式：x(k1)g(x(k)),(k0,1,2,)

如果limx(k)x*（存在），称{x(k)}为收敛。且当gi(x)为连续函数时，可得

kx*g(limx(k))g(x*)

k说明x*为方程组的解。又称为x=g（x）的不动点。 本实验中采用牛顿迭代法来求解非线性方程组。

2 实验目的和要求

运用matlab编写一个.m文件，要求用牛顿法非线性方程组：

1xcosx2012T,(取x(0)(0,0),要求x(k1)x(k)x1sinx0122103)

3 算法原理与流程图

1、算法原理

设有非线性方程组 F(x)=0

其中：F(x)(f1(x),f2(x),T ,fk(x))由fi(x)偏导数作成的矩阵记为J(x)或F'(x)称为F(x)的Jacobi矩阵

f1(x)f1(x)x2x1f(x)f(x)22J(x)F'(x)x1x2fn(x)fn(x)x2x1(k)k)设x*为F(x)=0的解，且设x(k)(x1,x(,2f1(x)xnf2(x)xn

fn(x)xnk)T,x()，为x*的近似解，现利用多元函数nfi(x) 在x(k)点的泰勒公式有

fi(x)fi(x)(x1x12(k)(k)1fi(x(k)))x12(xnxn(k)fi(x(k)))xnjlnk)(k)(xjx()(xx)jll,1fi(Ci)Pi(x)Rxjxl

其中，Ci在x(k)与x的所连的线段内。

如果用泰勒公式中的线性函数Pi(x)近似代替fi(x)，并将线性方程组

(k)f(x)(k)Pi(x)fi(x(k))(x1x1)ix1(k)f(x)k)i(xnx()0n xn(i1,2,,n)的解作为x*的第k+1次近似解记为x(k1) 将上述方程写成矩阵形式：

F(x(k))J(x(k))(xx(k))0

如果J(x(k))为非奇异矩阵，则得到牛顿迭代公式：

x(0)(初始向量) (k1)(k)(k)1(k)x[J(x)]F(x),(k0,1,2,)x求解非线性方程组F(x)=0牛顿法或为

x(0)J(x(k))x(k)F(x(k)) x(k1)x(k)x(k)用上式可知，每计算一步x(k)x(k1)，需要： （1） 计算矩阵J(x(k))及F(x(k))；

（2） 求解一个线性方程组：J(x(k))x(k)F(x(k)) （3） 计算x(k1)x(k)x(k)。 2、流程图见附图1

4 程序代码及注释

%牛顿法解非线性方程组 function [Z,P,k,e] = newton(P,e0) %用P输入初始猜想矩阵，不断迭代输出计算解 %Z为迭代结束后的F矩阵 %k为迭代次数，e为每次迭代后的无穷范数，e0为误差限 Z=F(P(1),P(2)); J=JF(P(1),P(2)); Q=P-J\\Z; e=norm((Q-P),inf); P=Q; Z=F(P(1),P(2)); k=1; while e>=e0 J=JF(P(1),P(2)); Q=P-J\\Z; e=norm((Q-P),inf); P=Q; Z=F(P(1),P(2)); k=k+1; end end %子函数一，用来求每一步的F（x） function [out]=F(x,y) syms x1 x2; f1=x1-0.5*cos(x2); f2=x2-0.5*sin(x1); Y=[f1;f2]; x1=x; x2=y; out=subs(Y); end %子函数二，用来求每一步的Jacobi矩阵 function [y]=JF(x,y) syms x1 x2 f1=x1-0.5*cos(x2); f2=x2-0.5*sin(x1); df1x=diff(sym(f1),'x1'); df1y=diff(sym(f1),'x2'); df2x=diff(sym(f2),'x1'); df2y=diff(sym(f2),'x2'); j=[df1x,df1y;df2x,df2y]; %j中的元素为一阶偏导数 x1=x; x2=y; y=subs(j); end 5算例分析

（1）首先输入系数矩阵初始猜想解和误差限e0 >> P=[0 0]'; >> e0=0.001; （2）输出结果 >> [Z P k e]=newton(P,e0) Z = 1.0e-009 * 0.1366 0.7575 P = 0.5713 0.8820 k = 3 e = 5.0920e-005 其中，P为计算解，k为迭代次数，Z为第k次迭代后的F矩阵，e为第k次迭代后的 x(k1)x(k)

6讨论与结论

1、 时间复杂度： >> tic;[Z P k e]=newton(P,e0);toc Elapsed time is 0. seconds. 2、 程序优化 在本次的程序设计中，我采用了C语言中子函数调用的思想，使得程序的可读性增强，条理清晰。

其次，在求每一步的x(k)时，要解方程J(x(k))x(k)F(x(k))，此语句可以用

Q=P-inv（J）*Z;但考虑到算法的时间复杂性，本语句采用了Q=P-J\\Z来实现，从而减少了计

算机时间。

参考文献

[1] 易大义，沈云宝，李有法. 计算方法(第2版)，浙江大学出版社. p.29-53. [2] 张琨 高思超 毕靖 编著 MATLAB2010从入门到精通 电子工业出版社

输入n,x,ε,N0 k=1,2，…，N0 计算 (1)F(i)fi(x),(i1,2,(2)AJ(x),n) 求解线性方程组 J(x)Y=-F(x) x(k1)计算 xxy计算Sfi(x) 2Y S<ε 输出x，S，k kN0 N 输出迭代N0次还没有达到精度要求信息 Stop 附图1 流程图