BOOST的LU分解求解线性方程以及求逆

原文网址：http://blog.sina.com.cn/s/blog_644dcdfc0100is9l.html

BOOST的LU分解求解线性方程以及求逆

1.LU分解

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/operation.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/lu.hpp>
int _tmain(int argc, _TCHAR* argv[])
{
 // 可不要与 std::vector 混在一起啊.
 // 可考虑 using namespace boost::numeric;
 // 然后用 ublas::vector
 using namespace boost::numeric::ublas;
 const int n = 4;
 permutation_matrix<double> P(n);
 matrix<double> m(n,n);
 vector<double> x(n); // 用来放解向量
 vector<double> v(n); // 常量向量

 v(0) = 5;
 v(1) = 3;
 v(2) = 2;
 v(3) = 2;
 m(0, 0) = 1;
 m(0, 1) = 2;
 m(0, 2) = 1;
 m(0, 3) = 1;
 m(1, 0) = 2;
 m(1, 1) = 1;
 m(1, 2) = 0;
 m(1, 3) = 0;
 m(2, 0) = 1;
 m(2, 1) = 0;
 m(2, 2) = 1;
 m(2, 3) = 0;
 m(3, 0) = 1;
 m(3, 1) = 0;
 m(3, 2) = 0;
 m(3, 3) = 1;

 try {

  lu_factorize(m,P);
  x = v;
  lu_substitute(m,P,x);
  std::cout << "解向量: " << x << std::endl;
 }
 catch(...) {
     std::cout << "奇异矩阵(无解或无穷解)." << std::endl;
 }

 getchar();
 return 0;
}

m=[1,2,1,1;2,1,0,0;1,0,1,0;1,0,0,1];

lu_factorize(m)的结果如下：

m=[1,2,1,1;2,-3,-2,-2;1,2/3,4/3,1/3;1,2/3,1/4,5/4];

其中L,U阵为

L=[1,0,0,0;2,1,0,0;1,2/3,1,0;1,2/3,1/4,1]; 即L=tril(m,-1)+eye(1);m的下三角阵与单位矩阵之和

U=[1,2,1,1;0,-3,-2,-2;0,0,4/3,1/3;0,0,0,4/5];即U=tril(m);m的上三角阵

可得m=L*U;

lu_factorize(m,P)的结果

P=[1,1,2,3];  注意P是一permutation_matrix

m=[2,1,0,0;1/2,3/2,1,1;1/2,-1/3,4/3,1/3;1/2,-1/3,1/4,5/4];

L=[1,0,0,0;1/2,1,0,0;1/2,-1/3,1,0;1/2,-1/3,1/4,1];

U=[2,1,0,0;0,3/2,1,1;0,0,4/3,1/3;0,0,0,5/4];

Pm=LU，P=[1,1,2,3]

Pm = A，A的0行为m的1行，A的1行为新的1行(m的第0行)，A的2行为m的2行，A的3的行为m的3行。