數(shù)值分析-上機實驗(精簡版)

《數(shù)值分析》上機實習報告學院：專業(yè)：姓名：學號：指導教師：使用語言：C語言使用軟件：VC++6.02013年12月20日第一題：列主元素法求方程組根1.題目用列主元消去法求解Ax=b。2.理論依據(jù)和應(yīng)用條件:所謂列主元消去法是，對矩陣作恰當?shù)恼{(diào)整，選取絕對值最大的元素作為主元素。然后把矩陣化為上三角陣，再進行回代，求出方程的解。此題采用宏定義將其一般化了，對于類似的方程都可以運用此法求解。3.計算程序：#include<stdio.h>#include<math.h>#defineN9/*宏定義N=9*/voidmain(){ intn,t,k,i,j; doublem,h,x[N],A[N][N+1]=/*A和b合并成A[N][N+1]*/ {{12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.74238,3.067813,-2.031743,2.1874369}, {2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124,33.992318}, {-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103,-25.173417}, {1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.037585,0.84671695}, {-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784317,1.784317}, {0.71819,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.10356,0.238417,-86.612343}, {1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.7138465,3.123789,-2.213474,1.1101230}, {3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782,4.719345}, {-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001,-5.6784392}}; n=N;n=n-1;for(k=0;k<n;k++) {m=fabs(A[k][k]);t=k;for(i=k+1;i<=n;i++)if(fabs(A[i][k])>m)/*選取列主元m*/ { m=fabs(A[i][k]); t=i;/*記錄行號t*/ }for(j=k;j<=n+1;j++)/*換行*/ { h=A[k][j]; A[k][j]=A[t][j]; A[t][j]=h; }for(i=k+1;i<=n;i++)/*消元并計算*/for(j=k+1;j<=n+1;j++)A[i][j]=A[i][j]-A[i][k]*A[k][j]/A[k][k]; } x[n]=A[n][n+1]/A[n][n];/*回帶計算x*/for(k=n-1;k>=0;k--) { x[k]=A[k][n+1];for(j=k+1;j<=n;j++)x[k]=x[k]-A[k][j]*x[j];x[k]=x[k]/A[k][k]; }for(i=0;i<=n;i++)/*輸出x*/printf("x[%d]=%f;\n",i+1,x[i]);}4.實驗結(jié)果5.問題討論〔1〕首先，我們將A和b的數(shù)值定義到一個數(shù)組中了，這方便的程序的調(diào)用；〔2〕其次，程序的實現(xiàn)，關(guān)鍵在于主元的選取和換行，需要仔細的思考；〔3〕再次，本人將程序一般化了，類似方程組，都可以通過此算法實現(xiàn)。第二題：三次樣條插值求近似值題目函數(shù)值如下表：x12345f(x)00.6931471.09861231.38629441.6094378x678910f(x)1.79175951.94591012.0794452.19722462.3025851f'(x)f'(1)=1f'(10)=0.1試用三次樣條插值求f(4.563)和f'(4.563)的近似值2.理論依據(jù)和應(yīng)用條件:依據(jù)三彎矩插值法列出相應(yīng)的方程，然后將x的值帶入〔在求解三彎矩方程的參數(shù)時，應(yīng)用追趕法求解方程組求出相應(yīng)參數(shù)值〕。運用到課本上4.7.2公式以及S〔x〕的一階導數(shù)的公式，4.7.4公式，4.7.5公式以及追趕法的相關(guān)公式。3.計算程序：#include<stdio.h>#defineN10/*計算10個節(jié)點*/intx[N]={1,2,3,4,5,6,7,8,9,10},i,h=1;/*定義全局變量x[N],i,h*/doublef[N]={0,0.69314718,1.0986123,1.3862944,1.6094378, 1.7917595,1.9459101,2.079445,2.1972246,2.3025851};voidmain(){doubleM(intk),s,s1;/*調(diào)用doubleM(intk)函數(shù)*/ floatx0,x1,y; intj;printf("X=");scanf("%f",&y);j=int(y-1);/*因為數(shù)組的第一個值是x[0],不是x[1]*/ x0=(float)x[j]; x1=(float)x[j+1];s=(x1-y)*(x1-y)*(x1-y)/(6*h)*M(j)+(y-x0)*(y-x0)*(y-x0)/(6*h)*M(j+1)/*求S(X)*/ +(f[j]-M(j)*h*h/6)*(x1-y)/h+(f[j+1]-M(j+1)*h*h/6)*(y-x0)/h;printf("S=%f\n",s);s1=-(x1-y)*(x1-y)*M(j)/(2*h)+(y-x0)*(y-x0)*M(j+1)/(2*h)/*求S'(X)即一階導數(shù)*/ -(f[j]-M(j)*h*h/6)/h+(f[j+1]-M(j+1)*h*h/6)/h;printf("S1=%f\n",s1);}doubleM(intk)/*定義doubleM(intk)函數(shù)*/{ doubleu[N-1],v[N-1],b[N],l[N-1],r[N],z[N],n[N];/*v[N-1]即書上；b[N]即書上d；n[N]即書上x[N]*/doublefd0=1.0,fd9=0.1;floatd=2;for(i=0;i<=N-2;i++)/*追趕法求值*/if(i<N-2)/*求u[N-1]*/u[i]=(float)h/(h+h);elseu[i]=1;for(i=0;i<=N-2;i++)/*求v[N-1]*/if(i==0)v[i]=1;elsev[i]=1-u[i-1];for(i=0;i<=N-1;i++)/*求b[N]*/ if(i==0)b[i]=6*((f[i+1]-f[i])/h-fd0)/h;else {if(i==N-1)b[i]=6*(fd9-(f[i]-f[i-1])/h)/h;elseb[i]=6*((x[i]-x[i+1])*f[i-1]+(x[i+1]-x[i-1])*f[i]+(x[i-1] -x[i])*f[i+1])/((x[i-1]-x[i])*(x[i]-x[i+1])*(x[i-1]-x[i+1]));}for(i=0;i<=N-2;i++)/*求l[N-1]*/if(i==0)l[i]=u[i]/d;elsel[i]=u[i]/(d-l[i-1]*v[i-1]);for(i=0;i<N;i++)/*求r[N]*/if(i==0)r[i]=d;elser[i]=d-l[i-1]*v[i-1];for(i=0;i<N;i++)/*求z[N]*/if(i==0)z[i]=b[i];elsez[i]=b[i]-l[i-1]*z[i-1];for(i=N-1;0<=i;i--)/*求n[N]*/if(i==N-1)n[i]=z[i]/r[i];elsen[i]=(z[i]-v[i]*n[i+1])/r[i];return(n[k]);/*返回值即M(intk)*/}4.計算結(jié)果5.問題討論在做數(shù)組相關(guān)的編程時，注意數(shù)組下標是從0開始的；編程的過程中，要慎用全局變量；算法的設(shè)計，盡量簡潔；追趕法的使用過程中，注意對應(yīng)項的轉(zhuǎn)換和實現(xiàn)。第三題：牛頓迭代法求方程近似根題目用Newton法求方程在(0.1,1.9)中的近似根(初始近似值取為區(qū)間端點，迭6次或誤差小于0.00001).2.理論依據(jù)和應(yīng)用條件:根據(jù)Newton迭代法的2.3.1公式，本算法設(shè)計一般化了，類似方程，皆可以通過此算法實現(xiàn)。3.計算程序：#include<stdio.h>#include<math.h>voidmain(){floatx,x1,f,f1;x=1.9;do/*設(shè)置一個循環(huán)*/ {x1=x; f=(x1*x1*x1-28)*x1*x1*x1*x1+14; f1=(7*x1*x1*x1-112)*x1*x1*x1; x=x1-f/f1; }while(fabs(x-x1)>1e-6);/*判斷是否滿足精度*/ printf("Therootis：%.6f\n",x);}4.計算結(jié)果5.問題討論〔1〕在運用Newton法時，確定初值的選?。弧?〕編程的過程要用到數(shù)學函數(shù)時，必須輸入#include<math.h>。第四題：Romberg算法求積分1.題目用Romberg算法求.2.理論依據(jù)和應(yīng)用條件:運用梯形公式5.1.2和Romberg算法的5.4.1公式、5.4.4公式、5.4.5公式、5.4.6公式求解。對于同類積分，皆可以通過這個算法求得。3.計算程序#include<stdio.h>#include<math.h>#defineN10/*先取最大值10，假設(shè)不行直接改N的值即可*/doubleb=3,a=1;voidmain(){ doublet[N],s[N-1],c[N-2],r[N-3];doublef(doublex1),sum;intk,i;for(i=0;i<=N-1;i++){if(i==0)/*開始求T[N]，從0開始考慮*/ t[i]=(f(a)+f(b))*(b-a)/2.0;else{for(k=0;k<=pow(2,i-1)-1;k++)if(k==0) sum=f(a+(b-a)/(pow(2,i-1))/2);elsesum=sum+f(a+(b-a)/(pow(2,i-1))/2+k*(b-a)/(pow(2,i-1)));/*前i項的求和*/t[i]=(t[i-1]+(b-a)/(pow(2,i-1))*sum)/2;}if(i>0) s[i-1]=(t[i]*4.0-t[i-1])/3.0;/*求s[N-1]*/if(i>1) c[i-2]=16*s[i-1]/15-s[i-2]/15;/*求c[N-2]*/if(i>2){r[i-3]=(64*c[i-2]-c[i-3])/63.0;/*求r[N-3]*/if((fabs(r[i-3]-r[i-4]))<(1e-5))break;/*誤差小于1e-5時，結(jié)束*/}}printf("Theresultis%f\n",r[i-3]);/*輸出結(jié)果*/}doublef(doublex1)/*定義函數(shù)f(x)*/{doublet2;t2=pow(3,x1)*pow(x1,1.4)*(5*x1+7)*(sin(x1*x1));return(t2);}4.計算結(jié)果5.問題討論由于是順序求值，不會產(chǎn)生冗余運算，先允許N為10，假設(shè)沒結(jié)果；那么可以直接改變N的值即可；注意數(shù)組的下標是從0開始的，注意轉(zhuǎn)換；不同數(shù)據(jù)類型的數(shù)據(jù)運算時，注意轉(zhuǎn)換；算方法的設(shè)計過程中，注意算法的通用型，不要就題論題。第六題：定步長四階Runge-Kutta法求方程組題目用定步長四階Runge-Kutta求解h=0.0005，打印yi(0.025)，yi(0.045)，yi(0.085)，yi(0.1)，(i=1，2，3)2.理論依據(jù)和應(yīng)用條件:運用四階Runge-Kutta公式6.3.8，對于定步長的方程組，求解時，h是定值。由于是方程組，可以將看成一個向量。3.計算程序#include<stdio.h>#include<math.h>#defineN4/*宏定義N*/voidmain(){doubleh=0.0005,x=0,k[5][4],Y[4],y1[4],f[4];doubleF(doublex1,doubley[4],doublef[4]);/*調(diào)用函數(shù)F*/ doubleb[N]={0.025,0.045,0.085,0.1};inti,j,t,c;for(t=0;t<=N-1;t++)/*Runge-kutta開始*/ {c=int(b[t]/h);{for(j=0;j<=3;j++)Y[j]=0;/*初始條件為0*/for(i=1;i<=c;i++) {for(x=0;x<=b[t];x=x+h) for(j=1;j<=3;j++)/*求K1值*/ y1[j]=Y[j]; F(x,y1,f); for(j=1;j<=3;j++) k[1][j]=f[j];for(j=1;j<=3;j++)/*求K2值*/ y1[j]=Y[j]+0.5*h*k[1][j]; F(x,y1,f); for(j=1;j<=3;j++) k[2][j]=f[j]; for(j=1;j<=3;j++)/*求K3值*/ y1[j]=Y