計算機(jī)數(shù)值方法試驗報告

1、v1.0可編輯可修改40TAIYUAN UNlVEBSITkF OF TECHNOLOGY本科實(shí)驗報告課程名稱：計算機(jī)數(shù)值方法實(shí)驗項目：方程求根線性方程組的直接解法線性方程組的迭代解法代數(shù)插值和最小二乘法擬合多項式實(shí)驗地點(diǎn)：逸夫302專業(yè)班級：軟件學(xué)號：學(xué)生姓名：指導(dǎo)教師：田華2013年4 月24 日學(xué)生姓名實(shí)驗成績實(shí)驗名稱實(shí)驗一方程求根實(shí)驗?zāi)康暮鸵?必填)熟悉使用、迭代法、牛頓法、割線法等方法對給定的方程進(jìn)行根的求解。選擇 上述方法中的兩種方法求方程：一分法 f(x)=x +4x-10=0在1,2內(nèi)的一個實(shí)根，*5且要求滿足精度|x -Xn|< X 10-實(shí)驗內(nèi)容和原理(必填)函數(shù)f

2、(x)在區(qū)間(x, y)上連續(xù)，先在區(qū)間(x, y)確定a與b，若f(a) , f(b)異號，說明在區(qū)間(a , b)內(nèi)存在零點(diǎn)，然后求 f(a+b)/2。 假設(shè) F(a)<0,F(b)>0,a<b, 如果f(a+b)/2=0 ，該點(diǎn)即為零點(diǎn)； 如果 f(a+b)/2<0，則區(qū)間(a+b)/2 , b)內(nèi)存在零點(diǎn)，(a+b)/2 > a; 如果 f(a+b)/2>0，則區(qū)間(a,(a+b)/2) 內(nèi)存在零點(diǎn)，(a+b)/2 < b;返回重新循環(huán)，不斷接近零點(diǎn)。通過每次把f(x)的零點(diǎn)所在區(qū)間收縮一半的方法，使區(qū)間內(nèi)的兩個端點(diǎn)逐步逼近函數(shù)零點(diǎn)，最終求得零

3、點(diǎn)近似值。主要儀器設(shè)備臺式或筆記本計算機(jī)實(shí)驗記錄(寫出實(shí)驗內(nèi)容中的程序代碼和運(yùn)行結(jié)果)(可分欄或加頁)代碼1.二分法：#in clude<>#in clude<>#in clude<> int mai n()double a=, b=; double x,s;while(1)x=(a+b)/2;s=pow(x,3)+4*x*x-10;if < s && s <break;else if(s < 0)a=x;else if(s > 0)b=x;prin tf("%ft%fn",a,b);prin tf

4、("%fn",x);prin tf("%fn",s);return 0;2.割線法：#i nclude""#i nclude""int mai n()float c,a=,b=;while(1)c=b-(b*b*b+4*b*b-10)*(b-a)/(b*b*b+4*b*b-(a*a*a+4*a*a); if(fabs(b-c)<*break;b=c;prin tf("%fn",b);prin tf("%fn",c);流程圖;運(yùn)行結(jié)果；1二分法1 2500001 5000

5、001 2500001.3750001-3125001 3750001 3437501.3750001 3593751 2750001.3593751.3671881363281136718813632811 3652341 3642581.3652341 3647461 3652341 3649901 36523413651121 36523413651731.3652341.3652041.3652341.3652191.3652341.3652271 3652341 3652273652311.3652291.3652311.3652301.365230 0 0000011.365231

6、13872561.360731 1366160 1 3650391 3652701 3652221 3652321 3652302割線法實(shí)驗結(jié)果和分析兩種方法均能求出方程的解，但割線法比二分法的收 斂速度更快，且程序的代碼更簡潔心得體會（遇到的問題和解決方法）通過實(shí)驗，加深了對方程求根方法的理解，加強(qiáng)了實(shí)踐操作能力，實(shí)現(xiàn)了理論和實(shí) 踐相結(jié)合。實(shí)驗名稱實(shí)驗二線性方程組的直接求解實(shí)驗?zāi)康暮鸵蠛侠砝肎auss消元法、LU分解法、追趕法求解下列方程組:123X114012X28241X3130.3 101559.1431X159.175.2916.13012X246.7811.2952X31121

7、1X424215X1287210X274836X371261120X4321X17121X25(n=5,10,100,)121Xn 1512Xn5實(shí)驗內(nèi)容高斯消元：I ik=aik/a kk3ij =3ij - l ik *3kj (k=1,2, ，n-1i=k+1,k+2,n j=k+1,k+2,n+1 )由回代過程求得原方程組的解：Xn = a nn+l/ a nn X k = ( a kn+1-刀 akj X j )/ a kk追趕法：當(dāng)矩陣A為三對角矩陣，在A的LU分解中，L取下三角陣，U取單位上二角陣，這樣求解方程組Ax=d的方法稱為追趕法LU分解法：將系數(shù)矩陣A轉(zhuǎn)化為A=L*U,

