非線性方程的數(shù)值計算方法實驗

1、數(shù)值計算方法非線性方程的數(shù)值計算方法實驗一、 實驗描述：在科學研究和工程實踐中，經(jīng)常需要求解大量的非線性方程。本實驗正是通過計算機的程序設計，使用迭代法、波爾查諾二分法、試值法、牛頓-拉夫森法和割線法，來實現(xiàn)非線性方程的求解。本實驗中通過對各種方法的實踐運用，可以比較出各種方法的優(yōu)缺點。并且，通過完成實驗，可加深對各種方法的原理的理解，熟悉掌握C語言在這些方法中的運用。二、實驗內(nèi)容：1、 求函數(shù)的不動點（盡可能多）近似值，答案精確到小數(shù)點后12位；2、 如果在240個月內(nèi)每月付款300美元，求解滿足全部年金A為500000美元的利率I，的近似值（精確到小數(shù)點后10位）。3、 利用加速牛頓-拉夫

2、森算法，用其求下列函數(shù)M階根p的近似值。(a)、f(x)=(x-2)5，M=5,p=2，初始值p0=1。(b)、f(x)=sin(x3)，M=3,p=0，初始值p0=1。(c)、f(x)=(x-1)ln(x),M=2,p=1，初始值p0=2。 4、 設投射體的運動方程為： y=f(t)=9600(1-e-t/15)-480t x=r(t)=2400(1-e-t/15)(a)求當撞擊地面時經(jīng)過的時間，精確到小數(shù)點后10位。 (b)求水平飛行行程，精確到小數(shù)點后10位。三、實驗原理： (1)、不動點迭代法：它是一種逐次逼近的方法,即用某個固定公式反復校正根的近似值,使之逐步精確化,最后得到滿足精度

3、要求的結果。它利用計算機運算速度快，適合做重復性操作的特點，讓計算機對一個函數(shù)進行重復執(zhí)行，在每次執(zhí)行這個函數(shù)時，都從變量的原值推出它的一個新值，直至推出最終答案為止。 迭代法一般可用于尋找不動點，即：存在一個實數(shù)P，滿足P=g(P)，則稱P為函數(shù)g(x)的一個不動點。且有定理：若g(x)是一個連續(xù)函數(shù)，且pnn=0是由不動點迭代生成的序列。如果limnpn=P,則P是g(x)的不動點。所以，不動點的尋找多用迭代法。(2)、波爾查諾二分法：起始區(qū)間a,b必須滿足f(a)與f(b)的符號相反的條件。由于連續(xù)函數(shù)y=f(x)的圖形無間斷，所以它會在零點x=r處跨過x軸，且r在區(qū)間內(nèi)。通過二分法可將

