不选主元的三角分解法.doc

2012-2013（1）专业课程实践论文 不选主元的三角分解法 葛禹泽，0818180109，R数学08-1班 徐金龙，0818180215，R数学08-2班 一、算法理论 设为非奇异矩阵，且有分解式，其中为单位下三角阵，为上三角阵， 即 =…………………………(1) 下面说明的元素可以由步直接计算定出，其中第步定出的第行和的第行元素。由式（1），有 于是得的第一行元素； 于是得得第一行元素。 总结上述讨论，得到第行，用直接三角分解法解（要求所有顺序主子式都不为零）的计算公式: 二、算法框图 三、算法程序 #includestdio.h #includestdlib.h #includemath.h #define N 20 double A[N][N], x[N], y[N], b[N], L[N][N], U[N][N], f1=0, f2=0; int n; int i,j,k; void main() { printf(请输入方程的维数n!\n n = ); scanf(%d, n); getchar(); printf(\n输入%d行%d列矩阵\n, n, n); for (int i=0; in; i++) { for (int j=0; jn; j++) scanf(%lf, A[i][j]); getchar(); } printf(\n输入线性方程组右端项b[%d]: , n); for (i=0; in; i++) { scanf(%lf, b[i]); } getchar(); for(i=0;in;i++) for(j=0;jn;j++) {U[i][j]=0;L[i][j]=0;} for(i=0;in;i++) {U[0][i]=A[0][i];L[i][0]=A[i][0]/U[0][0];L[i][i]=1;} for(i=1;in;i++) for(j=i;jn;j++) {for(k=0;k=i-1;k++){f1=f1+L[i][k]*U[k][j];f2+=L[j][k]*U[k][i];} U[i][j]=A[i][j]-f1;L[j][i]=(A[j][i]-f2)/U[i][i];f1=0;f2=0;} y[0]=b[0]; for(i=1;in;i++) {for(j=0;j=i-1;j++) f1+=L[i][j]*y[j];y[i]=b[i]-f1;f1=0;} x[2]=y[2]/U[2][2]; for(i=1;i=0;i--) {for(j=i+1;jn;j++) f2+=U[i][j]*x[j];x[i]=(y[i]-f2)/U[i][i];f2=0;} printf(输出L矩阵：\n); for(i=0;in;i++) {for(j=0;jn;j++) printf(%f ,L[i][j]); printf(\n);} printf(输出U矩阵：\n); for(i=0;in;i++) {for(j=0;jn;j++) printf(%f ,U[i][j]); printf(\n);} printf(输出求解结果：\n); for(i=0;in;i++) printf(%f ,x[i]); printf(\n); } 四、算法实现 例1. 用直接三角分解法解 解：运行程序 (1) 显示出 请输入方程的维数!：输入的值为3，回车。 (2) 显示出 输入3行3列矩阵：输入1 2 3回车2 5 2回车3 1 5回车。 (3) 显示出 输入线性方程组右端项（3）：输入14 18 20回车。 (4) 显示结果： = = = 例2. 用直接三角分解法解 解：运行程序 (1) 显示出 请输入方程的维数!：输入的值为3，回车。 (2) 显示出 输入3行3列矩阵：输入0.001 2.000 3.000回车-1.000 3.712 4.623回车-2.000 1.072 5.643回车。 (3) 显示出 输入线性方程组右端项（3）：输入1.000 2.000 3.000回车。 (4) 显示结果： = = =