crout分解法(1).doc

Crout 分解法解线性方程组的算法及程序设计【摘要】在自然科学和工程技术中很多问题的解决常常归结为解线性方程组,而方程组的系数矩阵大致分为两种,一种是低阶稠密矩阵(阶数不超过150),另一种是大型稀疏矩阵(矩阵阶数高且零元素较多).解低阶稠密矩阵和某些特殊形式的大型稀疏矩阵宜应用直接法,例如:高斯消元法、矩阵三角分解法。对于一般形式的系数矩阵，Crout分解法是一种有效的方法，本文就Grout分解法给出其计算公式的推导过程及算法程序。【正文】一、Crout 方法解线性方程组的算法给定线性方程组Ax=b ,其中系数矩阵A=(aij)nn 为可逆的,x=(x1 ,x2 ,xn)T ,b =( b1,b2,bn)T 为常数项列向量。Crout分解后A=LU,L和U的结构为(1) 公式推导： 由A=LU及矩阵乘积和相等概念，有ai1=(li1,li2,lii,0,0)=li1 (i=1,2,n)a1j=(l11,0,0,00)=l11u1j (j=1,2,n)得：li1=ai1;u1j=a1j/l11 (i,j=1,2,n) 当ij时aij=(l11,li2,lii,00)=当ij时aij=(l11,li2,lii,00)= 得：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 = ( 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,nb)再解单位上三角型方程组Ux=Y 。由UX=Y得x1+u12x2+u1nxn=y1，x2+u2nxn=y2xn-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=20int main()/主函数int i,j,k;int size;float aNN,lNN,uNN;float bN,xN,yN;/定义变量printf(tttCrout分解法解方程组n);printf(请输入方阵A的n：);scanf(%d,&size);printf(n);printf(请输入方程组的系数:n);for(i=0;isize;i+)for(j=0;jsize;j+)scanf(%f,&aij);/输入方程组系数矩阵aprintf(n请输入方程组的y:n);for(i=0;isize;i+)scanf(%f,&bi);/输入结果矩阵bprintf(n方阵A为：n);for(i=0;isize;i+)for(j=0;jsize;j+)printf(%f ,aij);/输出aprintf(n);printf(n方程组y为:n);for(i=0;isize;i+)printf(%f ,bi);/输出bprintf(n);for(i=0;isize;i+)uii=1;/定初始值 令uii=1for(i=0;isize;i+)for(j=i+1;jsize;j+)lij=0;/定初始值 令lij=0for(j=0;jsize;j+)for(i=j+1;isize;i+)uij=0;/定初始值 令uij=0l00=a00;for(i=1;isize;i+)li0=ai0;/计算第一行的lu0i=a0i/l00;/计算第一列的ufor(i=1;isize-1;i+)for(j=1;j=i;j+)/计算第2行到第size-1行的llij=aij;for(k=0;kj;k+)lij=lij-lik*ukj; printf(n);for(j=i+1;jsize;j+)/计算第2行到第size行的uuij=aij;for(k=0;k=i-1;k+)uij=uij-lik*ukj;uij=uij/lii;printf(n);for(j=1;jsize;j+)/计算第size行的llsize-1j=asize-1j;for(k=0;k=j-1;k+)lsize-1j=lsize-1j-lsize-1k*ukj;printf(n); printf(输出矩阵Lijn); for(i=0;isize;i+) for(j=0;jsize;j+) printf(%f,lij); printf( );/输出下三角矩阵l printf(n); printf(输出矩阵Uijn);for(i=0;isize;i+) for(j=0;jsize;j+) printf(%f,uij); printf( );/输出单位上三角矩阵u printf(n);y0=b0/l00;/给y0初始值for(i=1;isize;i+)/计算yi的值yi=bi;for(k=0;k=i-1;k+)yi=yi-lik*yk;/计算公式 yi=yi/lii;printf(n);p