crout分解法c语言编程,crout分解法.doc

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