三对角阵的LU分解和三对角方程组的求解(C语言)

原创 2012年03月30日 18:42:32

/*三对角阵的LU分解和三对角方程组的求解


-------------A=LU的分解算法-------
参考教材:《数值分析》李乃成,梅立泉,科学出版社
    《计算方法教程》第二版 凌永祥,陈明逵
*/
#include<stdio.h>
#include<stdlib.h>
#include<math.h>

int main(void)
{
 int i,j,n;
 int N;

 printf("请输入 N(10,20,30或任意值): ");
 scanf("%d",&N);
 float *a=(float *)malloc(sizeof(float)*N);
 float *b=(float *)malloc(sizeof(float)*N);
 float *c=(float *)malloc(sizeof(float)*(N-1));
 float *d=(float *)malloc(sizeof(float)*N);
 float *x=(float *)malloc(sizeof(float)*N);
 float *y=(float *)malloc(sizeof(float)*N);
 float *u=(float *)malloc(sizeof(float)*N);
 float *l=(float *)malloc(sizeof(float)*N);
 a[0]=0;
 for(n=0;n<N-1;n++)
 {
  b[n]=2;a[n+1]=1;c[n]=1;
 }
 b[N-1]=2;
 d[0]=1;d[N-1]=-1;
 for(i=1;i<N-1;i++)
 {
  d[i]=0;
  d[N-1]=d[N-1]*(-1);
 }
 //------追过程-------
 u[0]=b[0];y[0]=d[0];
 l[0]=0;
 for(i=1;i<N;i++)
 {
  l[i]=a[i]/u[i-1];
  u[i]=b[i]-l[i]*c[i-1];
  y[i]=d[i]-l[i]*y[i-1];
 }
 printf("三对角阵A进行LU分解后的结果:\n\n");
 //U的上次主对角线
 printf("U的上次主对角线元u[]:\n");
 for(i=0;i<N;i++)
  printf("%4.3f  ",u[i]);
 printf("\n\n");

 //L的下次主对角线
 printf("L的下次主对角线元l[]:\n");
 for(i=0;i<N;i++)
  printf("%4.3f  ",l[i]);
 printf("\n\n");
 //打印上三角阵U
 printf("上三角阵U:\n");
 for(i=0;i<N;i++)
 {
  for(j=0;j<N;j++)
  {
   if(j==i)
   {
    printf("%4.3f  ",u[i]);
   }
   else if(j==(i+1))
   {
    printf("%4.3f  ",c[i]);
   }
   else
   {
    printf("%4.3f  ",l[0]);
   }
  }
  printf("\n");
 }
 //打印单位下三角阵L
 printf("单位下三角阵L:\n");
 float dig=1.0;
 for(i=0;i<N;i++)
 {
  for(j=0;j<N;j++)
  {
   if(j==i)
   {
    printf("%4.3f  ",dig);
   }
   else if(i==(j+1))
   {
    printf("%4.3f  ",l[i]);
   }
   else
   {
    printf("%4.3f  ",l[0]);
   }
  }
  printf("\n");
 }
  printf("\n");
 //------赶过程-------
 x[N-1]=y[N-1]/u[N-1];
 for(i=N-2;i>=0;i--)
 {
  x[i]=(y[i]-c[i]*x[i+1])/u[i];
 }

 printf("中间解向量y[]:\n");
 for(i=0;i<N;i++)
  printf("%4.3f  ",y[i]);
 printf("\n\n");
 printf("解向量x[]:\n");
 for(i=0;i<N;i++)
  printf("%4.3f  ",x[i]);
 printf("\n\n");
 free(a);free(b);free(c);free(d);free(x);free(y);free(l);
 
 return 0;
}

相关文章推荐

LU分解法解线性方程组(C语言)

#include #include #define N 10 //矩阵大小范围 /** 使用已经求出的x,向前计算x(供getx()调用)* float a[][] 矩阵U* float x[] ...
  • linvo
  • linvo
  • 2009年02月15日 18:33
  • 5114

Delphi7高级应用开发随书源码

  • 2003年04月30日 00:00
  • 676KB
  • 下载

LU分解的实现

LU分解是将矩阵分解为一个下三角矩阵和一个上三角矩阵的乘积。矩阵可以不是NxN的矩阵 一个可逆矩阵可以进行LU分解当且仅当它的所有子式都非零。如果要求其中的L矩阵(或U矩阵)为单位三角矩阵,那么分解是...

带状对角矩阵的LU分解及回代求解算法实现

算法名称:带状对角矩阵的LU分解及回代求解   算法描述:        分解主要是使用笔者前面几篇文章提到过的Crout方法。因为不可能把一个带状对角矩阵A的LU分解也像其压缩形式本是一样紧凑...

追赶法解三对角方程组

  • 2010年07月23日 20:40
  • 2KB
  • 下载

Matlab实现——严格对角占优三对角方程组求解(高斯赛尔德Gauss-Seidel迭代、超松弛)

严格对角占优三对角方程组求解 对中等规模的n阶的(n)线性方程组,直接法的准确性和可靠性,所以常采用直接法 对于较高阶的方程组,特别是地于某些偏微分方程离散化后得到的大型稀疏方程组(系统矩 阵绝...

MATLAB追赶法求解三对角方程

  • 2017年03月28日 22:13
  • 350B
  • 下载

追赶法求解三对角矩阵

  • 2013年06月18日 21:12
  • 431B
  • 下载

三元组表示三对角矩阵

题目设计算法求三对角矩阵在压缩存储下的转置矩阵代码#include #include #define MAX 100typedef struct { int row, col; //非零...
内容举报
返回顶部
收藏助手
不良信息举报
您举报文章:三对角阵的LU分解和三对角方程组的求解(C语言)
举报原因:
原因补充:

(最多只允许输入30个字)