// 本程序用来解四维线性方程组。当然,将在最开始将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];
}
}