數(shù)值分析上機實習報告

《數(shù)值分析》上機實習題指導老師：姓名：學號：專業(yè)：結構工程院系：土建學院上海交通大學2006級研究生第一題一程序說明：1Householder變換方法是：1958年由A.S.Householder提出來的，它的乘法運算次數(shù)僅是Givenr方法的一半，且只需要作次開方運算。歸納起來，對換矩陣三對角化的算法步驟為：令，，已知，即。，。。，。2用超松馳法求解BX=b(取松馳因子,=0,迭代9次)。3用列主元素消去法求解BX=b：n階方程組的系數(shù)矩陣為：a11a12…a1nb1a21a22…a2nb………………an1an2…annbnGauss消去法的算法為：lij=aij/ajj(ajj!=0)j=1,2…n,i=j+1,j+2…n⑴aij=aij-lik-1ak-1i,j=k,k+1…n,k=2,3…nbi=bi-lik-1bk-1i=k,k+1…n,k=2,3…n⑵xi=(bi-∑aijxj)/aiii=n,n-1…1,j=i+1,i+2,…n列主元素消去法是在Gauss消去法的基礎上選主元，選取絕對值最大（或盡量大）的元素為主元，使lij絕對值很小。二計算程序:1househloder變換將A化為三對角陣B#include<stdio.h>#include<math.h>main(){inti,j,m,r,sign;doubleA0[9][9],s,z,p,n,v,h,l,y[9],u[9],k,q[9],X[9],x[9]={0,0,0,0,0,0,0,0,0},B[9][9],g[9];doubleA[9][9]={{12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743},{2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124},{-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103},{1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.037585},{-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784317},{0.718719,1.563849,1.836742,-3.741865,0.431637,9.789365,-0.103458,-1.103456,0.238417},{1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.713846,3.123789,-2.213474},{3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782},{-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001}};doubleb[9]={2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392};for(i=0;i<9;i++)/*給A0賦初值*/for(j=0;j<9;j++)A0[i][j]=A[i][j];for(r=0;r<7;r++)/*求*/{s=0;for(i=r+1;i<9;i++)s=s+A[i][r]*A[i][r];s=sqrt(s);l=s*s+fabs(A[r+1][r])*s;if(A[r+1][r]>0)sign=1;elseif(A[r+1][r]<0)sign=-1;for(i=0;i<9;i++){if(i<=r)u[i]=0; elseif(i==r+1)u[i]=A[r+1][r]+sign*s; elseu[i]=A[i][r];}for(i=0;i<9;i++)/*求*/{y[i]=0;for(j=0;j<9;j++)y[i]=y[i]+A[i][j]*u[j]; y[i]=y[i]/l;}k=0;for(i=0;i<9;i++)/*求*/{k=k+u[i]*y[i];}k=k/(2*l);for(i=0;i<9;i++)/*求*/{q[i]=y[i]-k*u[i];}for(i=0;i<9;i++)/*求*/for(j=0;j<9;j++)A[i][j]=A[i][j]-q[i]*u[j]-q[j]*u[i];}printf("\nHousehold變換矩陣B為:\n");for(i=0;i<9;i++){printf("\n"); for(j=0;j<9;j++) printf("%f,",A[i][j]);} }2超松弛法求解Bx=b#include"stdio.h"#include"math.h"#defineN9floatA[N][N]={{12.384120,-4.893077,0,0,0,0,0,0,0},{-4.893077,25.398417,6.494088,0,0,0,0,0,0}, {0,6.494088,20.611538,8.243930,0,0,0,0,0},{0,0,8.243930,23.422894,-13.880069,0,0,0,0}, {0,0,0,-13.880069,29.698231,4.534562,0,0,0},{0,0,0,0,4.534562,16.006021,4.881358,0,0}, {0,0,0,0,0,4.881358,26.013397,-4.503611,0},{0,0,0,0,0,0,-4.503611,21.254171,4.504458},{0,0,0,0,0,0,0,4.504458,14.534007}};floatB[N][N];floatx[N]={0,0,0,0,0,0,0,0,0};floatg[N]={2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392};main(){inti,k,j;floatn;for(i=0;i<N;i++)/*求gi=bi/aii*/if(A[i][i]!=0)g[i]=g[i]/A[i][i];for(i=0;i<N;i++)/*求bij=-aij/aii*/for(j=0;j<N;j++) {if(i==j) if(A[i][j]==0) B[i][j]=1; else B[i][j]=0; else if(A[i][i]==0) B[i][j]=-A[i][j]; else B[i][j]=-A[i][j]/A[i][i];}for(i=0;i<N;i++)for(j=0;j<N;j++) {n=0; if(j==0) for(k=1;k<N;k++) n=B[j][k]*x[k]+n; else {for(k=0;k<j;k++) n=B[j][k]*x[k]+n; for(k=j+1;k<N;k++) n=B[j][k]*x[k]+n;} x[j]=(1-w)*x[j]+w*(n+g[j]);/*求xi*/}for(i=0;i<N;i++)printf("%.2f",x[i]);printf("\n");}3用列主元素消去法求Bx=b#include"stdio.h"#include"math.h"#defineN9floatB[N][N]={{12.38412,-4.893077,0,0,0,0,0,0,0},{-4.893077,25.398417,6.494088,0,0,0,0,0,0},8,8.243930,0,0,0,0,0},{0,0,8.243930,23.422894,-13.880069,0,0,0,0},{0,0,0,-13.880069,29.698231,4.534562,0,0,0},{0,0,0,0,4.534562,16.006021,4.881358,0,0},{0,0,0,0,0,4.881358,26.013397,-4.503611,0},{0,0,0,0,0,0,-4.503611,21.254171,4.504458},{0,0,0,0,0,0,0,4.504458,14.534007}};floatx[N],l[N];floatg[N]={2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392};main(){inti,k,n,j;for(i=0;i<N;i++){floatmax=0,t,h; for(j=i;j<N;j++)/*找出每列最大值*/ if(abs(B[j][i])>=max) {max=abs(B[i][j]); n=j;} for(j=i;j<N;j++)/*交換B[i][j]與B[n][j]使主對角線上值最大*/ {t=B[i][j];B[i][j]=B[n][j];B[n][j]=t;} h=g[i];g[i]=g[n];g[n]=h;/*矩陣三角化*/ for(j=i+1;j<N;j++) l[j]=B[j][i]/B[i][i]; for(j=i+1;j<N;j++) if(B[j][i]!=0) {for(k=i;k<N;k++) B[j][k]=B[j][k]-l[j]*B[i][k]; g[j]=g[j]-l[j]*g[i];}}for(i=N-1;i>=0;i--)/*計算解值*/if(i==N-1) x[i]=g[i]/B[i][i];else {x[i]=g[i]/B[i][i]; for(j=N-1;j>i;j--) x[i]=x[i]-(B[i][j]*x[j])/B[i][i]; }for(k=0;k<N;k++)/*打印解值*/printf("%.2f",x[k]);printf("\n");}三計算結果:1．Householder變換三對角陣B:12.3841202.超松弛法求解Bx=b:1.072.27–2.862.292.11–6.421.360.63–3.用列主元素消去法求Bx=b:1.082.28–2.862.292.11–6.421.360.63–四問題討論：Householder變換方法可將對稱正定矩陣三對角化，則松弛法最佳松弛因子可用公式求得，選擇合適的松弛因子收斂速度很快。列主元素消去法算法穩(wěn)定。第三題一程序說明：對給定區(qū)間[a,b]均勻劃分π：a=x0<x1<L<xN-1<xN=b,x=a+ih,h=EQ(b-a)/N.由于三次樣條函數(shù)空間是N+3維的,故把分點擴展到X--1,XN+1.由插值條件來確定cj.，由s(xi)=yi,(i=0,1,2,…..n),sxn)=y′e可得s(x)=s1(x)=矩陣形式為：求解線性方程組可得Ci,從而得S(x).二計算程序:#include<stdio.h>#include<math.h>#defineN12floatB[N][N]={{-1,0,1,0,0,0,0,0,0,0,0,0},{1,4,1,0,0,0,0,0,0,0,0,0}, {0,1,4,1,0,0,0,0,0,0,0,0},{0,0,1,4,1,0,0,0,0,0,0,0}, {0,0,0,1,4,1,0,0,0,0,0,0},{0,0,0,0,1,4,1,0,0,0,0,0}, {0,0,0,0,0,1,4,1,0,0,0,0},{0,0,0,0,0,0,1,4,1,0,0,0}, {0,0,0,0,0,0,0,1,4,1,0,0},{0,0,0,0,0,0,0,0,1,4,1,0}, {0,0,0,0,0,0,0,0,0,1,4,1},{0,0,0,0,0,0,0,0,0,-1,0,1}};floatx[N],l[N],k1[N]={0,1,2,3,4,5,6,7,8,9,10,11};2944,6*1.6094378,6*1.7917595, 6*1.9459101,6*2.079445,6*2.1972246,6*2.3025851,2*0.1};floatfz1(floatk)/*定義fz1函數(shù)*/{floath;if(fabs(k)>=2)h=0;elseif(fabs(k)<=1)h=(k*k*fabs(k))/2-k*k+2.0/3;elseh=-(k*k*fabs(k))/6+k*k-2*fabs(k)+4.0/3;return(h);}floatfz2(floatk)/*定義fz2函數(shù)*/{floath;if(fabs(k)>=1.5)h=0;elseif(fabs(k)<=0.5)h=-k*k+3.0/4;elseh=(k*k)/2-(3*fabs(k))/2+9.0/8;return(h);}main(){inti,k,n,j;floats=0,sd=0,t1;for(i=0;i<N-1;i++)/*找出每列最大值*/{floatmax=0,t,h; for(j=i;j<N;j++) if(fabs(B[j][i])>max) {max=fabs(B[i][j]); n=j;} for(j=i;j<N;j++)/*交換B[i][j]與B[n][j]使主對角線上值最大*/ {t=B[i][j]; B[i][j]=B[n][j]; B[n][j]=t;} h=g[i];g[i]=g[n];g[n]=h;/*矩陣三角化*/ for(j=i+1;j<N;j++) l[j]=B[j][i]/B[i][i]; for(j=i+1;j<N;j++) { if(B[j][i]!=0) for(k=i;k<N;k++) B[j][k]=B[j][k]-l[j]*B[i][k]; g[j]=g[j]-l[j]*g[i]; }}for(i=N-1;i>=0;i--)/*計算C值*/if(i==N-1) x[i]=g[i]/B[i][i];else { x[i]=g[i]/B[i][i]; for(j=N-1;j>i;j--) x[i]=x[i]-(B[i][j]*x[j])/B[i][i]; }for(k=0;k<N;k++)/*打印C值*/printf("%.2f",x[k]);printf("\n");for(i=0;i<N;i++){ t1=4.563-k1[i]; s=s+x[i]*fz1(t1);/*調用fz1函數(shù)*/ sd=sd+x[i]*(fz2(t1+0.5)-fz2(t1-0.5));/*調用fz2函數(shù)*/}printf("s(4.563)=%.4f\n",s);printf("s'(4.563)=%.4f\n",sd);}三計算結果：-1.270.140.731.121.401.621.801.952.082.202.302.40s’四問題討論：樣條插值效果較好,由于樣條插值不必經(jīng)過所有點，所以沒有Runge現(xiàn)象.而且插值函數(shù)比較光滑。第四題一程序說明：牛頓迭代法又稱牛頓切線法。求x7-28*x4+14=0在（0）中的近似根（初始值為區(qū)間端點，迭代6次或誤差小于0.00001）用牛頓迭代法求7次方程的根，當N=1時，有：所以設其初始值為1.9，牛頓迭代公式為：二計算程序：#include<math.h>main(){doublea,b,x,y,k;x=0.1;y=1.9;k=1;while(fabs(k)>0.00001){x=y;a=x*x*x*x*x*x*x-28*x*x*x*x+14;b=7*x*x*x*x*x*x-28*4*x*x*x;y=x-a/b;k=x-y;}printf("%f\n",x);}三計算結果：0．845506四問題討論：Newton平方收斂，它是局部收斂。第五題一程序說明：Romberg算法是：把外推算法直接應用到復化梯形公式，可以證明，復化梯形公式的誤差的Euler-Maclaurin公式：這里，是Bernoulli數(shù)，可由確定。如，，，，，，，可用，，序列去逼近積分值，現(xiàn)用Richardson外推法構造一新序列，使它更快地收斂于積分值，此時，以，代入上式得：(Ⅰ(Ⅰ)且由逼近積分值時，有如下估計：算法（Ⅰ）就稱為數(shù)值積分的Romberg算法。我們可得具體的Romberg算法的計算步驟：先求出按梯形公式所得的積分值。把區(qū)間2等分，求出兩個小梯形面積之和，記為，即：。這樣由外推法可得，，這實際上是區(qū)間上的Simpson公式。把區(qū)間再等分（即等分），得到復化梯形公式，由與外推可得，，如此，若已算出等分的復化梯形公式，則由Richardson外推法，構造新序列最后求得。（4）或就停止計算，否則回到（3），計算，一般可用如下算法：其具體流程如下，并全部存入第一列（1）（3）（6）（2）（5）（4）二計算程序：#include"stdio.h"#include"math.h"#defineN10doubletx[N+1][N+1];/*定義tx函數(shù)*/doublehs(doublex)/*定義hs函數(shù)*/{doubley;y=pow(3,x)*pow(x,1.4)*(5*x+7)*sin(x*x);/*調用pow函數(shù)*/return(y);}intpow1(intm,intn)/*定義pow1函數(shù)*/{inti,h=1;/*求解mn的值*/if(n==0)h=1;elsefor(i=1;i<=n;i++) h=h*m;return(h);}main(){intl,i,j,k,kk;charch;inta=1,b=3,L;doublen,h,w,p;tx[1][0]=((b-a)*(hs(a)+hs(b)))/2;/*調用pow函數(shù)*/for(l=1;l<=N;l++)/*求的值*/{L=pow1(2,l-1); n=0; for(i=1;i<=L;i++) n=n+hs(a+((2*i-1)*(b-a))/pow(2,l)); w=(tx[1][l-1]+((b-a)*n)/pow(2,l-1))/2; tx[1][l]=w; for(k=l-1;k>=0;k--)/**/ {p=(pow(4,l-k)*tx[l-k][k+1]-tx[l-k][k])/(pow(4,l-k)-1); tx[1+l-k][k]=p;} if(tx[l][0]-tx[l+1][0]<0)/*求誤差*/ h=-(tx[l][0]-tx[l+1][0]); else h=tx[l][0]-tx[l+1][0]; if(h<0.00001)/*確定跳出求解的條件*/ {kk=l;break;}}printf("Rombergfa:%.6f\n",tx[kk+1][0]);}三計算結果：四問題討論：Romberge算法的優(yōu)點是:把積分化為代數(shù)運算,而實際上只需求T1(i),以后用遞推可得.算法簡單且收斂速度快,一般4或5次即能達到要求.Romberge算法的缺點是: