数值代数上机报告

一、 目的意义：把矩阵A 分解成一个下三角阵与一个上三角阵的乘积，即 A=LR ，其中L 为下三角阵，R 为上三角阵，这样原线性方程组就可以化为⎩⎨⎧==y Rx b Ly 的求解问题，方便求解。

二、 算法：1) 输入系数矩阵 A ； 2) 利用公式iji k ij ij ij r l a r ∑-=-=11和ki i k jk ji ji rr l a l /)(11∑-=-=交错进行，计算得出矩阵L 和R ；3) 回带到⎩⎨⎧==y Rx bLy 中得出原线性方程组的解。

三、 源程序：#include <stdio.h>#include <string.h> #include <math.h> #include <stdlib.h> #define N 100 main() {int i,j,k,s,n;printf("请输入系数矩阵A 的阶数：n= "); scanf("%d",&n); floata[N][N]={0},L[N][N]={0},R[N][N]={0},sigma1,sigma2,b[N],y[N],x[N];/*为L 主对角线元素赋1*/ for(i=0;i<n;i++) { L[i][i]=1; }printf("请输入系数矩阵A ：\n"); /*输入系数矩阵A*/报告1Doolittle 分解for(i=0;i<n;i++){for(j=0;j<n;j++)scanf("%f",&a[i][j]);}for(k=0;k<n;k++){for(j=k;j<n;j++) /*计算矩阵R*/{sigma1=0;for(s=0;s<=k-1;s++)sigma1+=L[k][s]*R[s][j];R[k][j]=a[k][j]-sigma1;}for(i=k;i<n;i++) /*计算矩阵L*/{sigma2=0;for(s=0;s<=k-1;s++)sigma2+=L[i][s]*R[s][k];L[i][k]=(a[i][k]-sigma2)/R[k][k];}}printf("\n A矩阵为：\n");/*输出矩阵L、R*/ for(i=0;i<n;i++){for(j=0;j<n;j++)printf("%5.1f ",a[i][j]);printf("\n");}printf("\n L矩阵为：\n");for(i=0;i<n;i++){for(j=0;j<n;j++)printf("%5.1f ",L[i][j]);printf("\n");}printf("\n R矩阵为：\n");for(i=0;i<n;i++){for(j=0;j<n;j++)printf("%5.1f ",R[i][j]);printf("\n");}printf("请输入b矩阵\n",i+1);for(i=0;i<n;i++)scanf("%f",&b[i]);for(i=0;i<n;i++)/*回代法求解方程组Ly=b*/{sigma1=0;for(k=0;k<=i-1;k++)sigma1+=L[i][k]*y[k];y[i]=b[i]-sigma1;}for(i=n-1;i>=0;i--){sigma2=0;for(k=i+1;k<n;k++)sigma2+=R[i][k]*x[k];x[i]=(y[i]-sigma2)/R[i][i];}printf("解得x为：\n");for(i=0;i<n;i++)printf("%5.1f ",x[i]);printf("\n");}四、计算结果与分析：分析：运行结果与预想的结果相近，误差对结果的影响不是很大，比较理想五、 参考文献：[1]刑志栋. 矩阵数值分析. 陕西： 陕西科学技术出版社， 2005 [2]谭浩强. C 语言程序设计. 北京：清华大学出版社，2005一、 目的意义：对称正定矩阵是实践中经常遇到的一种特殊矩阵类型矩阵，由于矩阵本身的流量好兴致，使得cholesky 分解在存储和运算量上较一般消去法节省一半左右，且解的精度高，cholesky 分解方法是目前计算机求解该类问题最有效的方法之一。

二、 算法：1) 输入系数矩阵 A ；2) 利用公式 iki k ik ii ii jjiij ij jkj k ik ij ij l s a d d s l l s a s ∑∑-=-=-==-=1111/交错进行，计算得出矩阵L 和D ；3) 回带到⎩⎨⎧==yD x L bLy TT 中得出原线性方程组的解 三、 源程序：#include <stdio.h>#include <stdlib.h> #include <math.h> #define EPS 1.0e-8 #define N 20double a[N][N], b[N], x[N];报告2cholesky 分解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; i<n; i++){for (int j=0; j<n; j++)scanf("%lf", &a[i][j]);getchar();}printf("\n输入线性方程组右端项b[%d]:\n", n);for (i=0; i<n; i++){scanf("%lf", &b[i]);}getchar();for (i=0; i<n-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; i<n; i++){printf("x%d=%lf ", i+1, x[i]);printf("\n");printf("\n");system("pause");}int zhuyuan(int row){double elem = a[row][row];int rowmark = row;for (int i=row+1; i<n; i++){if (elem<a[i][row]){elem = a[i][row];rowmark = i;}}if(fabs(elem) <= EPS){return -1;}return rowmark;}void hangjiaohuan(int row1, int row2) {double tmp;tmp = b[row1];b[row1] = b[row2];b[row2] = tmp;for (int j=0; j<n; j++){tmp = a[row1][j];a[row1][j] = a[row2][j];a[row2][j] = tmp;}}void xiaoyuan(int row){for (int i=row+1; i<n; i++){int j=row;double tmp = a[i][j]/a[row][row];b[i] -= tmp*b[row];for (a[i][j++] = 0; j<n; j++){a[i][j] -= tmp*a[row][j];}}}void huidai(){x[n-1] = b[n-1]/a[n-1][n-1];for (int i=n-2; i>=0; i--){double sum = 0.0;for (int j=i+1; j<n; j++){sum -= a[i][j]*x[j];}x[i] = (b[i]+sum)/a[i][i];}}四、计算结果与分析：分析：运行结果与预想的结果相近，误差对结果的影响不是很大，比较理想五、参考文献：[1]刑志栋. 矩阵数值分析. 陕西：陕西科学技术出版社，2005[2]谭浩强. C 语言程序设计. 北京：清华大学出版社，2005一、 目的意义：通过有限次的算术运算得到线性方程组的近似值二、 算法：1) 输入系数矩阵 A ；2) 输入迭代的初始向量；3) 利用公式()()ii ni j j k j ij i k i a x a b x 1*11⎪⎪⎪⎭⎫ ⎝⎛-=∑≠=+，计算得出原线性方程组的近似解。

三、 源程序：#include<stdio.h>#include<conio.h> #include<malloc.h> #include<math.h> #define EPS 1e-6 #define MAX 100 #define N 100float *Jacobi(float a[N][N],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;i<n;i++) scanf("%f",&x[i]); while(1)报告3Jacobi 迭代法{epsilon=0;k++;for(i=0;i<n;i++){s=0;for(j=0;j<n;j++){if(j==i) continue;s+=a[i][j]*x[j];}y[i]=(a[i][n]-s)/a[i][i];epsilon += fabs(y[i]-x[i]);}if(epsilon < EPS){printf("\n迭代次数为：%d次\n",k);return y;}if(k>=MAX){printf("The Method is disconvergent");return y;}for(i=0;i<n;i++)x[i]=y[i];}}main(){int i,j,n;float a[N][N];float *x;printf("请输入系数矩阵的阶数：n=");scanf("%d",&n);printf("请按行输入线性方程组的增广矩阵：\n");for(i=0;i<n;i++){for(j=0;j<n+1;j++)scanf("%f",&a[i][j]);}x=(float *)malloc(3*sizeof(float));x=Jacobi(a,n);printf("\n原线性方程组的解为：\n");for(i=0;i<n;i++)printf("x[%d]=%f\n",i,x[i]);getch();}四、计算结果与分析：分析：运行结果与预想的结果相近，误差对结果的影响不是很大，比较理想五、参考文献：[1]刑志栋. 矩阵数值分析. 陕西：陕西科学技术出版社，2005[2]谭浩强. C语言程序设计. 北京：清华大学出版社，2005报告4乘幂法一、目的意义：阶数超过四次的多项式的根一般不能用有限次运算求得，因而矩阵特征值的计算方法本质上总是采用迭代格式，而采用乘幂法可以计算最大的一个或按模最大的几个特征值和特征相应特征向量。