杜立特分解法

2013-2014（1）專業(yè)課程實(shí)踐論文題目：杜立特分解法、算法理論n階線性方程組的系數(shù)矩陣A非奇異且有分解式A=LR，其中L為單位下三角矩陣，R為上三角矩陣，即L=(1)，當(dāng)i<j時(shí)，l=0；R=(r)，當(dāng)i>jijijij時(shí)，r=0，矩陣A的這種分解方法為Doolittle的分解。ijaa…a■11■rr…rii121n11111121naa…a11r…r21???22???2n??????=21??????1222n??????aa…a11…11[rn1n2nn」n1n2—1nn比較等號(hào)兩邊的第i行和第j列的元素，得a丄億。k=1因?yàn)?=???=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從而當(dāng)rij二a一工1r，j=z,z+l,???n，j=z+1,z+2,?…n日寸,zjikkjk=1于是就得到了計(jì)算LR分解的一般計(jì)算公式。二、算法框圖三、算法程序#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;}四、算法實(shí)現(xiàn)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