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

1、數(shù)值計(jì)算方法非線性方程的數(shù)值計(jì)算方法實(shí)驗(yàn)一、 實(shí)驗(yàn)描述：在科學(xué)研究和工程實(shí)踐中，經(jīng)常需要求解大量的非線性方程。本實(shí)驗(yàn)正是通過(guò)計(jì)算機(jī)的程序設(shè)計(jì)，使用迭代法、波爾查諾二分法、試值法、牛頓-拉夫森法和割線法，來(lái)實(shí)現(xiàn)非線性方程的求解。本實(shí)驗(yàn)中通過(guò)對(duì)各種方法的實(shí)踐運(yùn)用，可以比較出各種方法的優(yōu)缺點(diǎn)。并且，通過(guò)完成實(shí)驗(yàn)，可加深對(duì)各種方法的原理的理解，熟悉掌握C語(yǔ)言在這些方法中的運(yùn)用。二、實(shí)驗(yàn)內(nèi)容：1、 求函數(shù)的不動(dòng)點(diǎn)（盡可能多）近似值，答案精確到小數(shù)點(diǎn)后12位；2、 如果在240個(gè)月內(nèi)每月付款300美元，求解滿足全部年金A為500000美元的利率I，的近似值（精確到小數(shù)點(diǎn)后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、 設(shè)投射體的運(yùn)動(dòng)方程為： y=f(t)=9600(1-e-t/15)-480t x=r(t)=2400(1-e-t/15)(a)求當(dāng)撞擊地面時(shí)經(jīng)過(guò)的時(shí)間，精確到小數(shù)點(diǎn)后10位。 (b)求水平飛行行程，精確到小數(shù)點(diǎn)后10位。三、實(shí)驗(yàn)原理： (1)、不動(dòng)點(diǎn)迭代法：它是一種逐次逼近的方法,即用某個(gè)固定公式反復(fù)校正根的近似值,使之逐步精確化,最后得到滿足精度

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

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

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

6、且存在數(shù)pa,b,滿足f(p)=0。如果f(p)0,則存在一個(gè)數(shù)>0，對(duì)任意初始近似值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。這樣，通過(guò)尋找函數(shù)的不動(dòng)點(diǎn)，可以實(shí)現(xiàn)尋找方程f(x)=0的根的牛頓-拉夫森迭代。附：定理(牛頓-拉夫森迭代的加速收斂)： 設(shè)牛頓-拉夫森算法產(chǎn)生的序列線性收斂到M階根x=p,其中M>1,則牛頓-拉夫森迭代公式： pk=p

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

8、精度取為0.000000000001。從迭代過(guò)程的誤差收斂速度可以看出不動(dòng)點(diǎn)迭代的誤差收斂較慢，所需迭代次數(shù)較多。2、滿足條件的利率I的計(jì)算(試值法)：(a)、計(jì)算結(jié)果：(b)、結(jié)果分析：此題經(jīng)過(guò)一定的計(jì)算分析，將書中所給公式結(jié)合題中的條件可得出最終的計(jì)算式：3600I1+I12240-1-500000=0,求解其中的I即為利率。又經(jīng)過(guò)一定的計(jì)算分析，取其初始計(jì)算區(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的近似值。其中每個(gè)函數(shù)的M、p0和最終結(jié)果p已給出。(a)、計(jì)算結(jié)果：(b)、結(jié)果分析： 由于C語(yǔ)言函數(shù)庫(kù)中沒有求導(dǎo)功能函數(shù)，而此題所要求用的加速牛頓-拉夫森算法的計(jì)算公式 pk=pk-1-M f(pk-1)f(pk-1) 中要用到題中所給函數(shù)的導(dǎo)函數(shù)，所以求出它們的導(dǎo)函數(shù)分別為f(x)=5(x-2)2, f(x)=3x2cos(x3)和f(x)=In(x)-1- 1x 。此題的誤差精度取為0.00000000001。從輸出結(jié)果中可以大致判斷出，加速牛頓-拉夫森法的收斂速度非常快。4、給出投射物體的運(yùn)動(dòng)方程y=f(t)=9600(1-e-t/15)-480t和x=r(

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

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

12、e pN,m,errN,relerrN;int k,n;p0=0.879;/給定一個(gè)合理的初始值以此來(lái)進(jìn)行迭代。/for(k=0;k<N;k+)m=pk-cos(pk);pk+1=pow(pk,m);/進(jìn)行迭代的方程，由書上給出。/errk=fabs(pk+1-pk);/對(duì)絕對(duì)誤差的運(yùn)算。/relerrk=errk/(fabs(pk+1)+FLT_EPSILON);/對(duì)相對(duì)誤差的運(yùn)算。/n=k;/把k賦給n,使下面輸出時(shí)能夠在合適的時(shí)候停下來(lái)。/if(errk<pre|relerrk<pre) break;/判斷循環(huán)在何處達(dá)到規(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);/對(duì)迭代點(diǎn)及誤差的輸出。/getch();return 0;第二題：#include<stdio.h>#include<math.h>#include<float.h>#define N 10000#define pre 0.000000000001/使循環(huán)能夠在達(dá)到某一誤差精度后停下來(lái)。/double main()double f(double x);double aN,bN,

14、IN,errN,relerrN;int i,n;a0=0.10,b0=0.20;/預(yù)估它們的初始區(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);/對(duì)絕對(duì)誤差的運(yùn)算。/ relerri=erri/(fabs(Ii)+FLT_EPSILON);/對(duì)相對(duì)誤差的運(yùn)算。/ n=i;/將i值賦給n,使輸出能夠在合

15、適的時(shí)候停下來(lái)。/ if(erri<pre|relerri<pre) break;/判斷循環(huán)在何處達(dá)到規(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);/輸出利率的計(jì)算過(guò)程。/ getch(); return 0;double f(doub

16、le x) return(3600/x)*(pow(1+x/12,240)-1)-500000);/此題的計(jì)算公式，由書上給出。/第三題：#include<stdio.h>#include<math.h>#include<float.h>/為了調(diào)用FLT_EPSILON，防止出現(xiàn)分母為零的情況。/#define pre 0.00000000001/使循環(huán)能夠在達(dá)到某一誤差精度后停下來(lái)。/#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;/定義結(jié)果數(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);/對(duì)絕對(duì)誤差的運(yùn)算。k加2，以保證下面輸出時(shí)其有足夠多的值輸出。/i=k;/將循環(huán)停止時(shí)的k值賦給i，便于在下面輸出時(shí)有恰當(dāng)?shù)闹凳馆敵鐾Ｖ埂Ｏ峦?if(errk<

18、pre)break;for(k=0;k<N;k+)bk+1=bk-3.0*fb(bk);errk=fabs(0.0-bk);/對(duì)絕對(duì)誤差的運(yùn)算。/j=k;if(errk<pre)break;for(k=0;k<N;k+)ck+1=ck-2.0*fc(ck);errk=fabs(1.0-ck);/對(duì)絕對(duì)誤差的運(yùn)算。/m=k;if(errk<pre)break; if(i>j) j=i;if(m>j) m=j;/對(duì)i,m,j進(jìn)行比較，并求出最大值，并用于輸出。/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);/書上所給的計(jì)算要求的公式。下同。/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ù)學(xué)函數(shù)的調(diào)用。/#include&l

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