杜立特分解法

2013-2014（1）专业课程实践论文题目：杜立特分解法、算法理论n阶线性方程组的系数矩阵A非奇异且有分解式A=LR，其中L为单位下三角矩阵，R为上三角矩阵，即L=(1)，当i<j时，l=0；R=(r)，当i>jijijij时，r=0，矩阵A的这种分解方法为Doolittle的分解。ijaa…a■11■rr…rii121n11111121naa…a11r…r21•••22•••2n••••••=21••••••1222n••••••aa…a11…11[rn1n2nn」n1n2—1nn比较等号两边的第i行和第j列的元素，得a丄亿。k=1因为1=•••=1=r=•••=r=0，所以a=工1r=工1r+r，z<j,i,i+1inj+l,jzjzjikQikQzjk=1k=1aij=工1rjkkjkaij=工1rjkkjk=1工1r+1r,jkkjjiiik=1从而1ji=(aji工1r)/r,j=z+1,•…,n,jkkiiik=1从而当rij二a一工1r，j=z,z+l,・・・n，j=z+1,z+2,•…n日寸,zjikkjk=1于是就得到了计算LR分解的一般计算公式。二、算法框图三、算法程序#include<iostream>#include<stdlib.h>#defineN3usingnamespacestd;intmain(){doubleA[N+1][N+1]={{0,0,0,0},{0,2,1,1},{0,1,3,2},{0,1,2,2}};doubleL[N+1][N+1]={0};doubleU[N+1][N+1]={0};doubleb[N+1]={0,4,6,5};doubley[N+1];doublex[N+1];inti,j,k,p;for(j=1;j<=N;j++)U[1][j]=A[1][j];L[1][1]=1;for(i=2;i<=N;i++){L[i][1]=A[i][1]/U[1][1];L[i][i]=1;}for(k=2;k<=N;k++){for(j=k;j<=N;j++){doublel=0;for(p=1;p<=k-1;p++)l+=L[k][p]*U[p][j];U[k][j]=A[k][j]-l;}for(i=k+1;i<=N;i++){doublel=0;for(p=1;p<=k-1;p++)l+=L[i][p]*U[p][k];L[i][k]=(A[i][k]-l)/U[k][k];}}y[1]=b[1];for(k=2;k<=N;k++){doublel=0;for(j=1;j<=k-1;j++)l+=L[k][j]*y[j];y[k]=b[k]-l;}x[N]=y[N]/U[N][N];for(k=N-1;k>=1;k--){doublel=0;for(j=k+1;j<=N;j++)l+=U[k][j]*x[j];x[k]=(y[k]-l)/U[k][k];}printf(”向量y为:”)；for(i=1;i<=N;i++)printf("%.1lf\t",y[i]);printf("\n");cout<<"L："<<endl；for(i=1；i<=N；i++){for(j=1；j<=N；j++)cout<<L[i][j]<<"\t"；cout<<endl；}cout<<endl；cout<<"U："<<endl；for(i=1；i<=N；i++){for(j=1；j<=N；j++)cout<<U[i][j]<<"\t"；cout<<endl；}cout<<endl；coutvv"方程组的解："vvendl;for(i=1；i<=N；i++)coutvv"x"vvivv"="vvx[i]vv"\t";return1;}四、算法实现2x+x+x=4TOC\o"1-5"\h\z例1.用杜立特分解法，求解三元方程组L+3x2+2X=6。123x+2x+2x=5V123解：2x+x+5x=11123例2.用杜立特分解法，求解三元方程组（4x+1x+12x=27123—2x—4x+5x=12V123解:「D:\Documents\Unttledl.exeTOC\o"1-5"\h