关于非线性方程组Newton 解法的研究综述
作者:谢玉婧 王培
一、研究现状
非线性代数方程组求解是一个基本而又重要的问题,这是因为在工程实践、经济学、信息安全和动力学等方面有大量的实际问题最终转化为代数方程组。这一类问题,我们不可能到它们的解析解,数值解是目前主要的研究方向。数值解法较为成熟,速度快,但其往往只能求出部分解,而且通常只能求出近似解。它一般用于解决大型问题。Newton 法是数值方法求解非线性方程组的一种基本方法。根据其特点和不足,对Newton 法进行简化和修正。并且与其他数值解法(如区间法)相结合,形成其他更加完善的求解非线性方程组的方法。
二、问题的提出
n 个变量个方程的非线性方程组, 其一般形式如下:
n 112n 212n n 12n f (x ,x ,,x )0f (x ,x ,,x )0f (x ,x ,,x )0
=⎧⎪=⎪⎨⎪⎪=⎩L L M L )1(式(1)中是定义在维欧式空间中开域上的实值函数。若用向量记,令:
)n ,,2,1i )(x x ,x (f n 21i L L =n n R D ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡====⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎣⎡=X)(f X)(f )X (f 0)x ,x ,(x f 0)x ,x ,(x f 0)x ,x ,(x f )X (F ,x x x X n 21n 21n n 212n 211n 21M L M L L M 则方程组(1)也可以表示为:
0)X (F = )2(其中:,n 0n n R )X (F ,R R :F ,R X ∈→∈n R 为赋值空间。
非线性数值方法种类较多,最常用的是Newton 迭代法。求解非线性方程组的Newton 法是一个最基本而且十分重要的方法,目前使用的很多有效的迭代法
都是以Newton 法为基础,或由它派生而来。下面首先介绍几种Newton 法及其原理,然后使用Mathwmatica 编写求解程序,最后通过一个计算实例说明几种牛顿法的差异。为方便起见
记:
T n 1T n 1]x ,x [X ,)]X (f ,),X (f [)X (F L L ==以及雅可比矩阵:
⎥⎥⎥⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂=′n n 2n 1n n 22212n 12111x )X (f x )X (f x )
X (f x )X (f x )X (f x )
X (f x )X (f x )X (f x )X (f )X (F L M L L 三、各种方法原理及math程序
1、Newton 法
(1) 方法原理简述
像单个方程的Newton 迭代法一样,采用逐次线性化的方法构造方程(1)的Newton 迭代法。在某个近似解处,将向量函数F(作Taylor 展开,则有:
k X X)()1F()F()F ()k k k k X X X X X +′≈+−
从而得到(2)的近似方程
()1F()F ()0k k k k X X X X +′+−=
令,将Newton 法的迭代公式改写为:
K 1k k X X X −=Δ+ ⎩⎨⎧−=Δ′Δ+=+)
X (F )X )(X (F X X X k k k k k 1k L ,2,1,0k = (3) 每一步迭代均需解Newdon 方程组:
)X (F )X )(X (F k k k −=Δ′
这是一个线性方程组,可用高斯消去法求解。 通常可用:
ε<Δk X 或δ<)X (F k
作为Newton 法的终止迭代准则,其中δ,ε为预先给定的精度要求,有时也可用
预先确定的最大迭代次数M 作为终止迭代条件。
(2) math 程序及程序说明
Clear[x,y,z]
f[1][x_,y_,z_]:=4x+y^2+Exp[-2*z]-8.03;
f[2][x_,y_,z_]:=Sin[x]-y*(z+10)+3.01;
f[3][x_,y_,z_]:=-2(x+0.3)^2-Cos[y]+10z+3Pi;
F1[x_,y_,z_]:=Table[{D[f[i][x,y,z],x],D[f[i][x,y,z],y],D[f[i][x,y,z],z]},{i,1,3}]; F1[x,y,z];
MatrixForm[%]
I1=IdentityMatrix[3];
DD=DiagonalMatrix[F1[x,y,z][[1,1]],F1[x,y,z][[2,2]],F1[x,y,z][[3,3]]]//N;
L=-{{0,0,0},{F1[x,y,z][[2,1]],0,0},{F1[x,y,z][[3,1]],F1[x,y,z][[3,2]],0}}//N;
DDLN=Inverse[DD-L]//N;
G=I1-DDLN.F1[x,y,z];
MatrixForm[%];
F[x_,y_,z_]:=Table[f[i][x,y,z],{i,1,3}];
f=DDLN.(-F[x,y,z]//N);
{x,y,z}={2,2,2};
For[k=1,k<=10,k++,
{tx,ty,tz}={1,1,1};
Do[{tx1,ty1,tz1}=G .{tx,ty,tz}+f//N;
{tx,ty,tz}={tx1,ty1,tz1},{i,1,10}];
{x,y,z}={x+tx1,y+ty1,z+tz1};
h=Max[Abs[tx1],Abs[ty1],Abs[tz1]]//N;
Print[k," ",x," ",y," ",z," ",h]]
程序2-6行是在定义雅克比矩阵(方程组的维数及各个方程的具体形式以后面出现的例题为参考),7-14行是高斯法求解方程组的矩阵分解准备,15-22行是程序的主体部分,是一个循环,是对方程组(3)的数学程序语言翻译。在这个循环中又嵌套了一个循环(18-19行),是高斯法求解)X (F )X )(X (F k k k −=Δ′的具体形式,求出k X Δ。其中h 是误差控制(ε<Δk X )。
2、简化的Newdon 法
(1) 方法原理简述
尽管Newton 法具有较高的收敛速度,但在每一步迭代中,要计算个函数值以及个导数值f 形成雅可比矩阵n f 2n ′k F (X )′而且要求k F (X )′的逆矩阵或者解一
个n 阶线性方程组,计算量很大。为了减少计算量。在上述Newton 法中,对于一切,取为,于是迭代公式便修正成为:
L ,2,1,0k =k k F (X )′0F (X )′ k 1k k 0k X X X F (X )(X )F(X )
+=+Δ⎧⎨′Δ=−⎩L ,2,1,0k = (4) 式(3)即为简化的Newton 法。该方法能使计算量大为减少,但却大大降低了收敛速度。简化的Newton 法的算法与Newton 法的算法区别就在于只由给定的初始近似值计算一次F (,以后在每一次迭代过程中不再计算F (X)′X)′, 保持初始计算值。
(2) math 程序及程序说明
Clear[x,y,z]
f[1][x_,y_,z_]:=4x+y^2+Exp[-2*z]-8.03;
f[2][x_,y_,z_]:=Sin[x]-y*(z+10)+3.01;
f[3][x_,y_,z_]:=-2(x+0.3)^2-Cos[y]+10z+3Pi;
F1[x_,y_,z_]:=Table[{D[f[i][x,y,z],x],D[f[i][x,y,z],y],D[f[i][x,y,z],z]},{i,1,3}]; F1[x,y,z];
N[%/.x->2,6];
N[%/.y->2,6];
F0=N[%/.z->2,6];
MatrixForm[%]
I1=IdentityMatrix[3];
DD=DiagonalMatrix[F0[[1,1]],F0[[2,2]],F0[[3,3]]]//N;
L=-{{0,0,0},{F0[[2,1]],0,0},{F0[[3,1]],F0[[3,2]],0}}//N;
DDLN=Inverse[DD-L]//N;
G=I1-DDLN.F0;
MatrixForm[%];
F[x_,y_,z_]:=Table[f[i][x,y,z],{i,1,3}];
f=DDLN.(-F[x,y,z]//N);
{x,y,z}={2,2,2};
{tx,ty,tz}={1,1,1};
Do[{tx1,ty1,tz1}=G .{tx,ty,tz}+f//N;
{tx,ty,tz}={tx1,ty1,tz1},{k,1,10}];
{x,y,z}={x+tx1,y+ty1,z+tz1};
h=Max[Abs[tx1],Abs[ty1],Abs[tz1]]//N;
Print[k," ",x," ",y," ",z," ",h]]
与上面Newton 法类似,区别在于整个工程中只计算了一次雅克比矩阵,2-10行是在计算。
0F (X )′3、修正的Newdon 法
(1) 方法原理简述
从以上两种方法可以看出, 它们各有利弊,如果吸取法收敛快与简化的
Newton 法工作量省的优点, 把m 步简化的Newton 步合并成一次Newton 步。则可以得到下列迭代程序:
(5)
⎪⎭
⎪⎬⎫=′−==+−−−m ,k 1k 1j ,k 1k 1j ,k j ,k k
0,k X X )f(X )X (f X X X X 式中:,,该式称为修正的Newton 法。在实际应用中,较为常用的是m ,,2,1j L =L ,2,1,0k =2m =的情形。此时,(5)式也可以简化为: L ,2,1,0k )]],X (f )]X (f -X [f )X ([f )]X (f [x x k 1k k k 1k k 1k =′+−=−−+ (6) 下面给出的程序是m 取一般值时的情形。
(2) math 程序及程序说明
Clear[x,y,z]
{x0,y0,z0}={2,2,2}//N;
f[1][x_,y_,z_]:=4x+y^2+Exp[-2*z]-8.03;
f[2][x_,y_,z_]:=Sin[x]-y*(z+10)+3.01;
f[3][x_,y_,z_]:=-2(x+0.3)^2-Cos[y]+10z+3Pi;
F1[x_,y_,z_]:=Table[{D[f[i][x,y,z],x],D[f[i][x,y,z],y],D[f[i][x,y,z],z]},{i,1,3}]; F1[x,y,z];
Clear[F1];
N[F1[x,y,z]/.x->x0,8];
N[%/.y->y0,8];
F0=N[%/.z->z0,8];
MatrixForm[%];
I1=IdentityMatrix[3];
DD=DiagonalMatrix[F0[[1,1]],F0[[2,2]],F0[[3,3]]]//N;
L=-{{0,0,0},{F0[[2,1]],0,0},{F0[[3,1]],F0[[3,2]],0}}//N;
DDLN=Inverse[DD-L]//N;
identity matrix是什么意思G=I1-DDLN.F0;
MatrixForm[%];
F[x_,y_,z_]:=Table[f[i][x,y,z],{i,1,3}];
f=DDLN.(-F[x,y,z]//N);
{x,y,z}={2,2,2};
For[k=1,k<=10,k++,
{tx,ty,tz}={1,1,1};
Do[{tx1,ty1,tz1}=G.{tx,ty,tz}+f//N;
{tx,ty,tz}={tx1,ty1,tz1},{k,1,10}];
{x,y,z}={x+tx1,y+ty1,z+tz1};
h=Max[Abs[tx1],Abs[ty1],Abs[tz1]]//N];
{x0,y0,z0}={x,y,z};
Print[x," ",y," ",z," ",h]]
该程序是随前两种程序的结合和改进。Newton法收敛速度较快,但是计算量大,计算时间较长。其主要原因是每次高斯法求解方程组都要计算雅克比矩阵,而简单Newton法避免了这个问题。这也是两者可以结合的主要契合点。程序中m是Newton步中简单步的子步数。即通过m控制雅克比矩阵计算次数的(最外层For循环)。整个程序相当于运行了m个简单Newton步。每一个简单步的结
X
果作为下一次运行的。
四、应用实例
版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198,我们将在24小时内删除。
发表评论