crout分解法1.docx
《crout分解法1.docx》由会员分享，可在线阅读，更多相关《crout分解法1.docx（11页珍藏版）》请在冰点文库上搜索。
crout分解法1
crout分解法
（1）
Crout分解法解线性方程组的算法及程序设计
【摘要】在自然科学和工程技术中很多问题的解决常常归结为解线性方程组,而方程组的系数矩阵大致分为两种,一种是低阶稠密矩阵（阶数不超过150）,另一种是大型稀疏矩阵（矩阵阶数高且零元素较多）.解低阶稠密矩阵和某些特殊形式的大型稀疏矩阵宜应用直接法,例如:
高斯消元法、矩阵三角分解法。
对于一般形式的系数矩阵，Crout分解法是一种有效的方法，本文就Grout分解法给出其计算公式的推导过程及算法程序。
【正文】
一、Crout方法解线性方程组的算法
给定线性方程组Ax=b,其中系数矩阵A=（aij）n×n为可逆的,x=（x1,x2,…,xn）T,b=（b1,b2,…bn）T为常数项列向量。
Crout分解后A=LU,L和U的结构为
aij=（l11,li2,···lii,0···0）·
=
当i
j时
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及矩阵的乘法原理可得:
lij=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=（y1,y2,…,yn）T=Ux，则LY=b,于是求解Ax=b就被化为求解下三角型方程组LY=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"//头文件
#defineN20//自定义N=20
intmain（）//主函数
{
inti,j,k;
intsize;
floata[N][N],l[N][N],u[N][N];
floatb[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;iscanf（"%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;iprintf（"%f",b[i]）;//输出b[]
printf（"\n"）;
for（i=0;i{
u[i][i]=1;//定初始值令u[i][i]=1
}
for（i=0;ifor（j=i+1;j{
l[i][j]=0;//定初始值令l[i][j]=0
}
for（j=0;jfor（i=j+1;i{
u[i][j]=0;//定初始值令u[i][j]=0
}
l[0][0]=a[0][0];
for（i=1;i{
l[i][0]=a[i][0];//计算第一行的l[][]
u[0][i]=a[0][i]/l[0][0];//计算第一列的u[][]
}
for（i=1;i{
for（j=1;j<=i;j++）//计算第2行到第size-1行的l[][]
{
l[i][j]=a[i][j];
for（k=0;k{
l[i][j]=l[i][j]-l[i][k]*u[k][j];
}
}
printf（"\n"）;
for（j=i+1;j{
u[i][j]=a[i][j];
for（k=0;k<=i-1;k++）
{
u[i][j]=u[i][j]-l[i][k]*u[k][j];
}
u[i][j]=u[i][j]/l[i][i];
}
printf（"\n"）;
}
for（j=1;j{
l[size-1][j]=a[size-1][j];
for（k=0;k<=j-1;k++）
{
l[size-1][j]=l[size-1][j]-l[size-1][k]*u[k][j];
}
}
printf（"\n"）;
printf（"输出矩阵L[i][j]\n"）;
for（i=0;ifor（j=0;j{
printf（"%f",l[i][j]）;printf（""）;//输出下三角矩阵l[][]
}
printf（"\n"）;
}
printf（"输出矩阵U[i][j]\n"）;
for（i=0;ifor（j=0;j{
printf（"%f",u[i][j]）;printf（""）;//输出单位上三角矩阵u[][]
}
printf（"\n"）;
}
y[0]=b[0]/l[0][0];//给y[0]初始值
for（i=1;i{
y[i]=b[i];
for（k=0;k<=i-1;k++）
{
y[i]=y[i]-l[i][k]*y[k];//计算公式
}
y[i]=y[i]/l[i][i];
}
printf（"\n"）;
printf（"y值:
\n"）;
for（i=0;iprintf（"y[%d]=%f",i+1,y[i]）;//输出y[i]的结果
printf（"\n\n"）;
printf（"x的值:
\n"）;
x[size-1]=y[size-1];//给x[size-1]赋值
for（i=size-2;i>=0;i--）
{
x[i]=y[i];
for（k=i+1;k{
x[i]=x[i]-u[i][k]*x[k];//计算x[]
}
}
for（i=0;i{
printf（"x[%d]=%f\n",i+1,x[i]）;//输出x[i]的结果
}
}
（2）程序应用
升级会员 