8、L為單位下三角矩陣，U為普通上三角矩陣，然后通過解方程組l*y=b，u*x=y,來求解x。主要儀器設(shè)備臺式或筆記本計算機(jī)實(shí)驗記錄（寫出實(shí)驗內(nèi)容中的程序代碼和運(yùn)行結(jié)果）（可分欄或加頁）咼斯消兀法#in clude<>#defi ne n 3main ()int i,j,k;float ann ,c nn ,b n,d n;for(i=0;i <n ;i+)for(j=0;j<n;j+)sca nf("%f", &aij); cij=aij;sca nf("%f",&bi);di=bi;for(k=0;k< n;

9、k+)bk=dk/ckk;for(i=0;i <n ;i+)if(i=k) continue; cik=cik/ckk;for(j=k+1;j <n ;j+)akj=ckj/ckk;aij=cij-cik*ckj; bi=di-cik*dk;for(i=0;i <n ;i+)di=bi;for(j=k+1;j <n ;j+)cij=aij;for(i=0;i <n ;i+)prin tf("b%d=%fn",i,bi);LU分解法：#in clude <>#in clude <>#defi ne L 30double a

10、LL,bL,lLL,uLL,xL,yL; int mai n()intn ,i,j,k,r;printf("請輸入矩陣元次：n");sea nf("%d",&n);printf(”請輸入矩陣各項：n");for(i=1;i<=n;+i)for(j=1;j<=n;+j)scan f("%lf",&aij);printf("請輸入方程組的常數(shù)項：n");for(i=1;i<=n;+i)scan f("%lf", &bi);for(i=1;i<

11、=n;+i)for(j=1;j<=n;+j)lij=0;uij=;for(k=1;k<=n;+k)for(j=k;j<=n;+j)ukj=akj;for(r=1;r<k;+r)ukj-=lkr*urj;for(i=k+1;i<=n ;+i)lik=aik; for(r=1;r<k;+r)lik-=lir*urk;lik/= ukk;lkk=;for(i=1;i<=n;+i)yi = bi;for(j=1;j<i;+j)yi-=lij*yj;for(i=n; i>O;-i)xi = yi;for(j=i+1;j<=n ;+j)xi-=u

12、ij*xj; xi/= uii;for(i=1;i<=n;+i)prin tf("%n ”,xi);return 0;追趕法#i nclude ""#defi ne n 5main ()float an ,b n,c n-1,d n,t;int i;scan f("%f%f%f",&b0, &c0, &d0); for(i=1;i< n-1;i+)scan f("%f%f%f%f",&ai, &bi,&ci, &di); scan f("%f%f%

13、f",&an-1,&b n-1，& d n-1); c0=c0/b0;dO=dO/bO;for(i=1;i< n-1;i+)t=bi-ci-1*ai;ci=ci/t;di=(di-di-1*ai)/t;dn-1=(d n-1-d n-2*a n-1)/(b n-1-c n-2*a n-1);for(i=n-2;i>=0;i-) di=di-ci*di+1;for(i=0;i< n;i+)prin tf("d%d=%fn",i,di);實(shí)驗結(jié)果和分析1高斯消元法2LU分解3追趕法實(shí)驗分析：高斯消元法，是先消元，再回帶的過程。

14、由程序段可以發(fā)現(xiàn)，始終消去 對角線下方的元素。從消元過程可以看出，對于n階線性方程組，只要各 步主元素不為零，經(jīng)過n-1步消元，就可以得到一個等價的系數(shù)矩陣為 上三角形陣的方程組，然后再利用回代過程可求得原方程組的解。LU分解法，分解矩陣為單位下三角陣 L與上三角陣U的乘積，然后解方 程組Ly=b，回代，解方程組Ux=yo其中的L為n階單位下三角陣、U為 上三角陣對于追趕法，追趕法是適用于三角矩陣的線性方程組的求解的方法，并 不適用于其他類型矩陣。心得體會（遇到的問題和解決方法）本次實(shí)驗難度比較大，在編譯時經(jīng)常出現(xiàn)各種錯誤，程序代碼也比較繁瑣，深深感 覺到自己的上機(jī)操作能力有限，應(yīng)加強(qiáng)自己的編

15、程能力，以后要繼續(xù)努力。實(shí)驗名稱實(shí)驗三線性方程組的迭代求解實(shí)驗?zāi)康暮鸵笫褂醚趴杀鹊ɑ蚋咚?賽德爾迭代法對下列方程組進(jìn)行求解。410碼七一2如=7.2 一土 *10花 一2益=8.3 a -=4.2實(shí)驗內(nèi)容設(shè)線性方程組Ax=b的系數(shù)矩陣A可逆，且主對角元素ai1,a22,ann均不為零，令D=diag( aii,a 22,a nn)并將A分解成A=(A-D)+D從而線性方程組可寫成Dx=(D-A)x+b則有迭代公式x(k+1)=B1x(k)+f1其中，B1=I-D-1A,f1=D-1b。主要儀器設(shè)備臺式或筆記本計算機(jī)實(shí)驗記錄(寫出實(shí)驗內(nèi)容中的程序代碼和運(yùn)行結(jié)果)(可分欄或加頁)雅可比迭代法

