#include<stdio.h>
#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
*/