#include <stdio.h>
#include <stdlib.h>
#include<stdbool.h>
#define MAX 10
double A[MAX][MAX];
double b[MAX];
double X[MAX];
double Y[MAX];
double U[MAX][MAX];
double L[MAX][MAX];
int NUM;
void Input_Matrix()//输入矩阵
{
int i,j;
printf("系数矩阵A的阶数:\n");
scanf("%d",&NUM);
for(i=1; i<=NUM; i++)
{
printf("系数矩阵A的第%d行元素:\n",i);
for(j=1; j<=NUM; j++)
scanf("%lf",&A[i-1][j-1]);
}
printf("右端项b:\n");
for(i=1; i<=NUM; i++)
{
scanf("%lf",&b[i-1]);
}
printf("\n");
printf("输入的系数矩阵A:\n");
for(i=0; i<NUM; i++)
{
for(j=0; j<NUM; j++)
printf("%.4lf\t",A[i][j]);
printf("\n");
}
printf("\n");
printf("输入的右端项b:\n");
for(i=0; i<NUM; i++)
printf("%.4lf\n",b[i]);
}
bool Direct()
{
bool flag = true;
int r,i,k;
double sum_u, sum_l;
for(r=0; r<NUM; r++)
U[0][r]=A[0][r];
for(r=1; r<NUM; r++)
L[r][0] = A[r][0]/A[0][0];
for(r=1; r<NUM; r++)
{
for(i=r; i<NUM; i++)
{
sum_u = 0;
for(k=0; k<r; k++)
sum_u += L[r][k]*U[k][i];
U[r][i] = A[r][i]-sum_u;
}
for(i=r+1; i<NUM&&r!=NUM-1; i++)
{
sum_l = 0;
for(k=0; k<r; k++)
sum_l += L[i][k]*U[k][r];
L[i][r] = (A[i][r]-sum_l)/U[r][r];
}
}
for(i=0; i<NUM; i++)
L[i][i]=1.0;
for(r=0; r<NUM; r++)
{
if(U[r][r]==0) //等价于判断矩阵的顺序主子式是否为0
{
flag = false;
return flag;
}
}
//回代求解
Y[0] = b[0];
double sum;
for(i=1; i<NUM; i++)
{
sum = 0;
for(k=0; k<i; k++)
sum += L[i][k]*Y[k];
Y[i] = b[i]-sum;
}
X[NUM-1] = Y[NUM-1]/U[NUM-1][NUM-1];
for(i=NUM-2; i>=0; i--)
{
sum = 0;
for(k=i+1; k<NUM; k++)
sum += U[i][k]*X[k];
X[i] = (Y[i]-sum)/U[i][i];
}
return flag;
}
int main()
{
int i,r;
Input_Matrix();
if (Direct())
{
printf("\n");
printf("Result:\n");
for (i=0; i<NUM; i++)
printf("%.4lf\n",X[i]);
printf("\n");
printf("L:\n");
for(i=0; i<NUM; i++)
{
for(r=0; r<NUM; r++)
printf("%.4lf\t",L[i][r]);
printf("\n");
}
printf("\n");
printf("U:\n");
for(i=0; i<NUM; i++)
{
for(r=0; r<NUM; r++)
printf("%.4lf\t",U[i][r]);
printf("\n");
}
}
else
printf("Failed.n");
return 0;
}
C语言实现LU三角分解法(计算方法)
最新推荐文章于 2021-05-23 23:05:10 发布