c语言报告 列主元lu分解,列主元LU分解法 C++程序

// 本程序用来解四维线性方程组。当然，将在最开始将N换成其他数字、系数矩阵和右侧矩阵换掉即可解更多维的线性方程组

//在运行本程序前，一定要确保方程组有解(通过线性代数里学的知识加以判断)

#include

#include

#include

#define N 4

//作函数声明:

void DirectLU(double a[N][N+1],double []);

//定义函数作LU分解并求得最终解

void swap(double &,double &);

//交换函数,用来交换两个元素

int main()

{

double x[N];

int i,j;

double a[N][N+1]={

1,-1,2,-1,-8,

2,-2,3,-3,-20,

1,1,1,0,-2,

1,-1,4,3,4

};  //输入系数矩阵和右侧矩阵y放在一个二维数组里面

cout<

for(i=0;i

for(j=0;j

cout<

cout<

}

cout<

for(i=0;i

cout<

cout<

cout<

DirectLU(a,x);  //调用LU分解函数

cout<

for(i=0;i

for(j=0;j

cout<

cout<

}

cout<

为: \n";

for(i=0;i

cout<

cout<

cout<

cout<

for(i=0;i

cout<

"<

cout<

return 0;

}

void swap(double &a,double &b)

//交换函数

{

a=a+b;

b=a-b;

a=a-b;

}

void DirectLU(double a[N][N+1],double x[])  //列主元LU分解函数

{

int i,r,k,j;

double s[N],t[N];//={-20,8,14,-0.8};

double max;

for(r=0;r

{

max=0;

j=r;

for(i=r;i

{

s[i]=a[i][r];

for(k=0;k

s[i]-=a[i][k]*a[k][r];

s[i]=s[i]>0?s[i]:-s[i]; //s[i]取绝对值

if(s[i]>max){

j=i;

max=s[i];

}

}

if(j!=r) //选出的主元所在行j若不是r,则对两行相应元素进行调换

{

for(i=0;i

swap(a[r][i],a[j][i]);

}

for(i=r;i

for(k=0;k

a[r][i]-=a[r][k]*a[k][i];

}

for(i=r+1;i

//记算第r列的元素

{

for(k=0;k

a[i][r]-=a[i][k]*a[k][r];

a[i][r]/=a[r][r];

}

}

for(i=0;i

t[i]=a[i][N];

for(i=N-1;i>=0;i--) //利用回代法求最终解

{

for(r=N-1;r>i;r--)

t[i]-=a[i][r]*x[r];

x[i]=t[i]/a[i][i];

}

}