杜立特尔分解法 MATLAB,杜立特分解法.doc

杜立特分解法杜立特分解法

2013-2014(1)专业课程实践论文

题目：杜立特分解法

一、算法理论

阶线性方程组的系数矩阵非奇异且有分解式，其中为单位下三角矩阵，为上三角矩阵，即，当时，；，当时，，矩阵的这种分解方法为Doolittle的分解。

比较等号两边的第行和第列的元素，得。

因为，所以，，

从而当，，时，，从而，

于是就得到了计算LR分解的一般计算公式。

二、算法框图

三、算法程序

#include

#include

#define N 3

using namespace std;

int main()

{

double A[N+1][N+1]={{0,0,0,0},{0,2,1,1},{0,1,3,2},{0,1,2,2}};

double L[N+1][N+1]={0};

double U[N+1][N+1]={0};

double b[N+1]={0,4,6,5};

double y[N+1];

double x[N+1];

int i,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++)

{double l=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++)

{double l=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++)

{ double l=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--)

{

double l=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<

for(i=1;i<=N;i++)

{for(j=1;j<=N;j++)

cout<

cout<

cout<

cout<

for(i=1;i<=N;i++)

{for(j=1;j<=N;j++)

cout<

cout<

cout<

cout<

for(i=1;i<=N;i++)

cout<

return 1;

}

四、算法实现

例1．用杜立特分解法，求解三元方程组。

解：

例2．用杜立特分解法，求解三元方程组

解：

开始

输入矩阵行数

输入矩阵A第各行元素

分解矩阵A得L

输入矩阵B第各行元素

分解矩阵A得R

求得Y

求得X