三角分解法

#include<stdio.h>
#include<math.h>
#define N 4 /*定义方程组系数矩阵的阶数N*/
static float a[N][N]={{1.5,0.5,2.5,-3},{6,0.6,0,-5},{1,1.2,14,1},{-3,-5,1,11.1}};/*输入被求解线性方程组*/
static float b[N]={14.4095,16.166,24.628,-5.8256};
float x[N]={0,0,0,0}; /*定义初始值*/
float U[N][N]={{1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1}};
void Lu(float a[N][N],float L[N][N],float U[N][N]) /* LU分解*/
{
         int i,j,k;
         for(k=0;k<N;k++)
         {
                  for(i=k;i<N;i++) /* 计算L矩阵的第k列元素*/
                  {
                          L[i][k]=a[i][k];
                          for(j=0;j<=k-1;j++)
                                   L[i][k]-=(L[i][j]*U[j][k]);
                  }
                  for(j=k+1;j<N;j++) /* 计算U矩阵的第k列元素*/
                  {
                          U[k][j]=a[k][j];
                          for(i=0;i<=k-1;i++)
                                   U[k][j]-=(L[k][i]*U[i][j]);
                          U[k][j]/=L[k][k];
                  }
         }
}
void solver(float b[N],float L[N][N],float U[N][N],float x[N])/*三角分解法计算核心*/
{
         int j,k;
         float y[N];
         for(k=0;k<N;k++) /* 计算Ly=b中的y */
         {
                  y[k]=b[k];
                  for(j=0;j<=k-1;j++)
                          y[k]=y[k]-L[k][j]*y[j];
                  y[k]=y[k]/L[k][k];
         }
         for(k=N-1;k>=0;k--) /* 计算Ux=y中的x */
         {
                  x[k]=y[k];
                  for(j=k+1;j<N;j++)
                          x[k]=x[k]-U[k][j]*x[j];
         }
}
void augm_matrix(float a[N][N],float b[N]) /* 输出增广矩阵*/
{
         int i,j;
         for(i=0;i<N;i++)
         {
                  for(j=0;j<N;j++)
                          printf("%10f",a[i][j]);
                  printf("%10f",b[i]);
                  printf("\n");
         }
         printf("\n");
}
main()
{
         float L[N][N];
         int i;
         printf("Linear equations:\n");
         augm_matrix(a,b);
         Lu(a,L,U);
         solver(b,L,U,x);
         printf("Results:\n");
         for(i=0;i<N;i++) /* 输出方程组的解*/
                  printf(" x%d=%f\n",i+1,x[i]);
         system("pause");
}

参考资料:
[1] 王汉青. 暖通空调流体流动数值计算方法与应用[M]. 北京: 科学出版社, 2013.10.

来源:,欢迎分享本文,转载请保留出处!