Doolittle分解

Doolittle分解

班级：计算062 姓名：王保翔 学号：3060811028

目的意义：在Gauss消元法中,对于n阶方程组,应用消去发经过n-1步消元之后,

得到一个与Ax=b等价的代数线性方程组A(n1)xb(n1),而且A(n1)为一个上三角矩阵.所以我们想是否能把矩阵A分解成一个下三角阵与一个上三角阵的乘积 A=LR,其中L为下三角阵,R为上三角阵.就变成

Lyb,了两个三角形方程组 的求解问题.

Rxy算法：

Setp1：利用for循环求出r[k][j]=a[k][j]-l[kp]r[kp]，l[ik]=（a[ik]-l[ip]r[ip]）

p1p1k1k1/r[k][k]。

Step2：y[i]=b[i]-l[ik]y[k],得出x[i]=(y[i]-k1i1ki1r[ik]x[k])/r[i][i].

n源程序：

#include #define N 10 void main() {

int i,j,k,n;

double a[N+1][N+1],l[N+1][N+1],r[N+1][N+1]; double x[N+1],y[N+1],b[N+1]; double p,q,t,s; printf(\"Input n:\"); scanf(\"%d\ for(i=1;i<=n;i++) {

printf(\"Input b[%d]:\ scanf(\"%lf\ }

for(i=1;i<=n;i++) for(j=1;j<=n;j++) {

printf(\"Input a[%d][%d]:\ scanf(\"%lf\ }

for(i=1;i<=n;i++) {

1 / 3

真诚为您提供优质参考资料，若有不当之处，请指正。

for(j=i;j<=n;j++) {

for(k=1,p=0;k<=i-1;k++) p=p+l[i][k]*r[k][j]; r[i][j]=a[i][j]-p; }

for(j=i+1;j<=n;j++) {

for(k=1,q=0;k<=i-1;k++) q=q+l[j][k]*r[k][i]; l[j][i]=(a[i][j]-q)/r[i][i]; }

l[i][i]=1; }

for(i=1;i<=n;i++) {

for(j=i;j<=n;j++)

printf(\"r[%d][%d]:%lf\\n\ for(j=i+1;j<=n;j++)

printf(\"l[%d][%d]:%lf\\n\ }

for(i=1;i<=n;i++) {

for(k=1,t=0;k<=i-1;k++) t=t+l[i][k]*y[k]; y[i]=b[i]-t; }

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

printf(\"Output y[%d],b[%d]:%lf,%lf\\n\ for(i=n;i>=1;i--) {

for(k=i+1,s=0;k<=n;k++) s=s+r[i][k]*x[k]; x[i]=(y[i]-s)/r[i][i]; }

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

printf(\"Output x[%d]:%lf\\n\}

算例及结果分析：

2 / 3

真诚为您提供优质参考资料，若有不当之处，请指正。

分析：输入矩阵A对应输出单位下三角矩阵L和上三角矩阵R，将L乘以R得到的就是A，说明结果可靠。输入b[i]解方程组Ax=b。

参考文献

[1]刑志栋，矩阵数值分析， 陕西： 陕西科学技术出版社， 2005。 [2]谭浩强，C语言程序设计，北京：清华大学出版社，2005。

[3]翁惠玉，c语言程序设计思想与方法，北京：人民邮电出版社，2008

3 / 3