LU分解双向代入法

LU分解法将系数矩A阵划为两个矩阵L与U的积，其中L为下三角，U为上三角，其中L的对角元素为1。

满足Ly=B，Ux=y。然后回代求解y，再回代解出x即可。

其代码如下：

---

#include<iostream.h> #include<math.h> #include<process.h> class lu { private: int i,j,k,n,*pivrow; double eps,pivot,sum,a1,*b,**a; public: void lu_input();//数据输入函数 void lu_decomposition();//分解与带入函数 void lu_output(); ~lu() { delete []b; delete []pivrow; for(i=0;i<n;i++) { delete[] a[i]; } delete []a; } }; void main() { lu solution; solution.lu_input();//输入系数矩阵 solution.lu_decomposition(); //执行LU分解，从前到后代入 solution.lu_output();//输出结果 } void lu::lu_input() { cout<<"输入方程的个数"; cin>>n; b=new double[n]; pivrow = new int[n]; a=new double*[n]; for(i=0;i<n;i++) { a[i] = new double[n]; } for(i=0;i<n;i++) //输入系数矩阵 for(j=0;j<n;j++) { cout<<"/n输入a["<<i<<"]["<<j<<"]="; cin>>a[i][j]; } for(i=0;i<n;i++) //输入右端向量的各元素 { cout<<"/n输入 b["<<i<<"]="; cin>>b[i]; } cout<<"/n输入最小主元素"; cin>>eps; //输入段结束 } void lu::lu_decomposition() { for(j=0;j<n;j++) { if(j>0) { for(i=0;i<j;i++) { if(i>0) { sum = 0.0; for(k=0;k<=(i-1);k++) { sum += a[i][k]*a[k][j]; } } else { sum = 0.0; } a[i][j] = a[i][j] - sum; } } pivot = 0.0; for(i=j;i<n;i++) { if(j>0) { sum = 0.0; for(k=0;k<=(i-1);k++) { sum += a[i][k]*a[k][j]; } } else { sum = 0.0; } a[i][j] = a[i][j] - sum; if(fabs(a[i][j])>pivot) { pivrow[j] = i; pivot = fabs(a[i][j]); } } if(pivot<eps) { cout<<"/n矩阵可能奇异.失败..."<<endl; exit(0); } if(pivrow[j]!=j) { for(k=0;k<n;k++) { a1 = a[pivrow[j]][k]; a[pivrow[j]][k] = a[j][k]; a[j][k] = a1; } } if(j!=(n-1)) { for(i=(j+1);i<n;i++) { a[i][j] /= a[j][j]; } } } for(i=0;i<n;i++) { sum = b[pivrow[i]]; b[pivrow[i]] = b[i]; if(i>0) { for(j=0;j<=(i-1);j++) { sum -= a[i][j] * b[j]; } } b[i] = sum; } for(i=(n-1);i>=0;i--) { sum = b[i]; if(i<(n-1)) { for(j=(i+1);j<n;j++) { sum -= a[i][j]*b[j]; } } b[i] = sum/a[i][i]; } } void lu::lu_output() { cout<<"/n结果是:"<<endl; for(i=0;i<n;i++) { cout<<"/nx["<<i<<"]="<<b[i]<<endl; } }