matlab非线性方程求解要点

1非线性方程的解法1引言数学物理中的许多问题归结为解函数方程的问题，即，0)(xf(1.1)这里，)(xf可以是代数多项式，也可以是超越函数。若有数*x为方程0)(xf的根，或称函数)(xf的零点。设函数)(xf在],[ba内连续，且0)()(bfaf。根据连续函数的性质知道，方程0)(xf在区间],[ba内至少有一个实根；我们又知道，方程0)(xf的根，除了极少简单方程的根可以用解析式表达外，一般方程的根很难用一个式子表达。即使能表示成解析式的，往往也很复杂，不便计算。所以，具体求根时，一般先寻求根的某一个初始近似值，然后再将初始近似值逐步加工成满足精度要求为止。如何寻求根的初始值呢？简单述之，为了明确起见，不妨设)(xf在区间],[ba内有一个实的单根，且0)(,0)(bfaf。我们从左端出点ax0出发，按某一预定的步长h一步一步地向右跨，每跨一步进行一次根的“搜索”，即检查每一步的起点kx和1kx(即，hxk)的函数值是否同号。若有：0)(*)(hxfxfkk(1.2)那么所求的根必在),(hxxkk内，这时可取kx或hxk作为根的初始近似值。这种方法通常称为“定步长搜索法”。另外，还是图解法、近似方程法和解析法。2迭代法2.1迭代法的一般概念迭代法是数值计算中一类典型方法，不仅用于方程求根，而且用于方程组求解，矩阵求特征值等方面。迭代法的基本思想是一种逐次逼近的方法。首先取一个精糙的近似值，然后2用同一个递推公式，反复校正这个初值，直到满足预先给定的精度要求为止。对于迭代法，一般需要讨论的基本问题是：迭代法的构造、迭代序列的收敛性天收敛速度以及误差估计。这里，主要看看解方程迭代式的构造。对方程(1.1)，在区间],[ba内，可改写成为：)(xx(2.1)取],[0bax，用递推公式：)(1kkxx,,2,1,0k(2.2)可得到序列：0210}{,,,,kkkxxxxx(2.3)当k时，序列0}{kkx有极限x~，且)(x在x~附近连续，则在式(2.2)两边极限，得，)~(~xx即，x~为方程(2.1)的根。由于方式(1.1)和方程(2.1)等价，所以，xx~*即，*limxxkk式(2.2)称为迭代式，也称为迭代公式；)(x可称为迭代函数。称求得的序列0}{kkx为迭代序列。2.2程序和实例下面是基于MATLAB的迭代法程序，用迭代格式)(1nnxgp，求解方程)(xgx，其中初始值为0p。**************************************************************************function[p,k,err,P]=fixpt(f1021,p0,tol,max1)%f1021是给定的迭代函数。%p0是给定的初始值。3%tol是给定的误差界。%max1是所允许的最大迭代次数。%k是所进行的迭代次数加1。%p是不动点的近似值。%err是误差。%P={p1,p2,…,pn}P(1)=p0;fork=2:max1P(k)=feval('f1021',P(k-1));k,err=abs(P(k)-P(k-1))p=P(k);if(errtol),break;endifk==max1disp('maximumnumberofiterationsexceeded');endendP=P;****************************************************************************例2.1用上述程序求方程0sin2xx的一个近似解，给定初始值5.00x，误差界为510。解：先用m文件先定义一个名为f1021.m的函数文件。functiony=f1021(x)y=sin(x)/x;建立一个主程序prog1021.mclcclearallfixpt('f1021',0.5,10^(-5),20)然后在MATLAB命令窗口运行上述主程序，即：prog1021计算结果如下。k=2err=0.45894k=3err=0.1052k=4err=0.0292k=5err=0.0078k=6err=0.0021k=7err=5.7408e-004k=8err=1.5525e-004k=9err=4.1975e-005k=10err=1.1350e-005k=11err=3.0688e-006P=Columns1through60.50000.95890.85370.88290.87510.8772Columns7through110.87660.87680.87670.87670.8767ans=0.87673二分法3.1二分法原理二分法是方程求解最直观、最简单的方法。二分法以连续函数的介值定理为基础的。由介值定理知道，若函数)(xf区间],[ba上连续，且0)(*)(bfaf，即)(af和)(bf负号相反，5则)(xf在],[ba内一定有实根。二分法的基本思想是：用对分区间的方法根据分点处函数)(xf的符号逐步将有限区间缩小，使在足够小的区间内，方程有且仅有一根。下面简述其基本步骤。首先记bbaa00,。用中点2000bax将区间],[00ba等分成2个小区间：],[00xa和],[00bx。然后分析可能存在的三种情况：如果0)(*)(0xfaf，则0x是零点，也就是方程的根。如果0)(*)(0xfaf，则区间],[00xa内存在零点。如果0)(*)(0bfxf，则区间],[00bx内存在零点。对有根的新区间施行同样的操作，于是得到一系列有空的区间：],[],[],[],[221100kkbabababa(3.1)其中每1个区间的长度都是前一区间长度的一半，最后1个区间的长度为：kkkabab2(3.2)如果取最后1个区间],[kkba的中点：2kkkabx(3.3)作为0)(xf根的近似值，则有误差估计式：1*22kkkkababxx(3.4)对于所给精度，若取k使得12kab(3.5)则有，kxx*(3.6)3.2程序与实例用二分法求解方程0)(xf在有根区间],[ba内的一个根，其中)(xf在],[ba只有一个根的情形。6***************************************************************************function[c,err,yc]=bisect(f1031,a,b,delta)%f1031是所要求解的函数。%a和b分别为有根区间的左右限。%delta是允许的误差界。%c为所求近似解。%yc为函数f在c上的值。%err是c的误差估计。ya=feval('f1031',a);yb=feval('f',b);ifyb==0,c=b；returnendifya*yb0,disp('(a,b)不是有根区间');returnendmax1=1+round((log(b-a)-log(delta))/log(2));fork=1:max1c=(a+b)/2;yc=feval('f',c);ifyc==0a=c;b=c;returnelseifyb*yc0b=c;yb=yc;elsea=c;ya=yc;endif(b-a)delta,return,endendk;c=(a+b)/2;err=abs(b-a);yc=feval('f',c);*****************************************************************************7例3.1用上述程序求方程013xx在区间]5.1,0.1[内的一个近似解，要求准确到小数点后第2位。解：先用m文件先定义一个名为f1031.m的函数文件。functiony=f1031(x)y=x^3-x-1;建立一个主程序prog1031.mclcclearallbisect('f1031',1.0,1.5,0.0005)然后在MATLAB命令窗口运行上述主程序，即：prog1031计算结果如下。ans=1.32474牛顿法4.1牛顿法原理从前面迭代法，我们知道，迭代函数)(x构造的好坏，不仅影响收敛速度，而且迭代格式有可能发散。怎样选择一个迭代函数才能保证迭代序列一定收敛呢？构代迭代函数的一条重要途径是用近似方程来代替原方程去求根。因此，如果能将非线性方程(1.1)用线性方程去代替，那么，求近似根问题就很容易解决，而且十分方便。牛顿(Newton)法就是一种将非线性方程线化的一种方法。设kx是方程(1.1)的一个近似根，把如果)(xf在kx处作一阶Taylor展开，即：))((')()(kkkxxxfxfxf(4.1)于是我们得到如下近似方程：0))((')(kkkxxxfxf(4.2)设0)('kxf，则方程(10.2.1)的解为：8()'()kkkfxxxfx(4.3)取x~作为原方程(1.1)的新近似根1kx，即令：)(')(1kkkkxfxfxx，,2,1,0k(4.4)上式称为牛顿迭代格式。用牛顿迭代格式求方程的根的方法就称为牛顿迭代法，简称牛顿法。牛顿法具有明显的几何意义。方程：))((')(kkkxxxfxfy(10.4.5)是曲线)(xfy上点))(,(kkxfx处的切线方程。迭代格式(4.4)就是用切线式(4.5)的零点来代替曲线(1.1)的零点。正因为如此，牛顿法也称为切线法。牛顿迭代法对单根至少是二阶局部收敛的，而对于重根是一阶局部收敛的。一般来说，牛顿法对初值0x的要求较高，初值足够靠近*x时才能保证收敛。若要保证初值在较大范围内收敛，则需对)(xf加一些条件。如果所加的条件不满足，而导致牛顿法不收敛时，则需对牛顿法作一些改时，即可以采用下面的迭代格式：)(')(1kkkkxfxfxx，,2,1,0k(4.6)式中，10，称为下山因子。因此，用这种方法求方程的根，也称为牛顿下山法。牛顿法对单根收敛速度快，但每迭代一次，除需计算)(kxf之外，还要计算)('kxf的值。如果)(xf比较复杂，计算)('kxf的工作量就可能比较大。为了避免计算导数值，我们可用差商来代替导数。通常用如下几种方法：(1)割线法。如果用11)()(kkkkxxxfxf代替)('kxf，则得到割线法的迭代格式为：)()()(111kkkkkkkxfxfxfxxxx(4.7)(2)拟牛顿法。如果用9)())(()(1kkkkxfxfxfxf代替)('kxf，则得到拟牛顿法的迭代格式为：))(()()(121kkkkkkxfxfxfxfxx(4.8)(3)Steffenson法。如果用)()())((kkkkxfxfxfxf代替)('kxf，则得到拟牛顿法的迭代格式为：)())(()(21kkkkkkxfxfxfxfxx(4.9)4.2程序与实例1.牛顿法的程序给定初值0p，用牛顿法格式)(')(111kkkkpfpfpp，,2,1k，求解非线性方程0)(xf。**************************************************************************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)10fork=1:max1p1=p0-feval('f1041',p0)/feval('