lu.cpp

解矩阵方程

#include <stdio.h> 
#include <stdlib.h> 
#define N 10 //矩阵大小范围 
void main() 
{ 
float l[N][N]={0}; //定义L矩阵 
float u[N][N]={0}; //定义U矩阵 
float y[N]={0}; //定义数组Y 
float x[N]={0}; //定义数组X 
float a[N][N]={0}; //定义系数矩阵 
float b[N]={0}; //定义右端项 
float sum=0; 
int i,j,k,s; 
int n=3; 
int flag=1; 

 
while(flag) //输入矩阵
{ 
printf("请输入系数矩阵的大小："); 
scanf("%d", &n); 
if(n>N){ 
printf("矩阵过大！\n"); 
continue; 
} 
flag=0; 
} 
printf("请输入系数矩阵值：\n"); 
for(i=0; i<n; i++) 
{ 
for(j=0; j<n; j++) 
{ 
printf("a[%d][%d]: ", i, j); 
scanf("%f", &a[i][j]); 
} 
} 
printf("请输入右端项数组：\n"); 
for(i=0; i<n; i++) 
{ 
printf("b[%d]: ", i); 
scanf("%f", &b[i]); 
} 


printf("\n原始矩阵：\n");      /*显示原始矩阵*/ 
for(i=0; i<n; i++) 
{ 
for(j=0; j<n; j++) 
printf("%0.3f ",a[i][j]); 
printf("\n"); 
} 
printf("\n\n"); 


for(i=0; i<n; i++)             /*初始化矩阵l*/ 
{ 
for(j=0; j<n; j++) 
{ 
if(i==j) l[i][j] = 1; 
} 
} 


for(i=0;i<n;i++)               /*开始LU分解*/
{  u[0][i]=(float)(a[0][i]);}  /*第一步：对矩阵U的首行进行计算*/

for(k=0;k<n-1;k++)             /*第二步：逐步进行LU分解*/
{  for(i=k+1;i<n;i++)          /*对L的第k列进行计算*/
{  for(s=0,sum=0;s<n;s++)
{  if(s!=k)
   sum+=l[i][s]*u[s][k];
}
   l[i][k]=(float)((a[i][k]-sum)/u[k][k]);
}

for(j=k+1;j<n;j++)             /*对U的第k+1行进行计算*/
{  for(s=0,sum=0;s<n;s++)
{  if(s!=k+1)
   sum+=l[k+1][s]*u[s][j];
}
u[k+1][j]=(float)((a[k+1][j]-sum));
}
}
 
printf("矩阵L：\n");           /*输出矩阵l*/
for(i=0; i<n; i++) 
{ 
for(j=0; j<n; j++) 
{ 
printf("%0.3f ", l[i][j]); 
} 
printf("\n"); 
} 


printf("\n矩阵U：\n");         /*输出矩阵u*/  
for(i=0; i<n; i++) 
{ 
for(j=0; j<n; j++) 
{ 
printf("%0.3f ", u[i][j]); 
} 
printf("\n"); 
} 

y[0]=b[0] ;                 /*回代法计算数组Y*/  
for(i=1;i<n;i++)
{ for(j=0,sum=0;j<i;j++)
{ sum+=y[j]*l[i][j];
}
  y[i]=(float)(b[i]-sum);
}

printf("\n\n数组Y：\n");      /*显示数组Y*/ 
for(i=0; i<n; i++) 
{ 
printf("y%d = %0.3f\n", i+1,y[i]); 
} 

x[n-1]=(float)(y[n-1]/u[n-1][n-1]);  /*回代法计算数组X*/ 
for(i=n-2;i>=0;i--)
{  for(j=n-1,sum=0;j>i;j--)
{  sum+=x[j]*u[i][j];}
x[i]=(float)((y[i]-sum)/u[i][i]);
}

 
printf("\n\n数组X：\n");     /*显示数组X*/
for(i=0; i<n; i++) 
{ 
printf("x%d = %0.3f\n", i+1,x[i]); 
} 
}