矩阵的三角分解的编程实现
一． 课程设计目的：
1． 将高斯消去法改写为紧凑形式，可以直接从矩阵A 的元素得到计算L,U 元素的递推公式，而不需要任何中间的步骤。
2． 理解三角分解的具体方法。
3． 掌握对于矩阵的运算的c 语言代码。了解二维数组的运用。
4． 掌握Doolittle 分解。
二． 基本知识回顾：
若n 阶方阵 n *n A ∈C n 的各阶顺序主子式不等于零，即：
a 11a 12 a 1k
a 21a 22 a 2k ∆k =≠0, (k =1, 2, , n ), （2）
ak 1ak 2 akk
则A 的LU 分解A =L ⨯U 存在且唯一。
⎡a 11 a 1r a 1n ⎤⎢ ⎥ ⎢⎥A =⎢a r 1 a rr a rn ⎥=⎢⎥ ⎢⎥⎢⎣a n 1 a nr a nn ⎥⎦
⎡1⎤⎡U 11 U 1r ⎢ ⎥⎢ ⎢⎥⎢⎢L r 1 1⎥⎢U rr ⎢⎥⎢ ⎢⎥⎢⎢⎣L n 1 L nr 1⎥⎦⎢⎣U 1n ⎤⎥⎥ U rn ⎥=LU ⎥ ⎥U nn ⎥⎦ （3）
由矩阵的乘法原理, 可推导出LU 分解的迭代算法
U 0j =a 0j ,(j =0,1,2, , n -1),
l i 0=（4） a i 0,(i =0,1,2, , n -1), （5） u 00
a rj -∑l ik u kj
k =1r -1u rj =l rr
r -1=a rj -∑l rk u kj , k =1r -1 （6） (r =0,1, 2, , n -1; j =r , , n -1),
a ir -∑l ik u kr
k =1l ir =u rr , （7）
(r =0,1, 2, , n -1; i =r +1, , n -1)
矩阵的LU 分解是一个循环迭代的过程, U矩阵是从第1行迭代到第n 行, 而L 矩阵则是从第1列迭代到第n 列, 且U 矩阵先于L 矩阵
一个节拍。
三． 基本问题：
编程实现Doolittle 分解。
四． 算法描述：
1. 输入矩阵A;
3. 若能分解，进行LU 分解。L 为单位下三角矩阵，U 为上三角矩阵。
五． 程序实现：
#include
#define N 10
int main()
{ float l[N][N]={0};
float u[N][N]={0};
float y[N]={0};
float x[N]={0};
float a[N][N];
float b[N];
float sum=0;
int i,j,k;
int n;
int flag=1;
while(flag)
{
printf("请输入系数矩阵的大小：");
scanf("%d", &n);
if(n>N){
printf("矩阵过大！\n");
continue;
}
flag=0;
}
printf("请输入系数矩阵值：\n");
for(i=0; i
{
for(j=0; j
{
printf("a[%d][%d]: ", i, j);
scanf("%f", &a[i][j]);
}
}
printf("\n原始矩阵：\n");
for(i=0; i
{
for(j=0; j
printf("%0.3f ",a[i][j]);
printf("\n");
}
printf("\n\n");
for(i=0; i
{
for(j=0; j
{
if(i==j) l[i][j] = 1;
}
}
for(i=0; i
{
u[0][i] = (float)(a[0][i]/l[0][0]);
}
for(i=0; i
{
for(j=i+1; j
{
for(k=0,sum=0; k
{
if(k != i) sum += l[j][k]*u[k][i]; }
l[j][i] = (float)((a[j][i]-sum)/u[i][i]); }
for(j=i+1; j
{
for(k=0,sum=0; k
{
if(k != i+1) sum += l[i+1][k]*u[k][j]; }
u[i+1][j] = (float)((a[i+1][j]-sum)); }
}
printf("矩阵L ：\n");
for(i=0; i
{
for(j=0; j
{
printf("%0.3f ", l[i][j]);
}
printf("\n");
}
printf("\n矩阵U ：\n"); for(i=0; i
{
for(j=0; j
printf("%0.3f ", u[i][j]); }
}
}
printf("\n");
矩阵的三角分解的编程实现
一． 课程设计目的：
1． 将高斯消去法改写为紧凑形式，可以直接从矩阵A 的元素得到计算L,U 元素的递推公式，而不需要任何中间的步骤。
2． 理解三角分解的具体方法。
3． 掌握对于矩阵的运算的c 语言代码。了解二维数组的运用。
4． 掌握Doolittle 分解。
二． 基本知识回顾：
若n 阶方阵 n *n A ∈C n 的各阶顺序主子式不等于零，即：
a 11a 12 a 1k
a 21a 22 a 2k ∆k =≠0, (k =1, 2, , n ), （2）
ak 1ak 2 akk
则A 的LU 分解A =L ⨯U 存在且唯一。
⎡a 11 a 1r a 1n ⎤⎢ ⎥ ⎢⎥A =⎢a r 1 a rr a rn ⎥=⎢⎥ ⎢⎥⎢⎣a n 1 a nr a nn ⎥⎦
⎡1⎤⎡U 11 U 1r ⎢ ⎥⎢ ⎢⎥⎢⎢L r 1 1⎥⎢U rr ⎢⎥⎢ ⎢⎥⎢⎢⎣L n 1 L nr 1⎥⎦⎢⎣U 1n ⎤⎥⎥ U rn ⎥=LU ⎥ ⎥U nn ⎥⎦ （3）
由矩阵的乘法原理, 可推导出LU 分解的迭代算法
U 0j =a 0j ,(j =0,1,2, , n -1),
l i 0=（4） a i 0,(i =0,1,2, , n -1), （5） u 00
a rj -∑l ik u kj
k =1r -1u rj =l rr
r -1=a rj -∑l rk u kj , k =1r -1 （6） (r =0,1, 2, , n -1; j =r , , n -1),
a ir -∑l ik u kr
k =1l ir =u rr , （7）
(r =0,1, 2, , n -1; i =r +1, , n -1)
矩阵的LU 分解是一个循环迭代的过程, U矩阵是从第1行迭代到第n 行, 而L 矩阵则是从第1列迭代到第n 列, 且U 矩阵先于L 矩阵
一个节拍。
三． 基本问题：
编程实现Doolittle 分解。
四． 算法描述：
1. 输入矩阵A;
3. 若能分解，进行LU 分解。L 为单位下三角矩阵，U 为上三角矩阵。
五． 程序实现：
#include
#define N 10
int main()
{ float l[N][N]={0};
float u[N][N]={0};
float y[N]={0};
float x[N]={0};
float a[N][N];
float b[N];
float sum=0;
int i,j,k;
int n;
int flag=1;
while(flag)
{
printf("请输入系数矩阵的大小：");
scanf("%d", &n);
if(n>N){
printf("矩阵过大！\n");
continue;
}
flag=0;
}
printf("请输入系数矩阵值：\n");
for(i=0; i
{
for(j=0; j
{
printf("a[%d][%d]: ", i, j);
scanf("%f", &a[i][j]);
}
}
printf("\n原始矩阵：\n");
for(i=0; i
{
for(j=0; j
printf("%0.3f ",a[i][j]);
printf("\n");
}
printf("\n\n");
for(i=0; i
{
for(j=0; j
{
if(i==j) l[i][j] = 1;
}
}
for(i=0; i
{
u[0][i] = (float)(a[0][i]/l[0][0]);
}
for(i=0; i
{
for(j=i+1; j
{
for(k=0,sum=0; k
{
if(k != i) sum += l[j][k]*u[k][i]; }
l[j][i] = (float)((a[j][i]-sum)/u[i][i]); }
for(j=i+1; j
{
for(k=0,sum=0; k
{
if(k != i+1) sum += l[i+1][k]*u[k][j]; }
u[i+1][j] = (float)((a[i+1][j]-sum)); }
}
printf("矩阵L ：\n");
for(i=0; i
{
for(j=0; j
{
printf("%0.3f ", l[i][j]);
}
printf("\n");
}
printf("\n矩阵U ：\n"); for(i=0; i
{
for(j=0; j
printf("%0.3f ", u[i][j]); }
}
}
printf("\n");