數(shù)值計算方法上機實驗

PAGE1PAGE16重慶交通大學學生實驗報告實驗課程名稱數(shù)值計算方法開課實驗室理學院實驗室學院土木建筑學院專業(yè)班級09級測繪工程專業(yè)2班學生姓名吳必坤學號09210205開課時間2010至2011學年第二學期評分細則評分考勤（10分）排版格式（10分）實驗程序代碼的正確性、完整性（40分）實驗結(jié)果及分析（30分）實驗創(chuàng)新（10分）總成績教師簽名張文忠目錄實驗一誤差——定積分 2實驗二非線性方程求根 4（一）牛頓迭代法 4（二）牛頓下山法 5實驗三解線性代數(shù)方程組——列主元消元法 7實驗四解線性代數(shù)方程組——追趕法 9實驗五插值與擬合——三次樣條差值 11實驗一誤差——定積分實驗內(nèi)容：用兩種方法求解以下積分：,初值2.In=1-In-1*n2實驗過程：1利用C語言編程求解，程序如下：#include<stdio.h>#include<math.h>voidmain(){doublex1;inti,n;x1=100;for(i=50;i>=0;i--){x1=(1-x1)/i;printf("i=%d,xi=%lf\n",i,x1);}}1.2.命令窗口中程序結(jié)果如下：2.1利用C語言編程求解，程序如下：#include<stdio.h>#include<math.h>voidmain(){doublex1=1-exp(-1);inti,n,x0;scanf("%d",&n);x0=1-exp(-1);printf("i=0,x0=0.6320205671\n");for(i=1;i<=n;i++){x1=1-i*x1;printf("i=%d,xi=%12.9lf\n",i,x1);}}1.2.命令窗口中程序結(jié)果如下：實驗結(jié)果分析：從程序運行的結(jié)果可以看出，當輸入n=14時，即運算到n=14時，結(jié)果與書上的結(jié)果并不一致，然而從n=1到n=10的運算結(jié)果是完全一致的，從n=11以后的運算結(jié)果可能是由于不同時期的計算機的處理數(shù)據(jù)的能力的不同或精度不同導致的。實驗二非線性方程求根實驗內(nèi)容：用牛頓迭代法和牛頓下山法求解方程f（x）=，取近似值。1．牛頓迭代法：，直到。2．牛頓下山法：，（λ=1至合適），且直到|實驗過程：1．用牛頓迭代法求解1.1利用C語言編程求解，程序如下：#include<stdio.h>#include<math.h>voidmain(){doublex=-0.99,y;intn=0;printf("%d,%lf\n",n,x);do{y=x;x=x-(x*x*x/3-x)/(x*x-1);n++;printf("n=%d,x=%lf\n",n,x);}while(fabs(y-x)>1e-5);}1.2.命令窗口結(jié)果截屏如下：2．用牛頓下山法求解2.1利用C語言編程求解，程序如下：#include<stdio.h>#include<math.h>doublef(doublex){ doublet;t=x*x*x/3-x;returnfabs(t);}voidmain(){doublex=-0.99,y;intk=0;printf("%d,%lf\n",k,x);do {y=x;intn=0;do{ x=y;x=x-pow(0.5,n)*(x*x*x/3-x)/(x*x-1); n++; }while(f(x)>=f(y));k++;printf("k=%d,n=%d,x=%lf\n",k,n,x);}while(fabs(x-y)>1e-5);}2.2.命令窗口結(jié)果截屏如下：實驗結(jié)果分析：從運算截圖中我們可以知道，一方面運用牛頓迭代法和牛頓下山迭代法的程序運算的結(jié)果是一致的，若出現(xiàn)小小的差異的話，可能是由于設(shè)置的變量類型不同導致，另一方面運用牛頓下山迭代法比運用一般牛頓迭代法計算的步驟要少一些，這證明對于同一個題，牛頓下山迭代法的優(yōu)越性比一般的牛頓迭代要高。實驗三解線性代數(shù)方程組——列主元消元法實驗內(nèi)容：用列主元法解線性方程組= 實驗過程：利用C語言編程求解，程序如下：#include<stdio.h>#include<math.h>voidmain(){doublet,m1,m2,m3,x1,x2,x3,x4,a[4][5]={1.1348,3.8326,1.1651,3.4017,9.5342,0.5301,1.7875,2.5330,1.5435,6.3941,3.4129,4.9317,8.7643,1.3142,18.4231,1.2371,4.9998,10.6721,0.0147,16.9237};inti,j;for(i=0;i<4;i++) {for(j=0;j<5;j++)printf("%lf",a[i][j]);printf("\n");}for(i=0;i<4;i++) for(j=0;j<5;j++) if(a[i][0]<a[i+1][0]) t=a[i][0]; a[i][0]=a[i+1][0];a[i+1][0]=t;for(i=1;i<4;i++){m1=a[i][0]/a[0][0];for(j=0;j<5;j++) a[i][j]=a[i][j]-m1*a[0][j];}printf("第一次消元a[][]=\n");for(i=0;i<4;i++) {for(j=0;j<5;j++) printf("%lf",a[i][j]);printf("\n");}printf("\n");for(i=1;i<4;i++) for(j=0;j<5;j++) if(a[i][0]<a[i+1][0]) t=a[i][0]; a[i][0]=a[i+1][0];a[i+1][0]=t;for(i=2;i<4;i++){m2=a[i][1]/a[1][1];for(j=0;j<5;j++) a[i][j]=a[i][j]-m2*a[1][j];}printf("第二次消元a[][]=\n");for(i=0;i<4;i++) {for(j=0;j<5;j++) printf("%lf",a[i][j]);printf("\n");}printf("\n"); for(i=0;i<4;i++)for(j=0;j<5;j++) if(a[i][0]<a[i+1][0]) t=a[i][0]; a[i][0]=a[i+1][0];a[i+1][0]=t;for(i=3;i<4;i++){m3=a[i][2]/a[2][2];for(j=0;j<5;j++) a[i][j]=a[i][j]-m3*a[2][j];}printf("第三次消元a[][]=\n");for(i=0;i<4;i++) {for(j=0;j<5;j++) printf("%lf",a[i][j]); printf("\n");} printf("\n"); x4=a[3][4]/a[3][3]; x3=(a[2][4]-x4*a[2][3])/a[2][2]; x2=(a[1][4]-x4*a[1][3]-x3*a[1][2])/a[1][1]; x1=(a[0][4]-x4*a[0][3]-x3*a[0][2]-x2*a[0][1])/a[0][0]; printf("x1=%lf\nx2=%lf\nx3=%lf\nx4=%lf",x1,x2,x3,x4);}2.命令窗口結(jié)果截屏如下：實驗結(jié)果分析：從程序中看，似乎程序很復雜，其實不然，中間有重復的語句，我這種方法雖然看起來比較過于冗長，但便于讀者閱讀，通俗易懂，我也曾想用循環(huán)體來代替中間的重復語句，但發(fā)現(xiàn)我編的循環(huán)體與程序中其他的循環(huán)體沖突，使程序無法運行出結(jié)果，怎么也理不清。既然計算機的運算能力很強，這點冗余度應該不是問題。實驗四解線性代數(shù)方程組——追趕法實驗內(nèi)容：，用追趕法解三角方程組Ax=d，A=，d=實驗過程：利用C語言編程求解，程序如下:#include<stdio.h>#include<math.h>voidmain(){inti;doublea[5]={0,-1,-1,-1,-1},b[5]={2,1,1,1,1},c[4]={2,2,2,2},f[5]={6,7,9,11,1};doubler[4],y[5],x[5];for(r[0]=c[0]/b[0],i=1;i<=3;i++)r[i]=c[i]/(b[i]-a[i]*r[i-1]);for(y[0]=f[0]/b[0],i=1;i<=4;i++)y[i]=(f[i]-a[i]*y[i-1])/(b[i]-a[i]*r[i-1]);for(x[4]=y[4],i=3;i>=0;i--)x[i]=y[i]-r[i]*x[i+1];for(i=0;i<=4;i++)printf("%lf\n",x[i]);}2.命令窗口結(jié)果截屏如下:實驗結(jié)果分析：運用追趕法與杜立特法有些類似，期間分解矩陣要運用到多個循環(huán)體，還有鑲嵌式的循環(huán)體，這樣很容易造成程序內(nèi)沖突，于是運用另外的方法，定義的數(shù)組不代表一行，而是代表矩陣對角線上有元素的對角斜線。這樣運算的過程就比較清晰了。實驗五插值與擬合——三次樣條差值實驗內(nèi)容：已知直升機旋轉(zhuǎn)機翼外形曲線輪廓線上的某些點及端點處的一階導數(shù)值，是S0=1.8648,S19=-0.046115試計算該曲線的橫坐標為2,4,6,12,16,30,60,110,180,280,400,515該點處的縱坐標（要求該曲線具有二階光滑度）。k01234560.523.18.017.9528.6539.6250.655.287949.413.8420.224.928.4431.1k78910111278104.6156.6208.6260.7312.53536.536.634.631.631.0k131415161718364.4416.346849450752020.914.87.83.71.50.2實驗過程：1.利用C語言編程求解，程序如下：#include<stdio.h>#include<math.h>voidmain(){inti;doublex[19]={0.52,3.1,8.0,17.95,28.65,39.62,50.65,78.0,104.6,156.6,208.6,260.7,312.5,364.4,416.3,468.0,494.0,507.0,520.0},y[19]={5.28794,9.4,13.84,20.2,24.9,28.44,31.1,35.0,36.5,36.6,34.6,31.6,31.0,20.9,14.8,7.8,3.7,1.5,0.2};doubleh[18],u[18],v[18],g[18],a[19],b[19],c[19];for(i=0;i<19;i++){h[i]=x[i+1]-x[i];u[i]=h[i]/(h[i]+h[i+1]); v[i]=1-u[i]; g[i+1]=(6/(h[i+1]+h[i+2]))*((y[i+2]-y[i+1])/h[i+2]-(y[i+1]-y[i])/h[i+1]);}g[0]=6/(h[1])*((y[1]-y[0])/h[1]-1.8648);g[18]=6/(h[19])*(-0.046115-(y[18]-y[17])/h[17]);for(i=0;i<19;i++)printf("x[%d]:%lf",i,x[i]);printf("y[%d]:%lf",i,y[i]);for(i=0;i<19;i++)b[i]=2;a[0]=0,a[18]=1;for(i=1;i<=17;i++)a[i]=u[i-1];c[0]=1;for(i=1;i<=17;i++)c[i]=v[i-1];doubleB[18],Y[19],M[19];for(B[0]=c[0]/b[0],i=1;i<=17;i++)B[i]=c[i]/(b[i]-a[i]*B[i-1]);for(Y[0]=g[0]/b[0],i=1;