crout分解法
Crout 分解法解线性方程组的算法及程序设计
【摘要】在自然科学和工程技术中很多问题的解决常常归结为解线性方程组,而方程组的系数矩阵大致分为两种,一种是低阶稠密矩阵(阶数不超过150),另一种是大型稀疏矩阵(矩阵阶数高且零元素较多).解低阶稠密矩阵和某些特殊形式的大型稀疏矩阵宜应用直接法,例如:高斯消元法、矩阵三角分解法。对于一般形式的系数矩阵,Crout分解法是一种有效的方法,本文就Grout分解法给出其计算公式的推导过程及算法程序。
【正文】
一、Crout 方法解线性方程组的算法
给定线性方程组A=b ,其中系数矩阵A=(aij)n×n ,x=(x1 ,x2 ,…,xn)T ,b =( b1,b2,…bn)T 为常数项列向量。Crout分解后A=LU,L和U的结构为
公式推导:
由A=LU及矩阵乘积和相等概念,有
ai1=(li1,li2···,lii,0,···0)·=li1 (i=1,2,···n)
a1j=(l11,0,0,0···0)·=l11·u1j (j=1,2,···,n)
得:li1=ai1;u1j=a1j/l11 (i,j=1,2,···,n)
当i
aij=(l11,li2,···lii,0···0)·=
当ij时
aij=(l11,li2,···lii,0···0)·=
得:lij=aij- (j=2,3,···,i, i=2,3,···,n)
uij=(aij-)/lii (i=2,3,···,n,j=i+1,i+2,···,n)
(2)具体算法如下: 对A 作LU 分解由A = LU及矩阵的乘法原理可得: ij = aij-( j = 1, 2 , …, i, i=1,2,…n)
uij = ( aij - ) / lii (j = i + 1, i + 2 , …, n,i=1,2,…n)
②解两个三角型方程组
由A = LU 及Ax = b 可得L (Ux) = b ,令 Y = ( y 1 , y 2 , …, y n) T = Ux,则L Y = b ,于是求解Ax = b 就被化为求解下三角型方程组L Y = b 及单位上三角型方程组Ux= Y 。
a) 先解下三角型方程LY=b。由 l11y1=b1, l21y1+l22y2=b2l11,……ln1y1+ln2y2+…+lnnyn=bn
所以y1=b1/l11, yi=( bi-)/lii, i=2,3,…,n
b)再解单位上三角型方程组Ux=Y 。由UX=Y得x1+u12x2+…+u1nxn=y1,x2+…+u2nxn=y2
………xn-1+un-1nxn=yn-1,xn=yn,利用回代解法可得方程组AX=b的解为
xn=yn, xi=yi-,i=n-1,…,2,1
二、Crout 方法解线性方程组的程序
(1)程序代码:
#include "stdio.h"
#include "math.h"//头文件
#define N 20//自定义N=20
int main()//主函数
{
int i,j,k;
int size;
float a[N][N],l[N][N],u[N][N];
float b[N],x[N],y[N];//定义变量
printf("\t\t\tCrout分解法解方程组\n");
printf("请输入方阵A的n:");
scanf("%d",&size);
printf("\n");
printf("请输入方程组的系数:\n");
for(i=0;i
{
for(j=0;j
{
scanf("%f",&a[i][j]);//输入方程组系数矩阵a[][]
}
}
printf("\n请输入方程组的y:\n");
for(i=0;i
scanf("%f",&b[i]);//输入结果矩阵b[]
printf("\n方阵A[][]为:\n");
for(i=0;i
{
for(j=0;j
{
printf("%f ",a[i][j]);//输出a[][]
}
printf("\n");
}
printf("\n方程组y为:\n");
for(i=0;i
printf("%f ",b[i]);//输出b[]
printf("\n");
for(i=0;i
{
u[i][i]=1;//定初始值 令u[i][i]=1
}
for(i=0;i
fo