非線性方程的數(shù)值計(jì)算方法實(shí)驗(yàn)

1、非線性方程的數(shù)值計(jì)算方法實(shí)驗(yàn)【摘要】在利用數(shù)學(xué)工具研究社會(huì)現(xiàn)象和自然現(xiàn)象，或解決工程技術(shù)等問題時(shí)，很多問題都可以歸結(jié)為非線性方程的求解問題，無論在理論研究方面還是在實(shí)際應(yīng)用中，求解非線性方程都占了非常重要的地位。綜合當(dāng)前各類非線性方程的數(shù)值解法，通過比較分析，二分法，迭代法，牛頓拉夫森方法，迭代法的收斂階和加速收斂方法，以上的算法應(yīng)用對(duì)某個(gè)具體實(shí)際問題選擇相應(yīng)的數(shù)值解法。關(guān) 鍵 詞 非線性方程；二分法；迭代法；牛頓-拉夫森法；割線法等。1、 實(shí)驗(yàn)?zāi)康?通過本實(shí)驗(yàn)的學(xué)習(xí)，應(yīng)掌握非線性方程的數(shù)值解法的基本思想和原理，深刻認(rèn)識(shí)現(xiàn)實(shí)中非線性方程數(shù)值的意義；明確代數(shù)精度的概念；掌握二分法、不動(dòng)點(diǎn)迭代法、

2、牛頓迭代法、割線法等常用的解非線性方程的方法；培養(yǎng)編程與上機(jī)調(diào)試能力。2、 實(shí)驗(yàn)原理二分法：?jiǎn)巫兞亢瘮?shù)方程： f（x）=0其中，f(x)在閉區(qū)間a，b上連續(xù)、單調(diào)，且f(a)*f(b)<0,則有函數(shù)的介值定理可知，方程f（x）=0在（a，b）區(qū)間內(nèi)有且只有一個(gè)解，二分法是通過函數(shù)在區(qū)間端點(diǎn)的符號(hào)來確定所在區(qū)域，將有根區(qū)間縮小到充分小，從而可以求出滿足給定精度的根的近似值。下面研究二分法的幾何意義： 設(shè)=1, =b, 區(qū)間，中點(diǎn)= 及，若=0，則，若 f()*f()<0，令=，=，則根 ,中，這樣就得到長(zhǎng)度縮小一半的有根區(qū)間,，若 f()*f()<0，令=，=，則根 ,中，這樣

3、就得到長(zhǎng)度縮小一半的有根區(qū)間,即f()f()<0,此時(shí)-=,對(duì)有根區(qū)間,重復(fù)上述步驟，即分半求中點(diǎn)，判斷中電處符號(hào)，則可得長(zhǎng)度有縮小一半的有根區(qū)間, 如圖所示： 重復(fù)上述過程，第n步就得到根的近似序列及包含的區(qū)間套，如下：（1）（2）（3）-=(4) 且|-| (n=1,2,3.)顯然lim，且以等比數(shù)列的收斂速度收斂于，因此用二分法求f（x）=0的實(shí)根可以達(dá)到任意指定精度。（二）迭代法：對(duì)給定的方程，將它轉(zhuǎn)達(dá)換成等價(jià)形式：。給定初值，由此來構(gòu)造迭代序列，如果迭代法收斂，即，有，則就是方程的根。在計(jì)算中當(dāng)小于給定的精度控制量時(shí)，取為方程的根。（三）牛頓迭代法：設(shè)方程f(x)=0在其根的某

4、個(gè)領(lǐng)域U(，)內(nèi)有一階連續(xù)導(dǎo)數(shù)，且f() 0。求f(x)=0的根，首先要將f(x)=0轉(zhuǎn)化為等價(jià)形式，并使 (x)滿足不動(dòng)點(diǎn)迭代的一般理論。 于是我們令 (x)=x+h(x)f(x),可由 ()=0來確定h(x)的結(jié)構(gòu)，根據(jù)(x)=1+h()f()+h()f(x1)=1+h()f()=0可得h()=-1/f() ，由于f(x) 0，且f(x) 連續(xù)，因此當(dāng)h(x)=-1/f(x) 時(shí)， h(x1)=0，即令 (x)=x-f(x)/f(x)， 從而有迭代格式 = （k=0,1,2，.）由于, , .都在U領(lǐng)域里，從而當(dāng)B比較小時(shí)，可用f()可近似代替f(),= - ，此方法稱為牛頓迭代法。（四）

