Newton迭代法求解非线性方程

Newton迭代法求解非线性方程一、Newton迭代法概述构造迭代函数的一条重要途径是用近似方程来代替原方程去求根。因此，如果能将非线性方程f(x)=0用线性方程去代替，那么，求近似根问题就很容易解决，而且十分方便。牛顿(Newton)法就是一种将非线性方程线化的一种方法。设kx是方程f(x)=0的一个近似根，把如果)(xf在kx处作一阶Taylor展开，即：)xx)(x('f)x(f)x(fkkk(1-1)于是我们得到如下近似方程：0)xx)(x('f)x(fkkk(1-2)设0)('kxf，则方程的解为：x̅=xk+f(xk)f(xk)́(1-3)取x~作为原方程(1.1)的新近似根1kx，即令：)x('f)x(fxxkkk1k,k=0,1,2,…(1-4)上式称为牛顿迭代格式。用牛顿迭代格式求方程的根的方法就称为牛顿迭代法，简称牛顿法。牛顿法具有明显的几何意义。方程：)xx)(x('f)x(fykkk(1-5)是曲线)x(fy上点))x(f,x(kk处的切线方程。迭代格式(1-4)就是用切线式(1-5)的零点来代替曲线的零点。正因为如此，牛顿法也称为切线法。牛顿迭代法对单根至少是二阶局部收敛的，而对于重根是一阶局部收敛的。一般来说，牛顿法对初值0x的要求较高，初值足够靠近*x时才能保证收敛。若要保证初值在较大范围内收敛，则需对)x(f加一些条件。如果所加的条件不满足，而导致牛顿法不收敛时，则需对牛顿法作一些改时，即可以采用下面的迭代格式：)x('f)x(fxxkkk1k,,2,1,0k(1-6)上式中，10，称为下山因子。因此，用这种方法求方程的根，也称为牛顿下山法。牛顿法对单根收敛速度快，但每迭代一次，除需计算)x(fk之外，还要计算)x('fk的值。如果)x(f比较复杂，计算)x('fk的工作量就可能比较大。为了避免计算导数值，我们可用差商来代替导数。通常用如下几种方法：1.割线法如果用1kk1kkxx)x(f)x(f代替)x('fk，则得到割线法的迭代格式为：)x(f)x(f)x(fxxxxk1kk1kkk1k(1-7)2.拟牛顿法如果用)x(f))x(fx(f)x(fk1kkk代替)x('fk，则得到拟牛顿法的迭代格式为：))x(fx(f)x(f)x(fxx1kkkk2k1k(1-8)3.Steffenson法如果用)x(f)x(f))x(fx(fkkkk代替)x('fk，则得到拟牛顿法的迭代格式为：)x(f))x(fx(f)x(fxxkkkk2k1k(1-9)二、算法分析1.割线法割线法的迭代公式为：)x(f)x(f)x(fxxxxk1kk1kkk1k,k=0,1,2,…割线法是超线性收敛，其程序流程图为：2.拟牛顿法牛顿拟迭代法迭代公式为：))x(fx(f)x(f)x(fxx1kkkk2k1k(1)对单根条件下，牛顿拟迭代法平方收敛，牛顿拟迭代法程序框图如下所示：(2)对重根条件下，此时迭代公式修改为：))x(fx(f)x(f)x(fmxx1kkkk2k1k,m为根的重数此时，牛顿迭代法至少平方收敛。3.Steffenson法Steffenson迭代法程序流程图与牛顿拟迭代法类似。三、牛顿法的程序给定初值0p，用牛顿法格式)p('f)p(fpp1k1k1kk，,2,1k，求解非线性方程0)x(f。*********************************************************************function[p1,err,k,y]=newton(f1041,df1041,p0,delta,max1)%f1041是非线性函数。%df1041是f1041的微商。%p0是初始值。%delta是给定允许误差。%max1是迭代的最大次数。%p1是牛顿法求得的方程的近似解。%err是p0的误差估计。%k是迭代次数。%y=f(p1)p0,feval('f1041',p0)fork=1:max1p1=p0-feval('f1041',p0)/feval('df1041',p0);err=abs(p1-p0);p0=p1;p1,err,k,y=feval('f1041',p1)if(errdelta)|(y==0),break,endp1,err,k,y=feval('f1041',p1)end*********************************************************************四、程序实例与计算结果例用上述程序求方程0233xx的一个近似解，给定初值为1.2，误差界为610。解：先用m文件先定义二个名为f1041.m和df1041.m的函数文件。functiony=f1041(x)y=x^3–3*x+2;functiony=df1041(x)y=3*x^2-3;建立一个主程序prog1041.mclearnewton('f1041','df1041',1.2,10^(-6),18)然后在MATLAB命令窗口运行上述主程序，即：prog1041计算结果如下：p0=1.2000ans=0.1280p1=1.1030err=0.0970k=1y=0.0329p1=1.1030err=0.0970k=1y=0.0329p1=1.0524err=0.0507k=2y=0.0084p1=1.0524err=0.0507k=2y=0.0084p1=1.0264err=0.0260k=3y=0.0021p1=1.0264err=0.0260k=3y=0.0021p1=1.0133err=0.0131k=4y=5.2963e-004p1=1.0133err=0.0131k=4y=5.2963e-004p1=1.0066err=0.0066k=5y=1.3270e-004p1=1.0066err=0.0066k=5y=1.3270e-004p1=1.0033err=0.0033k=6y=3.3211e-005p1=1.0033err=0.0033k=6y=3.3211e-005p1=1.0017err=0.0017k=7y=8.3074e-006p1=1.0017err=0.0017k=7y=8.3074e-006p1=1.0008err=8.3157e-004k=8y=2.0774e-006p1=1.0008err=8.3157e-004k=8y=2.0774e-006p1=1.0004err=4.1596e-004k=9y=5.1943e-007p1=1.0004err=4.1596e-004k=9y=5.1943e-007p1=1.0002err=2.0802e-004k=10y=1.2987e-007p1=1.0002err=2.0802e-004k=10y=1.2987e-007p1=1.0001err=1.0402e-004k=11y=3.2468e-008p1=1.0001err=1.0402e-004k=11y=3.2468e-008p1=1.0001err=5.2014e-005k=12y=8.1170e-009p1=1.0001err=5.2014e-005k=12y=8.1170e-009p1=1.0000err=2.6008e-005k=13y=2.0293e-009p1=1.0000err=2.6008e-005k=13y=2.0293e-009p1=1.0000err=1.3004e-005k=14y=5.0732e-010p1=1.0000err=1.3004e-005k=14y=5.0732e-010p1=1.0000err=6.5020e-006k=15y=1.2683e-010p1=1.0000err=6.5020e-006k=15y=1.2683e-010p1=1.0000err=3.2510e-006k=16y=3.1708e-011p1=1.0000err=3.2510e-006k=16y=3.1708e-011p1=1.0000err=1.6255e-006k=17y=7.9270e-012p1=1.0000err=1.6255e-006k=17y=7.9270e-012p1=1.0000err=8.1277e-007k=18y=1.9817e-012ans=1.0000这说明，经过18次迭代得到满足精度要求的值。以下是程序运行截图：