數值計算方法上機實驗報告

數值計算方法上機實驗報告實驗目的：復習和鞏固數值計算方法的基本數學模型，全面掌握運用計算機進行數值計算的具體過程及相關問題。利用計算機語言獨立編寫、調試數值計算方法程序，培養(yǎng)學生利用計算機和所學理論知識分析解決實際問題的能力。上機練習任務：利用計算機基本C語言編寫并調試一系列數值方法計算通用程序，并能正確計算給定題目，掌握調試技能。掌握文件使用編程技能，如文件的各類操作，數據格式設計、通用程序運行過程中文件輸入輸出運行方式設計等。一、各算法的算法原理及計算機程序框圖列主元高斯消去法算法原理：高斯消去法是利用現行方程組初等變換中的一種變換，即用一個不為零的數乘一個方程后加只另一個方程，使方程組變成同解的上三角方程組，然后再自下而上對上三角方程組求解。列選住院是當高斯消元到第k步時，從k列的a炒以下（包括a炒）的各元素中選出絕對值最大的，然后通過行交換將其交換到a^的位置上。交換系數矩陣中的兩行（包括常數項），只相當于兩個方程的位置交換了，因此，列選主元不影響求解的結果。

?計算機程序框圖如上?源程序：#defineN200#include"stdio.h"#include"math.h"FILE*fp1,*fp2;voidLZ()(intn,i,j,k=0,l;doubled,t,t1;staticdoublex[N],a[N][N];fp1二fopen(〃a1.txt〃，〃r〃)；fp2二fopen(〃b1.txt〃,〃w〃)；fscanf(fp1,〃%d〃,&n);for(i=0;i<n;++i)for(j=0;j<=n;++j)(fscanf(fp1,〃％lf〃，&a[i][j]);do(d=a[k][k];l二k;i=k+1;do(if(fabs(a[i][k])>fabs(d)) /*選主元*/(d=a[i][k];l=i；}i++;}while(i<n);if(d==0){printf("\n輸入矩陣有誤！\n");}else( /*換行*/if(l!=k)(for(j=k;j<=n;j++)(t=a[l][j];a[l][j]=a[k][j];a[k][j]=t;}}}for(j=k+1;j<=n;j++) /*正消*/a[k][j]/=a[k][k];for(i=k+1;i<n;i++)for(j=k+1;j<=n;j++)a[i][j]-=a[i][k]*a[k][j];k++;}while(k<n);if(k!=0)(for(i=n-1;i>=0;i--) /*回代*/(t1=0;for(j=i+1;j<n;j++)t1+=a[i][j]*x[j];x[i]=a[i][n]-t1;}

}for(i=0;i<n;i++)fprintf(fp2,"\n方程組的根為x[%d]=%lf”,i+1,x[i]);fclose(fp1);fclose(fp2);}main(){LZ();}?具體算例及求解結果：用列選主元法求解下列線性方程組%+2x22x1+x23x%+2x22x1+x23x+2x12-3x3+3x3+x3=22=28輸入32-3813222128?輸入變量、輸入變量:輸出結果：方程組的根為x[1]=6.000000方程組的根為x[2]=4.000000方程組的根為x[3]=2.000000輸出變量說明:a系數矩陣元素’4常向量元素輸出變量:杜里特爾分解法解線性方程?算法原理：求解線性方程組Ax=b時，當對A進行杜里特爾分解，則等價于求解LUx=b，這時可歸結為利用遞推計算相繼求解兩個三角形(系數矩陣為三角矩陣)方程組，用順代，由Ly=b求出y，再利用回帶，由Ux=y求出x。計算機程序框圖：源程序：#include"stdio.h"#include"math.h”FILE*fp1,*fp2;voidmain(){inti,j,k,N;doubles,A[200][200],B[200],x[200],y[200];staticdoubleL[200][200],U[200][200];fp1二fopen(〃a2.txt〃,〃r〃)；fp2二fopen(〃b2.txt〃,〃w〃)；fscanf(fp1,〃％d〃，&N);for(i=0;i<N;i++){for(j=0;j<N;j++)fscanf(fp1,〃％lf〃,&A[i][j]);}for(i=0;i<N;i++)fscanf(fp1,"%lf",&B[i]);for(i=0;i<N;i++) /*LU分解*/(for(j=i;j<N;j++)(s=0.0;for(k=0;k<i;k++)s+=L[i][k]*U[k][j];U[i][j]=A[i][j]-s;}for(j=i+1;j<N;j++)(s=0.0;for(k=0;k<i;k++)s+=U[k][i]*L[j][k];L[j][i]=(A[j][i]-s)/U[i][i];}}for(i=0;i<N;i++)for(j=0;j<N;j++)L[i][i]=1;fprintf(fp2,"\nU矩陣為:")；for(i=0;i<N;i++){fprintf(fp2,〃\n〃)；for(j=0;j<N;j++)fprintf(fp2,〃％10.3f〃，U[i][j]);}fprintf(fp2,"\nL矩陣為:");for(i=0;i<N;i++){fprintf(fp2,〃\n〃)；for(j=0;j<N;j++)fprintf(fp2,〃％10.3f〃，L[i][j]);for(i=0;i<N;i++){s=0.0;for(k=0;k<i;k++)s+=L[i][k]*y[k];y[i]=B[i]-s;}for(i=N-1;i>=0;i--){s=0.0;for(k=i+1;k<N;k++)s+=U[i][k]*x[k];x[i]=(y[i]-s)/U[i][i];}fprintf(fp2,"\n方程組解為：")；for(i=0;i<N;i++)fprintf(fp2,〃\nx%d二%10.3f〃，i+1,x[i]);fclose(fp1);fclose(fp2);}?具體算例及求解結果：用杜里特爾分解法求解方程組TOC\o"1-5"\h\z2x +3x +4x = 39<3x -2x +2x = 144x +2x +3x = 43輸入數據° 輸出結果：U矩陣為:2.0003.0004.0000.000-6.500-4.0000.0000.000-2.5383L矩陣為:2341.0000.0000.0003-221.5001.0000.0004232.0000.6151.000391443方程組解為:TOC\o"1-5"\h\zx1= 6.000x2= 5.000x3= 3.000?輸入變量、輸出變量說明：輸入變量：a^系數矩陣元素，^常向量元素輸出變量：b,b,b解向量元素1 2n拉格朗日插值法?算法原理：首先構造基函數l⑴=H三二，可以證明基函數滿足下列條件:k i=。Lii豐k10lk(Xi)=11對于給定n+1個節(jié)點，n次拉格朗日插值多項式由下式給出:L(x)=云七(x)yk

k=0其中回x-x

l(x)= i-ki=0Xk一氣i豐k由于lk(x)是一個關于x的n次多項式，所以L(x)為關于x的不高于n次的代數多項式。當x=x時，L(x)=y，滿足插值條件。

i ii?計算機程序框圖：開撕■入刀(xyt戶=0,■■■(1.爐？)始■i+lTA,出丁gg源程序：#include〃stdio.h"#include〃math.h〃intn,m,i,j;floatx2,x3,z1=0.0,z=0.0,t,x[50],y[50],c[50],A[50];main()(FILE*fp1,*fp2;fp1二fopen(〃a3.txt〃,〃r〃)；fp2二fopen(〃b3.txt〃,〃w〃)；fscanf(fp1,〃%d〃,&n);for(i=0;i<n;i++)fscanf(fp1,〃%f,%f〃,&x[i],&y[i]);fscanf(fp1,〃%d〃,&m);fscanf(fp1,〃%f〃,&x2);fscanf(fp1,〃%f〃,&x3);for(i=0;i<n;i++) /*選m個最接近的點*/c[i]=fabs(x[i]-x2);for(i=0;i<n;i++)for(j=i+1;j<n;j++)if(c[i]>c[j]) (t=c[i];c[i]=c[j];c[j]=t;t=x[i];x[i]=x[j];x[j]=t;t=y[i];y[i]=y[j];y[j]=t;}for(i=0;i<m;i++) /*求值*/{A[i]=1.0;for(j=0;j<m;j++)if(i!=j)A[i]=A[i]*(x2-x[j])/(x[i]-x[j]);z=z+A[i]*y[i];}for(i=0;i<m;i++) /*求值*/{A[i]=1.0;for(j=0;j<m;j++)if(i!=j)A[i]=A[i]*(x3-x[j])/(x[i]-x[j]);z1=z1+A[i]*y[i];}fprintf(fp2,"\nx二%10.7f處的函數值為：y二%10.7f”,x2,z);fprintf(fp2,"\nx二%10.7f處的函數值為：y二%10.7f”,x3,z1);fclose(fp1);fclose(fp2);}具體算例及求解結果：對于一組數據表進行二次數值插值編程，根據下面數值表計算f(0.49)和f(0.51)x0.20.40.60.8f(x)16201510輸入數據： 輸出結果：x=0.4900000處的函數值為：y=18.86375054 x=0.5100000處的函數值為：y=18.36375240.2,160.4,200.6,150.8,1030.490.51輸入變量、輸出變量說明：輸入變量：（氣,七）插值節(jié)點輸出變量：y插值所得到被插函數在插值點的近似值曲線擬合

?算法原理:對于給定的一組數據(x,y),i=1,2…，m，尋求做n次多項式y(tǒng)—axkkk—0使性能指標)=￡(y_—Ea^k)2為最小。i—1 k=0故上述由于性能指標J可以被看做關于ak，k=0，1，…，n的多元函數，擬合多項式的構造問題可轉化為多元函數的極值問題。令故上述從而有正則方程組da從而有正則方程組mExxEx2Exn]aPEy]乙iEx2iiEx3iiExn+1i0ai—Ex'yiiExn?1- iExn+1?i???Exn+2...?i …Ex2n?i 」an」Ex".k求解即得多項式系數。?計算機程序框圖： ...?源程序： …#include<math.h>#include<stdio.h>main(){inti,j,k,m,n,l,N,t,t1;doublemax,A[50][50],x[50],y[50],S[50],T[50],X[50];floatyb=0.0,xb,a1,a2,a0;FILE*fp1,*fp2;fp1二fopen(〃a4.txt〃,〃r〃)；fp2二fopen(〃b4.txt〃,〃w〃)；fscanf(fp1,"%d%d\n〃,&n,&m);for(i=0;i<n;i++){fscanf(fp1,"%lf%lf〃,&x[i],&y[i]);}fscanf(fp1,〃%f〃,&xb);for(i=0;i<=2*m;i++){S[i]=0.0;for(j=0;j<n;j++)S[i]+=pow(x[j],i);}for(i=0;i<=m;i++){T[i]=0.0;for(j=0;j<n;j++)T[i]+=y[j]*pow(x[j],i);}N=m+1;for(i=0;i<N;i++)(for(j=0;j<N;j++)(l=i+j;A[i][j]=S[l];}A[i][N]=T[i];}for(i=0;i<N-1;i++)(max二fabs(A[i][i]); /*選主*/for(j=i+1;j<N;j++)if(fabs(A[j][i])>max)(max=fabs(A[j][i]);m=j;}if(m!=i)for(k=0;k<=N;k++)(t=A[i][k];A[i][k]=A[m][k];A[m][k]=t;}for(j=i+1;j<N;j++) /*消元*/for(k=N;k>=0;k--)A[j][k]=A[j][k]-A[i][k]*A[j][i]/A[i][i];}for(i=N-1;i>=0;i--) /*回代*/for(j=i-1;j>=0;j--)for(k=N;k>=0;k--)A[j][k]=A[j][k]-A[i][k]*A[j][i]/A[i][i];fprintf(fp2,"\n解為:"); /*輸出結果*/for(i=0;i<N;i++)fprintf(fp2,〃\na%d二%10.7lf〃，i,A[i][N]/A[i][i]);fprintf(fp2,"\n擬合多項式為:\n");fprintf(fp2,〃P(x)二%10.7lf〃，A[0][N]/A[0][0]);for(i=1;i<N;i++)fprintf(fp2,〃+(%10.7lf)x"%d〃，A[i][N]/A[i][i],i);for(i=0;i<N;i++)yb+=(A[i][N]/A[i][i])*pow(xb,i);fprintf(fp2,〃\nP(%f)二%10.7f〃，xb,yb);fclose(fpl);fclose(fp2);}?具體算例及求解結果：對于一組數據表進行二次多項式曲線擬合，根據以下數據胡二次擬合曲線求y(5)Xi12345678910yi1.62.83.64.95.46.87.99.210.211.4試用最小二乘法求二次擬合多項式輸入數據： 102 輸出結果:11.6解為：22.8a0=-0.301245033.6a1=1.316733844.9a2=-0.016643955.4擬合多項式為：6.8 P(x)=-0.3012450+(1.3167338)x"1+(-0.0166439)x"27.9 P(5.000000)=5.8663259TOC\o"1-5"\h\z9.210.211.45.0輸入變量、輸出變量說明：輸入變量：(氣,七)已知數據點輸出變量：a.擬合多項式的系數改進歐拉法?算法原理：當A取值較小時，讓梯形法的迭代公式只迭代一次就結束。這樣先用歐拉公式求得一個初步近似值y”1(。)，稱之為預報值，預報值的精度不高，用它替代梯形法右端的y，再直接計算得出y，并稱之為校正值，這時得到預報-校正公n+1 n+1式。將預報-校正公式

