1、非線性方程的解法1引言數學物理中的許多問題歸結為解函數方程的問題，即，f(x)0(1.1)這里，f(x)可以是代數多項式，也可以是超越函數。若有數x*為方程f(x)0的根，或稱函數f(x)的零點。設函數f(x)在a,b內連續，且f(a)f(b)0。根據連續函數的性質知道，方程f(x)0在區間a,b內至少有一個實根；我們又知道，方程f(x)0的根，除了極少簡單方程的根可以用解析式表達外，一般方程的根很難用一個式子表達。即使能表示成解析式的，往往也很復雜，不便計算。所以，具體求根時，一般先尋求根的某一個初始近似值，然后再將初始近似值逐步加工成滿足精度要求為止。如何尋求根的初始值呢？簡單述之，為了明

2、確起見，不妨設f(x)在區間a,b內有一個實的單根，且f(a)0,f(b)0。我們從左端出點xa出發，按某一預定的步長h一步一步地0向右跨，每跨一步進行一次根的“搜索”，即檢查每一步的起點x和xkk1(即，xh)的函數值k是否同號。若有：f(x)*f(xh)0(1.2)kk那么所求的根必在(x,xh)內，這時可取x或xh作為根的初始近似值。這種方法通常稱kkkk為“定步長搜索法”。另外，還是圖解法、近似方程法和解析法。2迭代法2.1迭代法的一般概念迭代法是數值計算中一類典型方法，不僅用于方程求根，而且用于方程組求解，矩陣求特征值等方面。迭代法的基本思想是一種逐次逼近的方法。首先取一個精糙的近似

3、值，然后1/17用同一個遞推公式，反復校正這個初值，直到滿足預先給定的精度要求為止。對于迭代法，一般需要討論的基本問題是：迭代法的構造、迭代序列的收斂性天收斂速度以及誤差估計。這里，主要看看解方程迭代式的構造。對方程(1.1)，在區間a,b內，可改寫成為：x(x)(2.1)取xa,b，用遞推公式：0 xk1(x),k0,1,2,(2.2)k可得到序列：x,x,x,x,x012kkk0(2.3)kxx當k時，序列x有極限，且(x)在附近連續，則在式(2.2)兩邊極限，得，kk0 xx()x即，為方程(2.1)的根。由于方式(1.1)和方程(2.1)等價，所以，xx*即，limxx*k式(2.2)

