lu.h

数值计算利用fft分解法来求解线性方程组,简明方便

//LU分解法解线性方程组
#include"stdio.h"
void Desolute_LU(double *A,int N,double *B,double *Root)
{int i,j,k;double temp;
double *C=new double[N*(N+1)];
for(i=0;i<N;i++)
{
	for(j=0;j<N+1;j++)
	{
	if(j<N) *(C+i*(N+1)+j)=*(A+i*N+j);
	else *(C+i*(N+1)+j)=*(B+i);
	}
}
//紧凑格式
for(i=1;i<N;i++)
{       *(C+i*(N+1))=*(C+i*(N+1))/(*C);
             
	for(j=1;j<N+1;j++)
	{temp=0;
	if(i>j)
	{for(k=0;k<j;k++)
	{temp=temp+(*(C+k*(N+1)+j))*(*(C+i*(N+1)+k));}  
	*(C+i*(N+1)+j)=(*(C+i*(N+1)+j)-temp)/(*(C+j*(N+1)+j));}	
	else
	{for(k=0;k<i;k++)
	{temp=temp+(*(C+i*(N+1)+k))*(*(C+k*(N+1)+j));}
	*(C+i*(N+1)+j)=(*(C+i*(N+1)+j)-temp);}
	}
}
//求线性方程组的解
*(Root+N-1)=*(C+(N-1)*(N+1)+N)/(*(C+(N-1)*(N+1)+N-1));
for(i=N-2;i>=0;i--)
{temp=0;
for(k=N-1;k>i;k--)
{temp=temp+(*(C+i*(N+1)+k))*(*(Root+k));}
*(Root+i)=(*(C+i*(N+1)+N)-temp)/(*(C+i*(N+1)+i));
}
delete C;
}