yn+](0)=yn+hf(Xn,yn)hyn+1 n=y+-[_f(xn,yn)+f(xn+1，yyn+1 n稱為改進歐拉公式。?計算機程序框圖：?源程序：#include"stdio.h"#include"math.h"FILE*fp1,*fp2;floatfunc(floatx,floaty)(floatdy;dy=sqrt(2*x*x+3*y*y);/*定義函數的導/*定義函數的導*/}main()(inti;floath,yp,yc,y0,x1,x2;if((fp1=fopen(〃a5.txt〃,〃r〃))==NULL)(printf(〃cannotopenthisfile\n〃)；exit(0);}

fp2二fopen(〃b5.txt〃,〃w〃)；fscanf(fp1,〃％f,%f,%f,%f〃，&x1,&x2,&y0,&h);for(i=0;i<(x2-x1)/h;i++){yp=y0+h*func(x1+i*h,y0);yc=y0+h*func(x1+(i+1)*h,yp);y0=0.5*(yp+yc);fprintf(fp2,〃節(jié)點%6.2f處的值二%10.7f\n〃，x1+(i+1)*h,y0);}fclose(fp1);fclose(fp2);?}?具體算例及求解結果：求解初值問題。取h=0.2,用改進歐拉法求解下列初值問題y/=K+3y2<y(0)=5 ^0<x<20輸入數據：0,20,5,0.2輸出結果節(jié)點0.20處的值=7.03239390.60處的值=13.91481211.00處的值=27.5336847處的值=54.4735222處的值=107.75923542.20處的值=213.15559392.60處的值=421.62611393.00處的值=833.97662353.40處的值=1649.6000366處的值=3262.89355474.20處的值=6453.96997074.60處的值=12765.88525395.00處的值=25250.78710945.40處的值=49945.7949219處的值=98792.27343756.20處的值=195410.09375006.60處的值=386519.14062507.00處的值=764530.81250007.40處的值=1512233.87500007.80處的值=2991183.00000008.20處的值=5916528.25000008.60處的值=11702831.00000009.00處的值=23148077.0000000處的值=45786650.0000000處的值=90565512.0000000節(jié)點節(jié)點節(jié)點節(jié)點節(jié)點節(jié)點0.00處的值=127372260.0000000節(jié)點10.20處的值=179137632.0000000節(jié)點0.40處的值二9.8918471節(jié)點節(jié)點0.80處的值=19.5739155節(jié)點節(jié)點1.20處的值=38.7287178節(jié)點節(jié)點1.60處的值=76.6169243節(jié)點節(jié)點2.00處的值=151.5576096節(jié)點節(jié)點2.40處的值=299.7871094節(jié)點節(jié)點2.80處的值=592.9812927節(jié)點節(jié)點3.20處的值=1172.9145508節(jié)點節(jié)點3.60處的值=2320.0151367節(jié)點節(jié)點4.00處的值=4588.9670410節(jié)點節(jié)點4.4節(jié)點0.20處的值=7.03239390.60處的值=13.91481211.00處的值=27.5336847處的值=54.4735222處的值=107.75923542.20處的值=213.15559392.60處的值=421.62611393.00處的值=833.97662353.40處的值=1649.6000366處的值=3262.89355474.20處的值=6453.96997074.60處的值=12765.88525395.00處的值=25250.78710945.40處的值=49945.7949219處的值=98792.27343756.20處的值=195410.09375006.60處的值=386519.14062507.00處的值=764530.81250007.40處的值=1512233.87500007.80處的值=2991183.00000008.20處的值=5916528.25000008.60處的值=11702831.00000009.00處的值=23148077.0000000處的值=45786650.0000000處的值=90565512.0000000節(jié)點節(jié)點節(jié)點節(jié)點節(jié)點節(jié)點0.00處的值=127372260.0000000節(jié)點10.20處的值=179137632.0000000節(jié)點11.60處的值=1949713344.0000000節(jié)點11.80處的值=2742096640.000000012.00處的值=3856512512.0000000節(jié)點12.20處的值=5423838208.0000000節(jié)點 12.40 處的值=7628141056.0000000 節(jié)點 12.60 處的值=10728294912.0000000節(jié)點12.80處的值=.0000000節(jié)點 13.00 處的值=21220453376.0000000 節(jié)點 13.20 處的值=29844662272.0000000節(jié)點13.40處的值=41973835776.0000000節(jié)點 13.60 處的值=59032424448.0000000 節(jié)點 13.80 處的值=83023802368.0000000節(jié)點14.00處的值=116765528064.0000000節(jié)點 14.20 處的值=164220223488.0000000 節(jié)點 14.40 處的值=230960996352.0000000節(jié)點14.60處的值=324825874432.0000000節(jié)點 14.80 處的值=456838414336.0000000 節(jié)點 15.00 處的值=642502164480.0000000節(jié)點15.20處的值=9.0000000節(jié)點 15.40 處的值=1270862577664.0000000 節(jié)點 15.60 處的值=1787354021888.0000000節(jié)點15.80處的值=25.0000000節(jié)點 16.00 處的值=3535367438336.0000000 節(jié)點 16.20 處的值=4972176474112.0000000節(jié)點16.40處的值=6992919527424.0000000節(jié)點 16.60 處的值=9834913071104.0000000 節(jié)點 16.80 處的值=536.0000000節(jié)點17.00處的值=19453354967040.0000000節(jié)點 17.20 處的值=27359394660352.0000000 節(jié)點 17.40 處的值=38478530215936.0000000節(jié)點17.60處的值=54116592123904.0000000節(jié)點 17.80 處的值=76110125596672.0000000 節(jié)點 18.00 處的值=1456.0000000節(jié)點18.20處的值=8528.0000000節(jié)點18.40處的值=211728063266816.0000000節(jié)點18.60處的值=297776508305408.0000000節(jié)點18.80處的值=4312.0000000節(jié)點19.00處的值=588998829408256.0000000節(jié)點19.20處的值=828373899149312.0000000節(jié)點19.40處的值=11658.0000000節(jié)點19.60處的值=1638514923929600.0000000節(jié)點19.80處的值=23808.0000000節(jié)點20.00處的值=32428.0000000輸入變量、輸出變量說明：輸入變量：3°,*）處置點，h區(qū)間長度，N計算次數輸出變量：（氣,七）初值問題的數值解法結果四階龍格-庫塔法算法原理：用區(qū)間lxk,氣+1］內四個不同點上的函數值的線性組合就得到四階龍格-庫塔法。

四階龍格-庫塔法y=y+h(①k+必k+必k+wk)

n+1n 11223344k=f3,y)nn<k=f(尤+人h,y+pkh)n1n111k=f(尤+人h,y+pkh+pkh)TOC\o"1-5"\h\zn 2 n 211 222k=f(尤+人h,y+pkh+pkh+pkh)n 3 n 311 322 333其中，①,人,日均為待定系數。類似于前面的討論，把k,k,k分別在x點展開成h的幕級數，代入y并2 3 4 n n+1進行花間，然后與y(x)在x點上的泰勒展開式比較，使其兩式比較，使其兩n+1 n式右端直到h4的系數相等，經過復雜的數學演算可得到關于。K日的一組特解十氣=匕=%2=2%1=四31=%2=0<人=日=1］一、33_1廣4一6①=W從而得到下列常用的經典公式h12 3 4y=y+《(k+2k+2k+k)12 3 4h］、+—k)21h］、+2k2)h］、+—k)21h］、+2k2)k=f(x ,yn+1'2k=f(x,yn+1n2k=f(x,y+hk)4 n+1n3經典的龍格-庫塔法每一步需要4次計算函數值f(x,y)，它具有四階精度，即局部截斷誤差是O(h5)。計算機程序框圖：

源程序：#include"stdio.h"#include"math.h"FILE*fp1,*fp2;floatds(floatx,floaty)(floatd;d=sqrt(x*x+y*y);returnd;}voidmain()(inti;floatx[101],y[101],h,b,k1,k2,k3,k4;if((fp1=fopen(〃a6.txt〃,〃r〃))==NULL){printf(〃cannotopenthisfile\n〃)；exit(0);}fp2二fopen(〃b6.txt〃，〃w〃)；fscanf(fp1,〃％f,%f,%f,%f〃，&h,&x[0],&b,&y[0]);

for(i=1;i<=100;i++)if((x[0]+i*h)<=b+h){x[i]=x[0]+i*h;k1=ds(x[i-1],y[i-1]);k2=ds(x[i-1]+h/2,y[i-1]+h*k1/2);k3=ds(x[i-1]+h/2,y[i-1]+h*k2/2);k4=ds(x[i],y[i-1]+h*k3);y[i]=y[i-1]+h*(k1+2*k2+2*k3+k4)/6.0;}elsebreak;for(i=0;i<=100;i++)if((x[0]+i*h)<=b+h)(x[i]=x[0]+i*h;fprintf(fp2,〃\n%10.5f點處的值二%10.7f”,x[i],y[i]);}elsebreak;fclose(fp1);fclose(fp2);}具體算例及求解結果：求解初值問題，取h=0.2，用四階龍格-庫塔法求解下列初值問題<y(0)=6 >輸入數據：0.2,0,20,6輸出結果： 0.00000點處的值二0.40000輸入數據：0.2,0,20,6輸出結果： 0.00000點處的值二0.40000點處的值二8.95238110.80000點處的值=13.3631592點處的值=19.9501495點處的值=29.78152472.00000點處的值=44.4503784點處的值=66.3337555點處的值=98.9784317點處的值=147.6761780點處的值=220.32180794.00000點處的值=328.6933289點處的值=490.3615112點處的值=731.53887946.0000000 0.20000點處的值二7.32860180.60000點處的值=10.93721011.00000點處的值=16.3277931點處的值=24.3755608點處的值=36.3848495點處的值=54.3016205點處的值=81.02948763.00000點處的值=120.9007721點處的值=180.3790741點處的值=269.1072998點處的值=401.4711609點處的值=598.93170175.00000點處的值=893.5048828點處的值=1091.3298340點處的值=1628.07092296.00000點處的值點處的值=1091.3298340點處的值=1628.07092296.00000點處的值=2428.7895508點處的值=3623.3151855點處的值=5405.3295898點處的值=1988.5278320點處的值=2966.5251465點處的值=4425.51855477.00000點處的值=6602.0703125點處的值=8063.7695313點處的值=8063.7695313點處的值=12029.67773448.00000點處的值=17946.0898438點處的值=26772.3007813點處的值=39939.4023438點處的值=9849.0888672點處的值=14693.0488281點處的值=21919.3554688點處的值=32699.68945319.00000點處的值=48781.9882813點處的值=59582.3203125 9.40000點處的值=72773.8437500點處的值=88885.9765625點處的值=88885.976562510.00000點處的值=132601.7031250點處的值=197817.5937500點處的值=295107.8437500點處的值=440247.1875000點處的值=656768.687500012.00000點處的值=979779.3125000點處的值=1461652.3750000點處的值=2180519.2500000點處的值=3252937.7500000點處的值=108565.3359375點處的值=161959.7187500點處的值=241614.406250011.00000點處的值=360444.7187500點處的值=537717.9375000點處的值=802177.2500000點處的值=1196702.5000000點處的值=1785262.250000013.00000點處的值=2663286.2500000點處的值=3973138.2500000點處的值=4852791.000000013.80000點處的值=5927199.000000014.00000點處的值=7239481.000000014.20000點處的值=8842302.0000000點處的值=10799988.000000014.60000點處的值=1