16、#in clude<>#in clude<>#defi ne n 3void f(float *b,float x)float y n+1=0,0,0,1;int i,j,k;dok=0;for(i=0;i <n+1;i+)xi=yi；for(i=0;i <n ;i+)yi=0;for(j=0;j <n+1;j+) yi+=*(b+( n+1)*i+j)*xj; for(i=0;i <n ;i+) if(fabs(yi-xi)<= k+; if(k=3) break;while(1);for(i=0;i <n ;i+)prin tf(

17、"y%d=%fn",i,yi); main ()float bnn+1=0,0,0,;float x n+1=0,0,0,1;f(b0,x);高斯一賽德爾迭代法#i nclude "iostream"#in elude "ioma nip"using n amespace std;int mai n()int i,j,k=0,m,n;double t1,t2,e1,e2=;coutvv"請輸入精度e："cin> >e1;coutvv"請輸入系數(shù)矩陣行數(shù)："cin>>m;c

18、outvv"請輸入系數(shù)矩陣列數(shù)："cin>>n;coutvve ndl;double (*a)=new double *m;生成二維動態(tài)數(shù)組for(i=0;i<=m;i+)ai二new double n;double (*b)二new double m;double (*x)二new double n;cout«"請輸入系數(shù)矩陣："<<endl;cout«""<<e ndl;for(i nt nu ml二 0;nu m1< m;nu m1+)for(i nt nu m

19、2二 0;num2<n;nu m2+)cin> >a nu m1 nu m2;cout«e ndl;cout«"輸入的系數(shù)矩陣為："<<endl;for (int nu m3二 0;nu m3< m;nu m3+)for(i nt num4=0;num4<n;num 4+)cout<<a nu m3 nu m4<<""cout«e ndl;coutvv""<<e ndl;coutvv"請輸入矩陣 b: "&l

20、t;<endl;coutvv""vve ndl;for(i nt nu m5二 0;nu m5v m;nu m5+)cin>>bnu m5;coutvv"輸入的矩陣b為："vvendl;for(i nt nu m6二 0;nu m6v m;nu m6+)coutvvb nu m6vv""coutvve ndl;coutvv""vve ndl;for(i nt nu m7二 0;num7vn;nu m7+)xnu m7=;docoutvv"第"vvkvv"次迭代值：&q

21、uot;e2=;for(i=0;i<m;i+)double sum二；for(j=0;j <n ;j+)if(j!=i)sum+=aij*xj;t仁 xi;t2=e2;xi=(bi-sum)/aii;e2=(xi)-t1>=0(xi)-t1:t1-(xi);e2=(e2>=t2e2:t2);cout«setprecisi on( 8)<<xi<<"cout«e ndl;k+;while(e2>=e1 &&k<30);coutvv"共迭代了 "vvkvv"次&q

22、uot;deletea;deleteb;deletex;return 0 ;實(shí)驗結(jié)果和分析1雅克比迭代2高斯一一賽德爾迭代實(shí)驗分析：使用這兩種方法都可以求出方程的解，高斯賽德爾迭代法所 需的迭代次數(shù)比雅克比迭代少，能夠更早的達(dá)到精度要求。 但是雅克比的時效性要比咼斯賽德爾的好。心得體會（遇到的問題和解決方法）本次試驗，讓我對這兩種方法更加理解，在編程操作上也更加熟練，今后繼續(xù)努力, 不斷豐富自己的知識，增強(qiáng)操作能力。實(shí)驗名稱實(shí)驗四代數(shù)插值和最小二乘法擬合實(shí)驗?zāi)康暮鸵?使用拉格朗日插值法或牛頓插值法求解：已知f(x)在6個點(diǎn)的函數(shù)值如下表所 示，運(yùn)用插值方法，求f的近似值。Xf(x)2給定數(shù)據(jù)

