c语言 三角分解解矩阵方程,计算方法：矩阵三角分解法解线性方程组

#include#define N 100

int print(double a[N][N],int n)

{

int i,j;

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

{

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

printf("%.4f\t",a[i][j]);

printf("\n");

}

return 0;

}

int main()

{

double a[N][N+1]={0},l[N][N]={0},u[N][N]={0},y[N]={0},x[N]={0},sum;

int i,j,k,n;

printf("请输入方程的阶数:\n");

scanf("%d",&n);

printf("请输入方程增广矩阵:\n");

for(i=1;i<=n;i++)//输入方程增广矩阵

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

scanf("%lf",&a[i][j]);

for(i=1;i<=n;i++)//初始化l数组的对角线上的元素，全部为1

l[i][i]=1;

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

{

for(i=k;i<=n;i++)//用L的第K行去乘U的每一列，得到U的第K行，i是第K行的第i列

{

for(j=1,sum=0;j<=k-1;j++)//j是用来计算第i列之前的数的和

sum=sum+l[k][j]*u[j][i];

u[k][i]=a[k][i]-sum;

}

for(i=k+1;i<=n;i++)//用L的每一行去乘U的第K列，得到L的第K列，其中i是第K列的第i行

{

for(j=1,sum=0;j<=k-1;j++)//j是指第i列前面的数的和

sum=sum+l[i][j]*u[j][k];

l[i][k]=(a[i][k]-sum)/u[k][k];

}

}

printf("L阵为:\n");

print(l,n);

printf("U阵为:\n");

print(u,n);

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

{

for(j=1,sum=0;j<=i-1;j++)

sum=l[i][j]*y[j]+sum;//注意乘的是y[j]而不是y[i]，

y[i]=a[i][n+1]-sum;

}printf("中间向量y为:\n");

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

printf("y[%d]=%f\n",i,y[i]);

for(i=n;i>=1;i--)

{

for(j=n,sum=0;j>i;j--)

sum=u[i][j]*x[j]+sum;

x[i]=(y[i]-sum)/u[i][i];

}

printf("方程的根为:\n");

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

printf("x%d=%f\n",i,x[i]);

system("pause");

return 0;

}

/*3

11 -3 -2 3

-23 11 1 0

1 -2 2 -1

x3=1.001,x2=2.001,x1=1.000

3

1 4 -5 3

1 3 -2 2

6 -1 18 2

x1=1.333333.x2=0,x3=-0.33333

3

1 -1 1 -4

5 -4 3 -12

2 1 1 11

x1=3,x2=6,x3=-1

3

2 3 5 5

3 4 7 6

1 3 3 5

x1=-4,x2=1,x3=2(书上答案有误)

4

1.1161 0.1254 0.1397 0.1490 1.5471

0.1582 1.1675 0.1768 0.1871 1.6471

0.1968 0.2071 1.2168 0.2271 1.7471

0.2368 0.2471 0.2568 1.2671 1.8471

x1=1.0406,x2=0.9870,x3=0.9351,x4=0.8813

3

2 1 2 6

4 5 4 18

6 -3 5 5

x1=1,x2=2,x3=1

*/