4、稱為迭代式，也稱為迭代公式；(x)可稱為迭代函數。稱求得的序列x為迭代序kk0列。2.2程序和實例下面是基于MATLAB的迭代法程序，用迭代格式pn1g(x)，求解方程xg(x)，其中n初始值為p。0*functionp,k,err,P=fixpt(f1021,p0,tol,max1)%f1021是給定的迭代函數。%p0是給定的初始值。2/17%tol是給定的誤差界。%max1是所允許的最大迭代次數。%k是所進行的迭代次數加1。%p是不動點的近似值。%err是誤差。%P=p1,p2,pnP(1)=p0;fork=2:max1P(k)=feval(f1021,P(k-1);k,err=abs(P

5、(k)-P(k-1)p=P(k);if(errprog1021計算結果如下。k=2err=0.45893/17k=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.8767

6、ans=0.87673二分法3.1二分法原理二分法是方程求解最直觀、最簡單的方法。二分法以連續函數的介值定理為基礎的。由介值定理知道，若函數f(x)區間a,b上連續，且f(a)*f(b)0，即f(a)和f(b)負號相反，4/1722k則f(x)在a,b內一定有實根。二分法的基本思想是：用對分區間的方法根據分點處函數f(x)的符號逐步將有限區間縮小，使在足夠小的區間內，方程有且僅有一根。下面簡述其基本步驟。首先記aa,bb。用中點xa0b0將區間a,b等分成2個小區間：a,x和0000000 x,b。然后分析可能存在的三種情況：00如果f(a)*f(x)0，則x是零點，也就是方程的根。00如果f

7、(a)*f(x)0，則區間a,x內存在零點。000如果f(x)*f(b)0，則區間x,b內存在零點。000對有根的新區間施行同樣的操作，于是得到一系列有空的區間：a,ba,ba,ba,b(3.1)001122kk其中每1個區間的長度都是前一區間長度的一半，最后1個區間的長度為：baba(3.2)kk如果取最后1個區間a,b的中點：kk2xbkakk(3.3)x*xba22k1作為f(x)0根的近似值，則有誤差估計式：bakkk對于所給精度，若取k使得(3.4)ba2k1(3.5)則有，x*x(3.6)k3.2程序與實例用二分法求解方程f(x)0在有根區間a,b內的一個根，其中f(x)在a,b只

8、有一個根的情形。5/17*functionc,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=0

9、a=c;b=c;returnelseifyb*yc0b=c;yb=yc;elsea=c;ya=yc;endif(b-a)prog1031計算結果如下。ans=1.32474牛頓法4.1牛頓法原理從前面迭代法，我們知道，迭代函數(x)構造的好壞，不僅影響收斂速度，而且迭代格式有可能發散。怎樣選擇一個迭代函數才能保證迭代序列一定收斂呢？構代迭代函數的一條重要途徑是用近似方程來代替原方程去求根。因此，如果能將非線性方程(1.1)用線性方程去代替，那么，求近似根問題就很容易解決，而且十分方便。牛頓(Newton)法就是一種將非線性方程線化的一種方法。設x是方程(1.1)的一個近似根，把如果f(x)在x

10、處作一階Taylor展開，即：kkf(x)f(x)f(x)(xx)(4.1)kkk于是我們得到如下近似方程：f(x)f(x)(xx)0(4.2)kkk設f(x)0，則方程(10.2.1)的解為：k7/17xxkf(x)kf(x)k(4.3)x取作為原方程(1.1)的新近似根xk1，即令：f(x)xk1xf(xk)，k0,1,2,(4.4)kk上式稱為牛頓迭代格式。用牛頓迭代格式求方程的根的方法就稱為牛頓迭代法，簡稱牛頓法。牛頓法具有明顯的幾何意義。方程：yf(x)f(x)(xx)(10.4.5)kkk是曲線yf(x)上點(x,f(x)處的切線方程。迭代格式(4.4)就是用切線式(4.5)的零點

11、來代替kk曲線(1.1)的零點。正因為如此，牛頓法也稱為切線法。牛頓迭代法對單根至少是二階局部收斂的，而對于重根是一階局部收斂的。一般來說，牛頓法對初值x的要求較高，初值足夠靠近x*時才能保證收斂。若要保證初值在較大范圍內0收斂，則需對f(x)加一些條件。如果所加的條件不滿足，而導致牛頓法不收斂時，則需對牛頓法作一些改時，即可以采用下面的迭代格式：f(x)xk1xf(xk)，k0,1,2,kk(4.6)式中，01，稱為下山因子。因此，用這種方法求方程的根，也稱為牛頓下山法。牛頓法對單根收斂速度快，但每迭代一次，除需計算f(x)之外，還要計算f(x)的值。kk如果f(x)比較復雜，計算f(x)的

12、工作量就可能比較大。為了避免計算導數值，我們可用差k商來代替導數。通常用如下幾種方法：(1)割線法。如果用f(x)f(xkk1xxkk1)代替f(x)，則得到割線法的迭代格式為：kf(x)f(xxk1xxkxk1kkk1)f(x)(4.7)k(2)擬牛頓法。如果用8/17f(x)f(xf(xkkk1f(x)k)代替f(x)，則得到擬牛頓法的迭代格式為：kf(x)f(xf(x)xk1f2(x)kx(4.8)kkkk1(3)Steffenson法。如果用f(xf(x)f(x)kkkf(x)k代替f(x)，則得到擬牛頓法的迭代格式為：kf(xf(x)f(x)xk1f2(x)kx(4.9)kkkk4.

13、2程序與實例1.牛頓法的程序給定初值p，用牛頓法格式pp0kk1f(p)k1f(p)k1，k1,2,，求解非線性方程f(x)0。*functionp1,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)9/17fork=1:max1p1=p0-feval(f1041,p0)/feval(df1041,p0

14、);err=abs(p1-p0);p0=p1;p1,err,k,y=feval(f1041,p1)if(errprog1041計算結果如下。p0=1.2000ans=0.1280p1=1.1030err=0.0970k=1y=0.0329p1=1.103010/17err=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.

15、2963e-004p1=1.0133err=0.0131k=4y=5.2963e-004p1=1.0066err=0.0066k=511/17y=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-006p

16、1=1.0008err=8.3157e-004k=8y=2.0774e-006p1=1.000412/17err=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

17、=5.2014e-005k=12y=8.1170e-009p1=1.0001err=5.2014e-005k=1213/17y=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

18、-006k=15y=1.2683e-010p1=1.0000err=3.2510e-006k=16y=3.1708e-011p1=1.000014/17err=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次迭代得到滿足精度要求的值。2.割線法的MATLAB實現*functionp1,err,k,y=secant(f1042,p0,p1,delta,max1)%f1042是給定的非線性函數。%p0,p1為初始值。%delta為給定誤差界。%max1是迭代次數的上限。%p1為所求得的方程的近似解。%err為p1-p0的絕對值。%k為所需的迭代次數。%yf(p1)p0,p1,feval(f1042,p0),feval(f1042,p1),k=0;fork=1:max115/17p2=p1-feval(f104