数值代数上机报告.doc

Doolittle分解报告1一、 目的意义： 把矩阵A分解成一个下三角阵与一个上三角阵的乘积，即 A=LR，其中L为下三角阵，R为上三角阵，这样原线性方程组就可以化为的求解问题，方便求解。二、 算法：1) 输入系数矩阵 A；2) 利用公式和交错进行，计算得出矩阵L和R；3) 回带到中得出原线性方程组的解。三、 源程序：#include #include #include #include #define N 100main()int i,j,k,s,n;printf(请输入系数矩阵A的阶数：n= );scanf(%d,&n);float aNN=0,LNN=0,RNN=0,sigma1,sigma2,bN,yN,xN;/*为L主对角线元素赋1*/for(i=0;in;i+)Lii=1;printf(请输入系数矩阵A：n); /*输入系数矩阵A*/for(i=0;in;i+)for(j=0;jn;j+)scanf(%f,&aij);for(k=0;kn;k+)for(j=k;jn;j+) /*计算矩阵R*/sigma1=0;for(s=0;s=k-1;s+)sigma1+=Lks*Rsj;Rkj=akj-sigma1;for(i=k;in;i+) /*计算矩阵L*/sigma2=0;for(s=0;s=k-1;s+)sigma2+=Lis*Rsk;Lik=(aik-sigma2)/Rkk;printf(n A矩阵为：n);/*输出矩阵L、R*/for(i=0;in;i+)for(j=0;jn;j+)printf(%5.1f ,aij);printf(n);printf(n L矩阵为：n);for(i=0;in;i+)for(j=0;jn;j+)printf(%5.1f ,Lij);printf(n);printf(n R矩阵为：n);for(i=0;in;i+)for(j=0;jn;j+)printf(%5.1f ,Rij);printf(n);printf(请输入b矩阵n,i+1);for(i=0;in;i+)scanf(%f,&bi);for(i=0;in;i+)/*回代法求解方程组Ly=b*/sigma1=0;for(k=0;k=0;i-)sigma2=0;for(k=i+1;kn;k+)sigma2+=Rik*xk;xi=(yi-sigma2)/Rii;printf(解得x为：n);for(i=0;in;i+)printf(%5.1f ,xi);printf(n);四、 计算结果与分析：分析：运行结果与预想的结果相近，误差对结果的影响不是很大，比较理想五、 参考文献：1刑志栋. 矩阵数值分析. 陕西： 陕西科学技术出版社， 20052谭浩强. C语言程序设计. 北京：清华大学出版社，2005报告2cholesky分解一、 目的意义： 对称正定矩阵是实践中经常遇到的一种特殊矩阵类型矩阵，由于矩阵本身的流量好兴致，使得cholesky分解在存储和运算量上较一般消去法节省一半左右，且解的精度高，cholesky分解方法是目前计算机求解该类问题最有效的方法之一。二、 算法：1) 输入系数矩阵 A；2) 利用公式 交错进行，计算得出矩阵L和D；3) 回带到中得出原线性方程组的解三、 源程序：#include #include #include #define EPS 1.0e-8#define N 20double aNN, bN, xN;int n;int zhuyuan(int row); /* 选主元*/void hangjiaohuan(int row1, int row2); /* 行交换*/void xiaoyuan(int row); /*消元*/void huidai(); /* 回代*/void main()printf(请输入方程的维数n: n = ); scanf(%d, &n); getchar(); printf(输入%d行%d列系数矩阵A：n, n, n); for (int i=0; in; i+)for (int j=0; jn; j+)scanf(%lf, &aij);getchar();printf(n输入线性方程组右端项b%d:n, n); for (i=0; in; i+)scanf(%lf, &bi);getchar(); for (i=0; in-1; i+) double rowmark = zhuyuan(i);if (rowmark = -1) printf(多解!); system(pause); return;if (rowmark != i) hangjiaohuan(i, rowmark);xiaoyuan(i);huidai(); printf(n线性方程组的解为:n ); for (i=0; in; i+) printf(x%d=%lf , i+1, xi); printf(n); printf(n); system(pause);int zhuyuan(int row)double elem = arowrow; int rowmark = row; for (int i=row+1; in; i+) if (elemairow) elem = airow; rowmark = i;if(fabs(elem) = EPS) return -1;return rowmark;void hangjiaohuan(int row1, int row2)double tmp; tmp = brow1; brow1 = brow2; brow2 = tmp; for (int j=0; jn; j+) tmp = arow1j; arow1j = arow2j; arow2j = tmp;void xiaoyuan(int row)for (int i=row+1; in; i+)int j=row; double tmp = aij/arowrow; bi -= tmp*brow; for (aij+ = 0; j=0; i-)double sum = 0.0; for (int j=i+1; jn; j+) sum -= aij*xj;xi = (bi+sum)/aii;四、 计算结果与分析：分析：运行结果与预想的结果相近，误差对结果的影响不是很大，比较理想五、 参考文献：1刑志栋. 矩阵数值分析. 陕西： 陕西科学技术出版社， 20052谭浩强. C语言程序设计. 北京：清华大学出版社，2005Jacobi迭代法报告3一、 目的意义：通过有限次的算术运算得到线性方程组的近似值二、 算法：1) 输入系数矩阵 A；2) 输入迭代的初始向量；3) 利用公式 ，计算得出原线性方程组的近似解。三、 源程序：#include#include#include#include#define EPS 1e-6#define MAX 100#define N 100float *Jacobi(float aNN,int n)float *x,*y,s; double epsilon; int i,j,k=0; x=(float *)malloc(n*sizeof(float); y=(float *)malloc(n*sizeof(float); printf(请输入迭代的初始向量：n); for(i=0;in;i+) scanf(%f,&xi); while(1) epsilon=0; k+; for(i=0;in;i+) s=0; for(j=0;jn;j+) if(j=i) continue; s+=aij*xj; yi=(ain-s)/aii; epsilon += fabs(yi-xi); if(epsilon =MAX) printf(The Method is disconvergent); return y; for(i=0;in;i+) xi=yi;main()int i,j,n; float aNN; float *x;printf(请输入系数矩阵的阶数：n=);scanf(%d,&n);printf(请按行输入线性方程组的增广矩阵：n); for(i=0;in;i+)for(j=0;jn+1;j+)scanf(%f,&aij);x=(float *)malloc(3*sizeof(float);x=Jacobi(a,n);printf(n原线性方程组的解为：n); for(i=0;in;i+)printf(x%d=%fn,i,xi);getch();四、 计算结果与分析：分析：运行结果与预想的结果相近，误差对结果的影响不是很大，比较理想五、 参考文献：1刑志栋. 矩阵数值分析. 陕西： 陕西科学技术出版社， 20052谭浩强. C语言程序设计. 北京：清华大学出版社，2005报告4乘幂法一、 目的意义： 阶数超过四次的多项式的根一般不能用有限次运算求得，因而矩阵特征值的计算方法本质上总是采用迭代格式，而采用乘幂法可以计算最大的一个或按模最大的几个特征值和特征相应特征向量。二、 算法：1) 输入系数矩阵 A；2) 取出事向量为；3) 利用乘幂法的迭代格式 计算得出原矩阵按模最大特征值和相应特征的向量。三、 源程序：#include#include#define N 100#define M 100#define epsilon 0.00001int main()int n;int i,j,k;double xmax,oxmax;static double aNN;static double xN,nxN;printf(输入矩阵维数n: );scanf(%d,&n);if(nN)printf(the inpur n is larger than MAX_N,please redefine the N.n);return 1;if(n=0)printf(please input a number between 1 and %d.n,N);return 1;/输入A矩阵printf(请输入矩阵A(%d,%d):n,n,n);for(i=0;in;i+)for(j=0;jn;j+)scanf(%lf,&aij);for(i=0;in;i+) xi=1;oxmax=0;for(i=0;iM;i+)for(j=0;jn;j+)nxj=0;for(k=0;kn;k+)nxj+=ajk*xk;xmax=0.0;for(j=0;jxmax)xmax=fabs(nxj);for(j=0;jn;j+)nxj/=xmax;for(j=0;jn;j+)xj=nxj;if(fabs(xmax-oxmax)epsilon)printf(按模最大特征值=%.4fn,xmax);printf(相应的特征向量为:n);for(i=0;in;i+)printf(%.4fn,nxi);break;/return 0;oxmax=xmax;return 0;四、 计算结果与分析：分析：运行结果与预想的结果相近，误差对结果的影响不是很大，比较理想五、 参考文献：1刑志栋. 矩阵数值分析. 陕西： 陕西科学技术出版社， 20052谭浩强. C语言程序设计. 北京：清华大学出版社，2005报告5逆幂法一、 目的意义： 逆幂法用于计算非奇异矩阵A的按模最小特征值和相应的特征向量二、 算法：1) 输入系数矩阵 A和近似特征值；2) 将乘幂法用于便可求出A的按模最小的特征值和相应的特征向量。三、 源程序：#include#include#define e 0.0001#define N 100void main()int i,j,k,n,p,s,kk=1,flag=1; double q,m,lm,Lm,m0; double ANN,BNN,LNN,RNN,lNN,rNN;double z0N,z1N,y0N,y1N,yN;printf(请输入矩阵的阶数：n=);scanf(%d,&n); printf(n请输入矩阵A:n); for(i=1;i=n;i+)for(j=1;j=n;j+)scanf(%lf,&Aij);printf(n请输入近似特征值：); scanf(%lf,&lm); /* 令 B=A-lm*I ,用Doolittle分解B=LR */for(i=1;i=n;i+)for(j=1;j=n;j+)if(j=i)Bij=Aij-lm;else Bij=Aij; for(k=1;k=n;k+)for(i=k;i=n;i+)for(q=0,p=1;pk;p+)q=q+Lkp*Rpi;Rki=Bki-q;for(i=k+1;i=n;i+) for(q=0,p=1;pk;p+) q=q+Lip*Rpk; Lik=(Bik-q)/Rkk;for(i=1;i=n;i+)for(j=i;j=n;j+)if(i=j)Lij=1; else Lij=0;for(i=2;i=n;i+)for(j=1;ji;j+)Rij=0; /*第一步迭代，用Doolittle计算R*y1=y0中的y1*/for(i=1;i=n;i+)y0i=1;for(k=1;k=n;k+)for(i=k;i=n;i+)for(q=0,p=1;pk;p+)q=q+lkp*rpi;rki=Rki-q;for(i=k+1;i=n;i+)for(q=0,p=1;pk;p+) q=q+lip*rpk; lik=(Rik-q)/rkk;for(i=1;i=n;i+)for(q=0,k=1;k=1;i-)for(q=0,k=i+1;k=n;k+) q=q+rik*y1k; y1i=(yi-q)/rii; /*计算m,z1(存放在下列代码的z0中)*/for(m=y11,i=2;ifabs(y1i-1)m=y1i;for(i=1;i1)次迭代*/ kk=2; while(flag=1) /*用Doolittle解方程 LR*y2=z1中的y2，其中B=LR */ for(k=1;k=n;k+) for(i=k;i=n;i+)for(q=0,p=1;pk;p+)q=q+lkp*rpi;rki=Bki-q; for(i=k+1;i=n;i+) for(q=0,p=1;pk;p+) q=q+lip*rpk; lik=(Bik-q)/rkk; for(i=1;i=n;i+) for(q=0,k=1;k=1;i-) for(q=0,k=i+1;k=n;k+) q=q+rik*y1k; y1i=(yi-q)/rii; /*规范化*/ for(m=y11,i=2;ifabs(y1i-1)m=y1i; for(i=1;i=e，则进行下一轮迭代 ；否则得解。*/ for(s=0,i=1;i=n;i+) if(fabs(z0i)-fabs(z1i)e)s+; if(s=n)flag=0; else kk+; for(i=1;i=n;i+) y0i=y1i; z0i=z1i; m0=1/m; Lm=lm+m0; printf(最接近的特征值为：n); printf(%fnn,Lm); printf(相应的特征向量为：n); for(i=1;i=n;i+) printf(%lf,z1i); printf(nn); 四、 计算结果与分析：分析：运行结果与预想的结果相近，误差对结果的影响不是很大，比较理想五、 参考文献：1刑志栋. 矩阵数值分析. 陕西： 陕西科学技术出版社， 20052谭浩强. C语言程序设计. 北京：清华大学出版社，2005报告6古典Jacobi算法一、 目的意义： 寻找相似变换，是对称矩阵A经过变换之后所得矩阵的非对角线元素的平方和减少，对角线元素的平方和增大，且保持对称性不变，不断的施行这种正交变换，最终是非对角元素的平方和任意的接近与零，对角线元素平方的取极大值。二、 算法：1) 输入系数矩阵 A；2) 取，形成一个相似矩阵序列，k=1,2,3) 选取，使得元素=0，这时应满足 4) 令，因此可以得到： ，三、 源程序：#include stdio.h#includemath.h#define P 0.1void main()double a202020,b202020,r202020,l202020; double max,u,y,x,t,g,f,h,v,w; int i,j,k,n,c,d; printf(请输入矩阵的维数：n=); scanf(%d,&n); printf(请输入矩阵A:n); for(i=1;i=n;i+) for(j=1;j=n;j+)scanf(%lf,&aij1);printf(n所有的特征值为：); for(i=1;i=n;i+) for(j=1;j=n;j+) if(aij10) bij1=-aij1; else bij1=aij1; for(k=1;k=n;k+) for( max=0,i=1;i=n;i+)for(j=1;j=n;j+)if(imax+0.000001)max=bijk; c=i; d=j;u=bcck-bddk;if(u=0) y=u; else y=-u; if(u=0) x=2*bcdk; else x=-2*bcdk; t=sqrt(x*x+y*y); g=sqrt(0.5*1+0.5*y/t); f=0.5*x/t*(1/g);for(i=1;i=n;i+)for(j=1;j=n;j+) if(i!=j&i!=c,j!=d&i!=d,j!=c) rijk=0; if(i=c&j=c) rcck=g; if(i=d&j=d) rddk=g; if(i=c&j=d) rcdk=f; if(i=d&j=c) rdck=-f; if(i=j&i!=c&j!=d) rijk=1;printf(n);for(i=1;i=n;i+)for(j=1;j=n;j+)aick+1=bick*g+bidk*f; acik+1=bick*g+bidk*f; aidk+1=-bick*f+bidk*g; adik+1=-bick*f+bidk*g; if(i!=c&j!=d&i!=d&j!=c) aijk+1=bijk;h=(bcck-bddk)*f*g; v=bcdk*(g*g-f*f); acdk+1=h+v; adck+1=h+v; acck+1=bcck*g*g+2*bcdk*f*g+bddk*f*f; addk+1=bcck*f*f-2*bcdk*f*g+bddk*g*g;for(i=1;i=n;i+)for(j=1;j=n;j+)if(aij10) bijk+1=-aijk+1; else bijk+1=aijk+1;w=0;for(i=1;i=n;i+)for(j=1;j=n;j+)if(i!=j) w=w+bijk+1*bijk+1; else c