数值分析第四次实习报告.doc

学生学号实验课成绩 武汉理工大学 学 生 实 验 报 告 书 实验课程名称： 数 值 分 析( 第四次实习) 开 课 名 称：计算机科学与技术学院 指导老师姓名： 熊盛武学 生 姓 名： 学生专业班级：软件工程0803班 2009 2010学年 第 一 学期一、计算实习的内容：实习1 给定矩阵A与向量bA= b= (1)求A的三角分解（不必输出）；(2)利用A的三角分解解下列方程组： Ax=b A2x=b A3x=b对第题分析一下，如果先求M=A2,再解Mxb，有何缺点？实习2 追赶法的优点分别对n=5,100,300解下列方程组AxbA= b=再用现有的三角分解法程序解此方程组，比较一下两者的效率，哪种方法机时少，可求解的矩阵阶数高？实习3 Jacobi 迭代法与Gauss-Seidel迭代法的收敛性与收敛速度研究用Jacobi迭代法与Gauss-Seide迭代法解下列方程组Ax=b的收敛性，通过上机计算，验证分析是否正确，并观察右端项对迭代收敛是否有影响，比较两法的收敛速度。（1）A= b1= b2= (2) A= b1= b2=(3) A= b=实习4 松弛因子对超松弛迭代法收敛速度的影响编写一个超松弛迭代法解方程组Axb的计算机程序，其中A= b= 分别对不同的阶数（例如n=10，n=100）w=1.1，1.2，1.9进行迭代，记录近似解e=时所用的迭代次数k，观察松弛因子对收敛速度的影响。实习5 用欧拉公式求解二、计算实习：（一）三角分解法程序源代码#includeusing namespace std;int main()const int MAX_N=20; static double aMAX_NMAX_N,bMAX_N,xMAX_N,yMAX_N;static double uMAX_NMAX_N,lMAX_NMAX_N;int i,j,k,n,r; coutn; if(nMAX_N) coutThe input n is larger than MAX_N,please redefine the MAX_Nendl; return 1; if(n=0) coutplease input n between 1 andMAX_N; return 1; /输入 aij,bi coutnow input the matrix aij,i,j=0.n-1 :; for(i=0;in;i+) for(j=0;jaij; coutnow input the matrix bi,i=0.n-1 :; for(i=0;ibi; /分解过程for(i=0;in;i+)lii=1;/l矩阵对角元素为1 for(k=0;kn;k+) for(j=0;jn-1;j+)/计算u矩阵 ukj=akj; for(r=0;r=k-1;r+)ukj-=lkr*urj; for(i=k+1;kn;i+)/计算l矩阵 lik=aik; for(r=0;r=k-1;r+) lik-=lir*urk; lik/=ukk; /ly=b for(i=0;in;i+) yi=bi; for(j=0;j=0;i-) xi=yi; for(j=i+1;jn;j+) xi-=uij*xj; xi/=uii; cout方程组的解为：;/输出x for(i=0;in;i+) coutxi ; return 0;（二）追赶法程序源代码#includeusing namespace std;int main()const int MAX_N=300;int aMAX_N,bMAX_N,cMAX_N,dMAX_N;float lMAX_N,rMAX_N,yMAX_N,xMAX_N;int i,n; coutn; if(nMAX_N) coutThe input n is larger than MAX_N,please redefine the MAX_Nendl; return 1; if(n=0) coutplease input n between 1 andMAX_N; return 1; /输入 an,bn,cn,dn coutnow input the matrix an,n=2.n :; for(i=2;iai;coutnow input the matrix bn,n=1.n :; for(i=1;ibi; coutnow input the matrix cn,n=1.n-1 :; for(i=1;ici; coutnow input the matrix dn,n=1.n :; for(i=1;idi; /追的过程 r0=0; y0=0; a1=0; cn=0; for(i=1;i=n;i+) li=bi-ai*ri-1; ri=ci/li; yi=(di-ai*yi-1)/li; /赶的过程 cout=1;i-) xi=yi-ri*xi+1; cout方程组的解为：;for(i=1;i=n;i+)coutxi ; coutendl; return 0;计算得解为：2 3 3 3 2（三）雅可比迭代法和高斯赛德尔迭代法程序源代码#include #include using namespace std; int a,b,m; double *x0; void Jacobi(double *c,double *d,int n,double eps); void Gauss(double *c,double *d,int n,double eps); void main() int n; double *A,*B; double e; cout请选择求方程组的迭代方法!雅可比选0，高斯-赛德尔选1!n; cout输入方程组的变量的个数以及方程的个数!ab; A=new double*b; for(int i=0;ib;i+) Ai=new doublea; B=new doubleb; x0=new doublea; cout输入每个方程组的变量的系数以及方程右端的值!endl; for(int k=0;kb;k+) for(int j=0;jAkj; cinBk; cout输入方程组迭代的次数及所要求的精度!me; cout输入方程组迭代的初值!endl; for(int j=0;jx0j; switch (n) case 0:Jacobi(A,B,m,e); break; case 1:Gauss(A,B,m,e); break; default:cout你没有选择求解方程组的一种方法!endl; break; void Jacobi(double *c,double *d,int n,double eps) int k,i; double *y = new doublea,*x=new doublea,s,temp=0.0; k=1; while(1) temp = 0.0; for(i=0;ia;i+) s=0.0; for(int j=0;jtemp) temp=fabs(x0i-s); if(tempeps) cout迭代成功!迭代结果为：endl; for(i=0;ia;i+) coutyi =yiendl; break; if(k=m) cout迭代失败!endl; break; k+=1; for(i=0;ia;i+) x0i=yi; void Gauss(double *c,double *d,int n,double eps) int k,i; double *y=new doublea,*x=new doublea,s,temp=0.0; for(i=0;ia;i+) xi=x0i; yi=xi; k=1; while(1) temp=0.0; for(i=0;ia;i+) s=0.0; for(int j=0;jtemp) temp=fabs(xi-s); if(tempeps) cout迭代成功!迭代结果为：endl; for(i=0;ia;i+) coutyi =yiendl; break; if(k=m) cout迭代失败!endl; break; k+=1; for(i=0;ia;i+) xi=yi; 用雅可比迭代法，程序运行结果如下:用高斯赛德尔迭代法程序运行结果如下:通过运行结果可以看出高斯赛德尔迭代法比雅可比迭代法所求结果精度高，收敛性好。（四）超松弛迭代法程序源代码#includeusing namespace std;#include#define MAX_N 20#define MAXREPT 100#define epsilon 0.00001int main() int n;int i,j,k;double err,w;static double aMAX_NMAX_N,bMAX_NMAX_N,cMAX_N,gMAX_N;static double xMAX_N,nxMAX_N;coutn; if(nMAX_N) coutThe input n is larger than MAX_N,please redefine the MAX_Nendl; return 1; if(n=0) coutplease input n between 1 andMAX_N; return 1; /输入aij,ci coutnow input the matrix aij,i,j=0.n-1 :; for(i=0;in;i+) for(j=0;jaij; coutnow input the matrix ci,i=0.n-1 :; for(i=0;ici; coutw; if(w=2) coutw must between 1 and 2.endl; return 1; for(i=0;in;i+)/形成x_k+1=bx_k+g迭代矩阵bfor(j=0;jn;j+)bij=-aij/aii;gi=ci/aii;/为了简化程序，假设aii!=0/否则要附加对aii的处理for(i=0;jMAXREPT;i+)for(j=0;jn;j+)nxj=gj;for(j=0;jn;j+)for(k=0;kn;k+)if(j=k)continue;nxj+=bjk*nxk;/迭代for(k=j+1;kn;k+)nxj+=bjk*xk;nxj=(1-w)*xj+w*nxj;err=0;for(j=0;jn;j+)if(errfabs(nxj-xj)err=fabs(nxj-xj);for(j=0;jn;j+)xj=nxj;if(errepsilon)coutthe solve are:;for(i=0;in;i+)coutxi ;coutendl;return 0;coutafter MAXREPTrepeat,no result.;/输出return 1;1.2时 解为：0232179 0.0662066 0.4941941.4时 解为：023218 0.0662066 0.4941941.6时 解为：023218 0.0662074 0.494193由以上3组结果可知，越小，精度越大（五）欧拉公式程序源代码#include#include#define f(x,y) (x+y)int main()int m;int i;double a,b,y0;double xn,yn,xn1,yn1,yn1b;double h;couta; cinb; couty0; coutm; if(m=0) coutplease input a number lager than 1.; return 1; h=(b-a)/m; xn=a;yn=y0; for(i=1;i=m;i+) xn1=xn+h; yn1b=yn+h*f(xn,yn); yn1=yn+h/2*(f(xn,yn)+f(xn1,yn1b); coutxi xn1 yi yn1; xn=xn1;yn=yn1; return 0;调试结果： x10.2 y11.24 x2=0.4 y2=1.5768x30.6 y32.031x40.8 y42.63067x51 y53.40542第四次实验小结及体会1、任何非奇异矩阵都可分解为三角形式。2、如果需要重复