23、點(diǎn)(Xi , yi)，用最小二乘法擬合數(shù)據(jù)的多項式，并求平方誤差。Xioyi1實(shí)驗內(nèi)容1設(shè)函數(shù)在區(qū)間a,b上n+1互異節(jié)點(diǎn)Xo, Xi,Xn上的函數(shù)值分別為yo, yi,yn,求n次插值多項式Pn(x),滿足條件Pn(Xj)=yj, j=o,1,n令Ln(x)=y ol o(X)+y il i(x)+ +yn(x)= 刀 yi l i(x)其中l(wèi) o(x),l l(x),l n(x)為以Xo,X 1,X n為節(jié)點(diǎn)的n次插值基函數(shù)，則Ln(x)是一次數(shù)不超過n的多項式，且滿足Ln(Xj)=yj, L=o,1,n再由插值多項式的唯一性，得Pn(x)三 Ln(x)2建立正規(guī)方程組：j+kjE(E X

24、i ) ak=EXi y , j=0,1,n平方誤差：-k2I =E(E akXi -y i)對給定數(shù)據(jù)點(diǎn)(Xi , Yi)(i=0,1，m),在取定的函數(shù)類 中，求p(x)，使誤差的平方和 E最小，E2二刀p(Xi) -Yi 2。從幾何意義上講，就是尋求與給定點(diǎn)(Xi ，Yi)(i=0,1，m)的距離平方和為最小的曲線y=p(x)。函數(shù)p(x )稱為擬合函數(shù)或最小二乘解，求擬合函數(shù)p(x )的方法稱為曲線擬合的最小二乘法。得到的兩個關(guān)于a0、a1為未知數(shù)的兩個方程組，解這兩個方程組得出：a0 =(刀 Y i) / m - a1 (刀 Xi) / ma1 = m 刀 Xi Yi -(刀 Xi

25、刀 Yi) / m 刀Xi2 -(刀 Xi)2 )即最終的擬合多項式各項系數(shù)主要儀器設(shè)備臺式或筆記本計算機(jī)實(shí)驗記錄(寫出實(shí)驗內(nèi)容中的程序代碼和運(yùn)行結(jié)果)(可分欄或加頁)1#in elude <>#in elude <>#in elude <>#in elude <>void differe nce(float *x,float *y,i nt n)float *f;int k,i;f=(float *)malloc( n*sizeof(float);for(k=1;k<=n ;k+)f0=yk;for(i=0;i<k;i+)fi+1=

26、(fi-yi)/(xk-xi);yk=fk;return;int mai n()int i,n;float x20,y20,xx,yy;printf("請輸入數(shù)據(jù)個數(shù)n: ”);scan f("%d",&n);prin tf("n");for(i=0;i<=n _1;i+)prin tf("x%d=",i);sca nf("%f", &xi);prin tf("y%d=",i);sca nf("%f", &yi);printf(&quo

27、t;n");prin tf("n ”);differe nce(x,(float *)y,n);printf("請輸入插值xx:");scan f("%f', &xx);yy=y20;for(i=n-1;i>=0;i-)yy=yy*(xx-xi)+yi;prin tf("n近似值為：(f)=%fn",xx,yy);2#i nclude<>#in clude<>#defi ne N 15double power(double &a,i nt n)double b=1;for

28、(i nt i=0;i <n ;i+)b*=a;return b;void Gauss();double XN,YN,sumXN,sumYN,aNN,bN,lNN,xN;int mai n()ofstream outdata;ifstream in data;double s;int i,j,k, n,in de x;cout<<"請輸入已知點(diǎn)的個數(shù) n="cin»n;cout<<e ndl;cout<<"請輸入 X和Y:"«endl;/輸入給定數(shù)據(jù)for(i=0;i< n;i+)cou

29、t<<"X"<<i<<"="ci n> >Xi;sumX1+=Xi;cout<<"Y"<<i<<"="ci n»Yi;sumY1+=Yi; cout<<e ndl;cout<<"sumX1="<<sumX1<<"t"<<"sumY1="<<sumY1<<e ndl;cout<

30、<"請輸入擬合次數(shù) index="cin»in dex;cout<<e ndl;i=n;sumX0=i;for(i=2;i<=2* in dex;i+)sumXi=0;for(j=0;j< n;j+)sumXi+=power(Xj,i);cout<<"sumX"<<i<<"="<<sumXi<<e ndl;for(i=2;i<=i ndex+1;i+)sumYi=O;for(j=0;j< n;j+)sumYi+=power(

31、Xj,i-1)*Yj;cout<<"sumY"<<i<<"="<<sumYi<<e ndl;for(i=1;i<=i ndex+1;i+)/建立正規(guī)方程組for(j=1;j<=i ndex+1;j+)aij=sumXi+j-2;bi=sumYi;k=1;/用高斯消元法解方程組dofor(j=k+1;j<=i ndex+1;j+)ljk=ajk/akk;for(i=k+1;i<=i ndex+1;i+)for(j=k+1;j<=i ndex+1;j+)aij=aij-lik*akj;bi=bi-lik*bk;if(k=i ndex+1)break;k+;while(1);xin dex+1=bi ndex+1/ai ndex+1 in dex