4、區(qū)間內(nèi)的端點逐步逼近零點，直到得到一個任意小的包含零點的間隔。二分法定理：設fC(a,b),且存在數(shù)ra,b滿足f(r)=0。如果f(a)和f(b)的符號相反，且cnn=0為二分法生成的中點序列，則：r-cnb-a2n+1 其中n=0,1, (1)這樣，序列cnn=0收斂到零點x=r即可表示為：limncn=r (2)(3)、試值法：假設一個函數(shù)中，有f(a)和f(b)符號相反。二分法使用區(qū)間a,b的中點進行下一次迭代。如果找到經(jīng)過點(a,f(a)和(b,f(b)的割線L與x軸的交點(c,0),則可得到一個更好的近似值。為了尋找值c，定義了線L的斜率m的兩種表示方法，一種表示方法為：m=fb-

5、f(a)b-a (3)這里使用了點(a,f(a)和(b,f(b)。另一種表示方法為：m=0-f(b)c-b (4)這里使用了點(c,0)和(b,f(b)。 使式(3)和式(4)的斜率相等，則有：fb-f(a)b-a=0-f(b)c-b (5)為了更容易求解c，可進一步表示為：c=b-fb(b-a)fb-f(a) (6)這樣會出現(xiàn)3種可能性： 如果f(a)和f(c)的符號相反，則在a,c內(nèi)有一個零點。 如果f(c)和f(b)的符號相反，則在c,b內(nèi)有一個零點。如果f(c)=0，則c是零點。然后，可按二分法的方法進行下一步運算。(4)、牛頓-拉夫森法： 此法根據(jù)，牛頓-拉夫森定理：設fC2a,b,

6、且存在數(shù)pa,b,滿足f(p)=0。如果f(p)0,則存在一個數(shù)>0，對任意初始近似值p0p-,p+，使得由如下迭代定義的序列pkk=0收斂到p: pk=g(pk-1)= pk-1f(pk-1)f(pk-1) 其中k=1,2, (7)其中，函數(shù)g(x)由如下定義：g(x)=x-f(x)f(x) (8)且被稱為牛頓-拉夫森迭代函數(shù)。由于f(p)=0,顯然g(p)=p。這樣，通過尋找函數(shù)的不動點，可以實現(xiàn)尋找方程f(x)=0的根的牛頓-拉夫森迭代。附：定理(牛頓-拉夫森迭代的加速收斂)： 設牛頓-拉夫森算法產(chǎn)生的序列線性收斂到M階根x=p,其中M>1,則牛頓-拉夫森迭代公式： pk=p

7、k-1-M f(pk-1)f(pk-1) (9)(5)、割線法： 割線法包含的公式與試值法的公式一樣，只是在關于如何定義每個后續(xù)項的邏輯判定上不一樣。需要兩個靠近點(p,0)的初始點(p0,f(p0)和(p1,f(p1)。 可根據(jù)兩點迭代法公式，得到一般項: pk+1=g(pk, pk-1)= pk-fpk(pk-pk-1)fpk-f(pk-1) (10)四、結果計算及分析：1、函數(shù)g(x)=xx-cos(x)的不動點的迭代(迭代法)： (a)計算結果：(b)結果分析:此題經(jīng)過matlab圖形仿真，可以看出其解的大致范圍在0.8,1.2之間。 此結果是在取p0=0.879的情況下得出的，其誤差

8、精度取為0.000000000001。從迭代過程的誤差收斂速度可以看出不動點迭代的誤差收斂較慢，所需迭代次數(shù)較多。2、滿足條件的利率I的計算(試值法)：(a)、計算結果：(b)、結果分析：此題經(jīng)過一定的計算分析，將書中所給公式結合題中的條件可得出最終的計算式：3600I1+I12240-1-500000=0,求解其中的I即為利率。又經(jīng)過一定的計算分析，取其初始計算區(qū)間為0.1,0.2，誤差精度取為0.000000000001。從其誤差收斂的速度可以看出，試值法的誤差收斂速度快，所需的迭代次數(shù)少。3、利用加速牛頓-拉夫森算法，求函數(shù)f(x)=(x-2)5、f(x)=sin(x3)和f(x)=(x

9、-1)ln(x)的M階根p的近似值。其中每個函數(shù)的M、p0和最終結果p已給出。(a)、計算結果：(b)、結果分析： 由于C語言函數(shù)庫中沒有求導功能函數(shù)，而此題所要求用的加速牛頓-拉夫森算法的計算公式 pk=pk-1-M f(pk-1)f(pk-1) 中要用到題中所給函數(shù)的導函數(shù)，所以求出它們的導函數(shù)分別為f(x)=5(x-2)2, f(x)=3x2cos(x3)和f(x)=In(x)-1- 1x 。此題的誤差精度取為0.00000000001。從輸出結果中可以大致判斷出，加速牛頓-拉夫森法的收斂速度非常快。4、給出投射物體的運動方程y=f(t)=9600(1-e-t/15)-480t和x=r(

10、t)=2400(1-e-t/15)求物體落地時經(jīng)過的時間和水平飛行程。(割線法) (a)計算結果：(b)、結果分析：此題用了割線法的思路進行程序設計，由于割線法的公式要求，需要給出t0,t1 兩個初始值。經(jīng)過一定的分析，給出t0=5.0，t1=15.0，誤差精度設為0.0000000001。從輸出結果中可以看出，割線法的誤差收斂速度比較快，所需的迭代次數(shù)比較少。五、實驗結果分析： 經(jīng)過以上四個實驗的比較分析，結合書上的實例與原理分析，可以得出。不動點迭代法比起其他的迭代法，迭代速度慢、次數(shù)多。試值法的迭代次相對少一些，但若對其初始區(qū)間估算得太大的話，會導致迭代步數(shù)增多，甚至會導致發(fā)散。割線法比

11、試值法要快得多，可以達到相對較高的精度，且在其迭代過程中每步只需一次新的函數(shù)賦值，其收斂階一般能達到1.1618。加速牛頓-拉夫森法再求多重根時速度很快，收斂階能達到2。附件：第一題：#include<stdio.h>#include<math.h>#include<float.h>/為了調(diào)用FLT_EPSILON，防止出現(xiàn)分母為零的情況。/#define N 1000/保證有足夠多的運算次數(shù)能達到比較精確的結果，且便于修改。/#define pre 0.0000000000001/使循環(huán)能夠在達到某一誤差精度后停下來。/double main()doubl

12、e pN,m,errN,relerrN;int k,n;p0=0.879;/給定一個合理的初始值以此來進行迭代。/for(k=0;k<N;k+)m=pk-cos(pk);pk+1=pow(pk,m);/進行迭代的方程，由書上給出。/errk=fabs(pk+1-pk);/對絕對誤差的運算。/relerrk=errk/(fabs(pk+1)+FLT_EPSILON);/對相對誤差的運算。/n=k;/把k賦給n,使下面輸出時能夠在合適的時候停下來。/if(errk<pre|relerrk<pre) break;/判斷循環(huán)在何處達到規(guī)定精度，并跳出循環(huán)。/for(k=0;k<

13、n;k+)printf("nX%d=%14.12f err%d=%14.12f relerr%d=%14.12f",k,pk,k,errk,k,relerrk);/對迭代點及誤差的輸出。/getch();return 0;第二題：#include<stdio.h>#include<math.h>#include<float.h>#define N 10000#define pre 0.000000000001/使循環(huán)能夠在達到某一誤差精度后停下來。/double main()double f(double x);double aN,bN,

14、IN,errN,relerrN;int i,n;a0=0.10,b0=0.20;/預估它們的初始區(qū)間/for(i=0;i<N;i+)Ii=bi-(f(bi)*(bi-ai)/(f(bi)-f(ai);/試值法公式，有書上給出。/ if(f(ai)*f(Ii)>0)ai+1=Ii,bi+1=bi;/根據(jù)試值法，判斷區(qū)間。/else ai+1=ai,bi+1=Ii;/根據(jù)試值法，判斷區(qū)間。/ erri=fabs(Ii-Ii-1);/對絕對誤差的運算。/ relerri=erri/(fabs(Ii)+FLT_EPSILON);/對相對誤差的運算。/ n=i;/將i值賦給n,使輸出能夠在合

15、適的時候停下來。/ if(erri<pre|relerri<pre) break;/判斷循環(huán)在何處達到規(guī)定精度，并跳出循環(huán)。/ err0=I0-0;/定義初始誤差，便于下面輸出。/ relerr0=err0/(fabs(I0)+FLT_EPSILON);/定義初始誤差，便于下面輸出。/ for(i=0;i<n;i+) printf("n I%d=%12.10f err%d=%12.10f relerr%d=%12.10f",i,Ii,i,erri,i,relerri);/輸出利率的計算過程。/ getch(); return 0;double f(doub

16、le x) return(3600/x)*(pow(1+x/12,240)-1)-500000);/此題的計算公式，由書上給出。/第三題：#include<stdio.h>#include<math.h>#include<float.h>/為了調(diào)用FLT_EPSILON，防止出現(xiàn)分母為零的情況。/#define pre 0.00000000001/使循環(huán)能夠在達到某一誤差精度后停下來。/#define N 1000/迭代次數(shù)，以保證能得到最精確的值。/double main() double fa(double x);/聲明fa函數(shù)。/ double fb(

17、double x);/聲明fb函數(shù)。/ double fc(double x);/聲明fc函數(shù)。/ double aN,bN,cN,errN;/定義結果數(shù)組，與誤差數(shù)組。/ int k,i,j,m; a0=1.0,b0=1.0,c0=2.0;/賦初始值，由書上給出。/ for(k=0;k<N;k+)ak+1=ak-5.0*fa(ak);/加速牛頓-拉夫森公式，見書上P64。下同。/errk=fabs(2.0-ak+2);/對絕對誤差的運算。k加2，以保證下面輸出時其有足夠多的值輸出。/i=k;/將循環(huán)停止時的k值賦給i，便于在下面輸出時有恰當?shù)闹凳馆敵鐾Ｖ埂Ｏ峦?if(errk<

18、pre)break;for(k=0;k<N;k+)bk+1=bk-3.0*fb(bk);errk=fabs(0.0-bk);/對絕對誤差的運算。/j=k;if(errk<pre)break;for(k=0;k<N;k+)ck+1=ck-2.0*fc(ck);errk=fabs(1.0-ck);/對絕對誤差的運算。/m=k;if(errk<pre)break; if(i>j) j=i;if(m>j) m=j;/對i,m,j進行比較，并求出最大值，并用于輸出。/for(k=0;k<m;k+)printf("na%d=%12.10f b%d=%12

19、.10f c%d=%12.10f",k,ak,k,bk,k,ck);getch();return 0;double fa(double x)return(1.0/5.0)*(x-2.0);/書上所給的計算要求的公式。下同。/double fb(double x)return(1.0/(3.0*x*x)*tan(x*x*x);double fc(double x)return(x-1.0)*log(x)/(log(x)+1.0-1.0/x);第四題：#include<stdio.h>#include<math.h>/方便各種數(shù)學函數(shù)的調(diào)用。/#include&l

20、t;float.h>/為了調(diào)用FLT_EPSILON，防止出現(xiàn)分母為零的情況。/#define N 1000/保證有足夠多的運算次數(shù)能達到比較精確的結果，且便于修改。/#define pre 0.0000000001/使循環(huán)能夠在達到某一誤差精度后停下來。/double main()double f(double t);/聲明f函數(shù)。/double r(double t);/聲明r函數(shù)。/double tN,errN,relerrN,difN;/使結果能夠以數(shù)組形式輸出，便于分析。/char ch1="t",ch2="dif",ch3="err",ch4="relerr"/對輸出進行標注。/int k,i;/定義i，以便在達到所需精度之后，對所有結果以及對最終結果進行輸出。/t0=5.0,t1=15.0;/給出任意合法的初始值。/err0=fabs(t1-t0),relerr0=err0/fabs(t1+FLT_EPSILON),dif0=t1-t0;/對各個結果的初始值進行運算。/for(k=1;k<N;k+)tk+1=tk-(f(tk)*(tk-tk-1)/(f(tk)-f(tk