5、割線法：設(shè),為方程f(x)=0的兩個(gè)近似根。用差商得：f()-f()/ - ，代替牛頓迭代公式中的導(dǎo)數(shù)f()，于是得到如下的迭代公式： =-。下面研究割線法的幾何意義： 經(jīng)過點(diǎn)（,f()）及點(diǎn)（,f()）兩點(diǎn)作割線，其點(diǎn)斜式方程為： Y=f（）- ，其零點(diǎn)為X=- 把X用表示即得到迭代格式，需要兩個(gè)初值此割線與 X 軸交點(diǎn)的橫坐標(biāo)就是新的近似值，如圖所示。 三、實(shí)驗(yàn)內(nèi)容（1） 、實(shí)驗(yàn)描述 1.P40.1:參照程序2.1求解出單調(diào)收斂的不動(dòng)點(diǎn)。 2.P49.1:已知初值，時(shí)間和末值，求解匯率。 4.P69.1：已知運(yùn)動(dòng)方程求解運(yùn)動(dòng)時(shí)間和距離。 （二）、實(shí)驗(yàn)題目 1.使用程序2.1求解下面每個(gè)函數(shù)

6、的不動(dòng)點(diǎn)（盡肯能多）近似值，答案精確到 小數(shù)點(diǎn)后12位。同時(shí)，構(gòu)造每個(gè)函數(shù)的圖和直線y=x來顯示所有不動(dòng)點(diǎn)。(d) 2. 如果在240個(gè)月內(nèi)每月付款300美元，求解滿足全部現(xiàn)金A為500000美元的 匯率I的近似值（精確到小數(shù)點(diǎn)后10位）。4.設(shè)投射方的運(yùn)動(dòng)程為 （a） 求當(dāng)撞擊地面時(shí)經(jīng)過的時(shí)間，精確到小數(shù)點(diǎn)后10位。（b） 求水平飛行行程，精確到小數(shù)點(diǎn)后10位。四、實(shí)驗(yàn)結(jié)果及分析 1. P40.1（d）： 算法：（1）輸入函數(shù)g，p0,tol,max1,令k=2。 （2）判斷k>max1是否成立，如成立輸出結(jié)果，如不成立，執(zhí)行（3）。 （3）令p(k)=g(p(k-1),err=|p(

7、k)-p(k-1)|。 (4)判斷err<tol是否成立，如成立，輸出結(jié)果，如不成立，令k=k+1,跳轉(zhuǎn)到（2）。startInput g,p0,tol;k=2k>max1YNLet p(k)=g(p(k-1);err=|p(k)-p(k-1)|k=k+1err<to1outputYNoutputend令，簡(jiǎn)單畫出與的圖像，可令x=1,。圖4. 與當(dāng)輸入k,p,err,P=fixpt(g,1,1e-12,100)輸出k =2 p =1 err =0P =1 1 通過圖形我們知道不動(dòng)點(diǎn)有2個(gè)，X=1是吸引不動(dòng)點(diǎn).因而函數(shù)只有一個(gè)吸引迭代點(diǎn)。算法：（1） 輸入f,a,b,delt

8、a,max1,令k= 1。（2） f(a)=feval(f,a);f(b)=feval(f,b);k=1.（3）判斷f(a)*f(b)>0是否成立，若不成立輸出結(jié)果，若成立執(zhí)行步驟（4）。（4）判斷k>max1是否成立,若成立輸出結(jié)果，若不成立執(zhí)行步驟（5）。（5）令dx=yb*(b-a)/(yb-ya),c=b-x,yc=feval(f,c). (6) 判斷yc=0是否成立，若成立輸出結(jié)果，若不成立執(zhí)行步驟（7）。（7）判斷yc*yb>0是否成立，若成立，令b=c,yb=yc.若不成立,令a=c,ya=yc.( 8)令err=|b-a|.(9)判斷err<delta是

9、否成立，若成立，輸出結(jié)果。若不成立，令k=k+1,執(zhí)行（4）。start流程圖：Input f,a,b,delta,max1，k=k+1Yf(a)*f(b)>0NoutputYk>max1Ndx=yb*(b-a)/(b-a)c=b-xyc=feval(f,c)Yyc=0NYNyc*yb>0b=cyb=yca=cya=ycerr=|b-a|Nerr<deltaYoutputend當(dāng)輸入c,err,yc=regula(f,0.15,0.16,1e-10,100)輸出err =-1.257285475730896e-008yc =2.775557561562891e-017算

10、法：（1） 輸入h,r,p0,p1,delta,max1.令k=1.（2） 判斷k>max1是否成立，若成立，輸出；若不成立，執(zhí)行（3）。（3） 令p2=p1-feval(f,p1)*(p1-p0)/(feval(f,p1)-feval(f,p0),err1=|p1-p0|,err2=|feval(r,p1)-feval(r,p0)|,p0=p1,p1=p2 (4) 判斷err1<delta與err2<delta是否成立，若成立，輸出；若不成立，令k=k+1,執(zhí)行步驟（2）。start流程圖：Input f,r,p0,p1,delta max1 k=1k>max1P2=

11、p1-feval(f,p1)*(p1-p0)/(feval(h,p1)-feval(h,p0)Err1=|p1-p0| err2=|feval(r,p1)-feval(r,p0)|C=feval(r,p1) p0=p1 p1=p2Err2<delta&err2<deltaoutputend輸入p1,err1,err2,k,c=secant(h,r,9,9.5,1e-10,100)輸出p1 = 9.087899668785003err2 =2.523847797419876e-011k = 6c = 1.090547960254200e+003五、實(shí)驗(yàn)結(jié)論1.采用的不動(dòng)點(diǎn)迭代

12、法，只能計(jì)算出吸引不動(dòng)點(diǎn)，不能算出排斥迭代點(diǎn)。 2.計(jì)算得所求的利率為。4. 計(jì)算得投射體撞擊地面經(jīng)過的時(shí)間p1 = 9.087899668785003飛行行程c = 1.090547960254200e+003。附件（代碼）：1(d). plot the figure of x.(x-cos(x) with y=xx=-2:0.05:2.4;y=x;g=x.(x-cos(x);plot(x,g,x,y)xlabel('x');ylabel('x.(x-cos(x) y=x');%definite function g(x)function b=g(x)b=x.

13、(x-cos(x);function k,p,err,P=fixpt(g,p0,tol,max1)%input-g is iteratin function input as string 'g'% -p0 is initial guess for the fixed point% -tol is the tolerance% -max1 is the maximum number of iterations%output-k is the number of iteratins that were carried out% -p is approximation to the

14、 fixed point% -err is the error in the approximation % -P contains the sequencepnP(1)=p0;for k=2:max1 P(k)=feval(g,P(k-1); err=abs(P(k)-P(k-1); relerr=err/(abs(P(k)+eps); p=P(k); if (err<tol)|(relerr<tol),break,endendif k=max1 disp('maximum number of iterations exceeded')endP=P'2.

15、%definite the function f(x)function A=f(x)p=300;n=240;A=12*p*(1+x/12)n-1)/x-500000;function c,err,yc=regula(f,a,b,delta,max1)%input-f is the function input as a string 'f'% -a and b are the left and right points% -delta is the tolerance for the zero% -max1 is the maximum number of iterations

16、%output-c is the zero% -yc=f(c)% -err is the error estimate for cya=feval(f,a);yb=feval(f,b);if ya*yb>0 disp('ya,yb are not suitable ') returnend for k=1:max1 dx=yb*(b-a)/(yb-ya); c=b-dx; yc=feval(f,c); if yc=0,break; elseif yc*yb<0 a=c;ya=yc; else b=c;yb=yc; end err=abs(b-a); if err&l

17、t;delta break,endend4.%definite function r(t)function x=r(t)x=2400*(1-exp(-t/15);%definite function h(t)function y=h(t)y=9600*(1-exp(-t/15)-480*t;function p1,err1,err2,k,c=secant(h,r,p0,p1,delta,max1)%input-h is the function of altitude inpput as string 'h'% -r is the function of path input as a string 'r'% -p0 and p1 is the initial approximations to a zero% -delta is the tolerance for p1 and y% -max1 is the maximum number of iterations%output-p1 is the secant method approximation to the zero% -err1 is estimate for p1%