数值分析上机实习报告

《数值分析》上机实习题指导老师：姓名：学号：专业：结构工程院系：土建学院上海交通大学2006级研究生第一题一程序说明：1Householder变换方法是：1958年由A.S.Householder提出来的，它的乘法运算次数仅是Givenr方法的一半，且只需要作次开方运算。归纳起来，对换矩阵三对角化的算法步骤为：令，，已知，即。，。。，。2用超松驰法求解BX=b(取松驰因子,=0,迭代9次)。3用列主元素消去法求解BX=b：n阶方程组的系数矩阵为：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变换方法可将对称正定矩阵三对角化，则松弛法最佳松弛因子可用公式求得，选择合适的松弛因子收敛速度很快。列主元素消去法算法稳定。第三题一程序说明：对给定区间[a,b]均匀划分π：a=x0<x1<L<xN-1<xN=b,x=a+ih,h=EQ(b-a)/N.由于三次样条函数空间是N+3维的,故把分点扩展到X--1,XN+1.由插值条件来确定cj.，由s(xi)=yi,(i=0,1,2,…..n),sxn)=y′e可得s(x)=s¹(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函数*/{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函数*/{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函数*/ sd=sd+x[i]*(fz2(t1+0.5)-fz2(t1-0.5));/*调用fz2函数*/}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’四问题讨论：样条插值效果较好,由于样条插值不必经过所有点，所以没有Runge现象.而且插值函数比较光滑。第四题一程序说明：牛顿迭代法又称牛顿切线法。求x7-28*x4+14=0在（0）中的近似根（初始值为区间端点，迭代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数，可由确定。如，，，，，，，可用，，序列去逼近积分值，现用Richardson外推法构造一新序列，使它更快地收敛于积分值，此时，以，代入上式得：(Ⅰ(Ⅰ)且由逼近积分值时，有如下估计：算法（Ⅰ）就称为数值积分的Romberg算法。我们可得具体的Romberg算法的计算步骤：先求出按梯形公式所得的积分值。把区间2等分，求出两个小梯形面积之和，记为，即：。这样由外推法可得，，这实际上是区间上的Simpson公式。把区间再等分（即等分），得到复化梯形公式，由与外推可得，，如此，若已算出等分的复化梯形公式，则由Richardson外推法，构造新序列最后求得。（4）或就停止计算，否则回到（3），计算，一般可用如下算法：其具体流程如下，并全部存入第一列（1）（3）（6）（2）（5）（4）二计算程序：#include"stdio.h"#include"math.h"#defineN10doubletx[N+1][N+1];/*定义tx函数*/doublehs(doublex)/*定义hs函数*/{doubley;y=pow(3,x)*pow(x,1.4)*(5*x+7)*sin(x*x);/*调用pow函数*/return(y);}intpow1(intm,intn)/*定义pow1函数*/{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函数*/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算法的优点是:把积分化为代数运算,而实际上只需求T1(i),以后用递推可得.算法简单且收敛速度快,一般4或5次即能达到要求.Romberge算法的缺点是: