牛顿迭代法求解非线性方程组的解

数值计算实习报告——非线性方程组的牛顿迭代解法姓名：吴健学号：139084186班级：数132非线性方程组的数值解法摘要本文着重介绍了求解单变量非线性方程()0fx的牛顿迭代法理论和MATLAB求解过程，其中二分法、不动点迭代等也是常用的解非线性方程组的重要工具，本文把牛顿迭代法推广到非线性方程组，对于方程变量及个数相当大时，很难通过运算得出结果，此时采用MATLAB中的numjac命令很好的避免牛顿迭代法中所遇到的jacobi矩阵难求的问题。关键词：非线性方程组、牛顿迭代法、MATLAB、jacobi矩阵一、前言非线性方程组在实际问题中经常出现，并且在科学与工程计算中的地位越来越来重要，很多常见的线性模型都是在一定条件下由非线性问题简化得到的，为得到更符合实际的解答，往往需要直接研究非线性模型，然而从线性到非线性是一个质的飞跃，方程的性质的不同，所以求解方法也有很大差别。本文主要介绍关于非线性方程及方程组的数值解法，先分析非线性方程的数值解法，然后再延伸到方程组的解法。本文以如下多元非线性方程组为例，尝试寻找一种数值方法求解该方程组。1212322123313cos()0281(0.1)sin1.0601020103xxxxxxxxex二、线性方程与非线性方程线性方程的形式如ykxb(2-1)而对于这类问题已经完全地被解决了，而诸如函数()fx是多项式函数，即1011()nnnnfxaxaxaxa(2-2)其中00,(0,1,,)iaain为实数，当()0fx称为n次代数方程。这类方程当5n时就不能直接用公式表示方程的根，所以只好选用数值解代替。还有一类函数成为超越函数，如/10sin100xex(2-3)它在整个x轴上有无穷多个解，若x的取值范围不同，解也不同，所以在讨论非线性方程的求解问题时，必须强调x的定义域，即x的求解区间[,]ab。超越函数都可以通过多项式函数插值得到，故解非线性方程归结为求解多项式函数的解或者是直接求解超越函数。三、求解非线性方程的其他方法1、二分法考察有根区间[,]ab，取中点0()/2xab将它分成两半，假设中点0x不是()fx的零点，然后进行根的搜索，即检查0()fx与()fa是否同号，如果确系同号，则说明所求根*x在0x的右侧，这时令101,axbb;否则*x必在0x的左侧，这时令110,aabx.对于压缩了的有根区间11[,]ab进行如上用样的操作，直到找到方程的根。2、不动点迭代法将方程()0fx(3-1)改写成等价的形式()xx,若要求*x满足*()0fx,则**()xx(3-2)反之亦然.称*x为函数()x的一个不动点.求()fx的零点就等价于求()x的不动点，选择一个初始值0x，将它代入()xx的右端，即可求得10()xx，如此反复迭代计算1(),0,1,kkxxk(3-3)这时的()x称为迭代函数.如果对0[,]xab,有*limkkxx(3-4)成立，则称迭代方程1()kkxx收敛,且(3-2)为()x的不动点，故称为不动点迭代法。其基本思想是将隐式方程()xx归结为一组显示的计算公式(3-3)，就是说，迭代过程实质上是一个逐步显示化的过程。3、线性近似法线性近似法中比较突出的有牛顿法、弦截法等，牛顿法的基本思想是将非线性方程()0fx逐步归结为某种线性方程来求解。设kx,1kx是()0fx的近似根，我们利用()kfx,1()kfx构造一次插值多项式1()px，并用1()0px的根作为()0fx的新的近似根1kx.由于111()()()()()kkkkkkfxfxpxfxxxxx(3-5)因此有111()()()()kkkkkkkfxxxxxfxfx(3-6)由下文中的牛顿迭代法会知道，这里的迭代公式可以看作牛顿公式1'()()kkkkfxxxfx(3-7)中的导数'()kfx用差商11()()kkkkfxfxxx取代的结果，这种算法成为弦截法。4、上述算法的缺点二分法：收敛速度慢，且无法知道是否有重根。不动点迭代法：收敛与否或收敛快慢取决于函数()xx的性质。弦截法：在牛顿法中只需要两个节点，弦截法需要找三个节点，收敛速度慢，要求初始值的选取靠近方程的根，否则也可能会不可能。四、牛顿法1、牛顿法思想设已知方程()0fx有近似根kx（假定'()0kfx），将函数()fx在点kx展开，有'()()()()kkkfxfxfxxx(4-1)于是方程()0fx可以近似地表示为'()()()0kkkfxfxxx(4-2)式（4-2）是线性方程，记其根为1kx，则1kx的计算公式为1'()0,1,()kkkkfxxxkfx(4-3)这就是牛顿法，亦称为切线法。牛顿法的几何意义如下图所示，方程()0fx的根*x可解释为曲线()yfx与x轴的交点的横坐标。设kx是根*x的某个近似值，过曲线()yfx上横坐标为kx的点kP引切线，并将该切线与x轴的交点的横坐标1kx作为*x的新的近似值。注意到切线方程为'()()()kkkyfxfxxx(4-4)这样求得的值1kx比满足'()()()0kkkfxfxxx，从而就是牛顿公式1'()()kkkkfxxxfx(4-5)的计算结果。由几何意义，牛顿法亦称为切线法。图4-12、牛顿法的算法步骤步骤一：准备，选定初始近似值0x,计算00()ffx,''00()ffx；令k=0步骤二：迭代，按公式'1000/xxff迭代一次，得新的近似值1x，计算11()ffx，''11()ffx；再分别赋值给01ff，''01ff。步骤三：判断，如果10||15xxe，则终止迭代，以1x作为所求的根；否则转步骤四步骤四：修改，如果迭代次数达到预先指定的次数N，或者'10f，则方法失败；否则以'111(,,)xff代替'000(,,)xff转步骤二继续迭代。例：求方程3()10fxxx在01.5x附近的根.解：结果如下表所示：迭代次数0123456初始值1.50001.25001.35801.31041.33101.32201.3258迭代值1.25001.35801.31041.33101.32201.32591.3242误差值0.25000.10800.04760.02060.00890.00390.00173、牛顿法的改进牛顿法的缺点：一是每次迭代都要计算()kfx和'()kfx，计算量较大且有时'()kfx计算较困难；二是初始近似0x只在根*x附近才能保证收敛，如果0x给的不合适可能不收敛。为克服这两个缺点，通常有下述方法。01简化牛顿法简化牛顿法又称平行弦法，其迭代公式为1(),0,0,1,kkkxxCfxCk(4-7)从不动点迭代法的角度看，简化牛顿法的迭代函数()()xxCfx，下面讨论简化牛顿法的收敛性。若''|()||1()|1xCfx,即取'0()2Cfx.在根*x附近成立，则迭代法1()kkkxxCfx局部收敛，在1()kkkxxCfx中取'01()Cfx，则称为简化牛顿法。其几何意义是用斜率为'0()fx的平行弦与x轴的交点作为*x的近似，如图所示：图4-22o牛顿下山法牛顿法收敛性依赖于初值0x的选取，如果0x偏离所求根*x较远，则牛顿法可能发散。为了防止迭代发散，对迭代过程再附加一个要求，即具有单调性：1|()||()|kkfxfx(4-8)满足此要求的算法称为下山法。将牛顿法和下山法一起使用时，即在下山法保证函数值稳定下降的前提下，用牛顿法加快收敛速度。为此，为此将牛顿法的计算结果1'()()kkkkfxxxfx(4-9)与前一步的近似值kx的适当加权平均作为新的改进值11(1)kkkxxx(4-10)其中(01)称为下山因子，11(1)kkkxxx即为1'(),0,1,()kkkkfxxxkfx(4-11)称为牛顿下山法。选择下山因子时从1开始，逐次将减半进行试算，直到能使下降条件(4-8)成立。五、算法实施例1、求方程3()10fxxx在01.5x附近的根.解：（1）二分法(1)10f，(2)50f，故取区间[1,2]，则由二分法的算法得MATLAB程序：函数值[1,2][1,1.5][1.25,1.5][1.25,1.375][1.3125,1.375][1.3125,1.3438][1.3125,1.3281]()fa-1-1-0.2969-0.2969-0.0515-0.0515-0.0515()fb50.87500.87500.22460.22460.08260.0147误差值10.50.250.1250.06250.031250.015625（2）不动点迭代法：根据不动点迭代法的思想，取迭代函数311kkxx，则有结果如下表：迭代次数1234567迭代函数(x+1)^(1/3)1.35721.33091.32591.32491.32481.32471.3247误差值0.14280.02630.00500.00100.00010.00010.0000六、推广牛顿法到非线性方程组求解考虑方程组11221212(,,,)0(,,,)0(,,,)0nnnnfxxxfxxxfxxx(6-1)其中12,,,nfff均为12(,,,)nxxx的多元函数.若向量记号12(,,,)TnnXxxx，12(,,,)TnFfff，方程组就可以写成()0FX(6-2)当2n，且(1,2,,)ifin中至少有一个是自变量(1,2,,)ixin的非线性函数时，则称方程组(6-1)为非线性方程组.非线性方程组是非线性科学的重要组成部分，关于非线性方程组的解法至关重要的，其基本思想是：将单个方程的牛顿法直接用于方程组(6-2)则可得到解非线性方程组的牛顿迭代法(1)()'()1()()(),0,1,kkkkXXFXFXk(6-3)其中，'1()FX是非线性方程组的Jacobi矩阵的逆矩阵，具体计算时记(1)()()kkkXXX先解线性方程组'()()()()()kkkFXXFX，求出向量()kX，再令(1)()()kkkXXX.注：Jacobi矩阵为11112222'1212()()()()()()()()()()nnnnnnfxfxfxxxxfxfxfxxxxFxfxfxfxxxx(6-4)七、非线性方程组牛顿法实例例2：应用牛顿迭代法解论文前言中给出的非线性方程组即(1-1)式，1212322123313cos()0281(0.1)sin1.0601020103xxxxxxxxex解：由MATLAB实现程序结果如下：初始值X0=[0;0;0]X1=[1;1;1]X2=[-5;-6;-7]迭代结果X=[0.5;0;-0.5236]X=[0.5;0;-0.5236]无结果迭代次数68100注：对于不同的初值可能得到的结果会不一样，这里不再一一说明，因为前文已经论述了为什么会出现这种情况。八、大型非线性方程组求解非线性方程组所含未知元个数和方程个数较大时，求方程组的jacobi矩阵会计算繁琐，效率低，考虑借助于MATLAB去寻找jacobi矩阵的数值解，通过MATLAB已有的numjac命令去处理。程序实现：clear;clc;formatlongg;F=@(t,x)[3*x(1)-cos(x(2)*x(3))-1/2;x(1)^2-81*(x(2)+0.1)^2+sin(x(3))+1.06;exp(-x(1)*x(2))+20*x(3