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

/*三对角阵的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
=2;a[n+1]=1;c
=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;

}