LUP分解法求解线性方程组

本节我们讨论如何用LUP分解法求解线性方程组，对于含有n个未知变量x1,x2,x3,…,xn的线性方程组:

同时满足方程组中所有方程的一个数值集：x1,x2,…,xn称为方程组的解。

将方程组改写成矩阵向量等式：

记为：

Ａx=b

如果A为非奇异矩阵，那么A存在逆矩阵，亦即方程组有解。

LUP分解的思想是找出三个n*n的矩阵，L、U和P，满足

PA=LU

其中，L是一个单位下三角矩阵，U是一个上三角矩阵，P是一个置换矩阵。则称满足上述等式的L、U和P为矩阵Ａ的LUP分解。

等式Ax=b两边同乘以P，变形为PAx=Pb，亦即：

LUx=Pb

设y=Ux,其中x就是我们要求的向量解。我们首先通过正向替换求解下三角系统：

Ly=Pb

得到未知变量y。然后，通过正向替换求解上三角系统：

Ux=y

得到未知变量x。

由于置换矩阵P是可逆的，等式PA=LU两边同乘以P^-1，得：

A＝P^-1LU

＝P^-1Ly

＝P^-1Pb

＝b

于是x就是Ax=b的解。

调用方法：

    linear::CLinearEqualtion l(3);
    l.init_A("A.txt");
    l.init_b("b.txt");
    l.lu_decomposition();
    l.lup_solve();
    l.show_y();
    l.show_x();
    l.save_solution();

Ｃ++实现：

#include <iostream>
#include <vector>
#include <cassert>
#include <fstream>

using namespace std;

namespace linear
{
#define array_1(__type) std::vector<__type>
#define array_2(__type) std::vector<array_1(__type)>

    class CLinearEqualtion
    {
        /*
        使用方法:
        1. 定义计算方程组的类对象,并初始化其系数矩阵的大小
        linear::CLinearEqualtion l(3);

        2. 读取系数阵 A
        l.init_A("A.txt");

        3. 读取 b
        l.init_b("b.txt");

        4. 对系数矩阵进行lu分解
        l.lu_decomposition();

        5. 求解方程组
        l.lup_solve();

        6. 显示反向替换的解
        l.show_y();

        7. 显示正向替换的解
        l.show_x();

        8. 保存方程的解
        l.save_solution();

        A.txt
        1    5    4
        2    0    3
        5    8    2

        b.txt
        12
        9
        5

        x[0] = -0.157895
        x[1] = -0.0526316
        x[2] = 3.10526
        */
    private:
        array_2(double) A;
        array_2(double) L;
        array_2(double) U;

        array_1(double) b;
        array_1(int) p;

        array_1(double) x;
        array_1(double) y;

    public:
        CLinearEqualtion(int n)
        {
            assert(n>1);

            A = array_2(double)(n);
            L = array_2(double)(n);
            U = array_2(double)(n);
            for (int i= 0; i < n;i++)
            {
                A[i] = array_1(double)(n);
                L[i] = array_1(double)(n);
                U[i] = array_1(double)(n);
            }
            b = array_1(double)(n);
            x = array_1(double)(n);
            y = array_1(double)(n);
            p = array_1(int)(n);
        } // !CLinearEqualtion(int n)
        virtual ~CLinearEqualtion(){} // !virtual ~CLinearEqualtion()

    public:
        void init_A(array_2(double) _A)
        {
            for (int i = 0; i < _A.size(); i++)
                A[i].assign(_A[i].begin(), _A[i].end());
        } // !void init_A(array_2(double) _A)

        void init_A(std::string _fileName)
        {
            std::ifstream in(_fileName.c_str());
            int i = 0;
            int j = 0;
            while(!in.eof())
            {
                double e = 0;
                in>>e;
                std::cout<<e<<std::endl;
                A[i][j] = e;
                if (j==A.size() - 1)
                {
                    i++;
                    j = 0;
                }
                else
                {
                    j++;
                }
            }
        } // !void init_A(std::string _fileName)

        void init_b(std::string _fileName)
        {
            std::ifstream in(_fileName.c_str());
            int i = 0;
            while(!in.eof())
            {
                double e = 0;
                in>>e;
                std::cout<<e<<std::endl;
                b[i++] = e;
            }
        } // !void init_b(std::string _fileName)

        void lu_decomposition()
        {
            int n = A.size(); // get array size
            for (int i = 0; i < n; i++)
                p[i] = i;
            for (int k = 0; k < n; k++)
            {
                int max = 0;
                int k_ = k;
                // get max e in col k.
                for (int i = k; i < n; i++)
                {
                    if (fabs(A[i][k]) > max)
                    {
                        max = fabs(A[i][k]);
                        k_ = i;
                    }
                }
                // make sure not all is zero.
                if (max == 0)
                    assert("singular matrix");

                // swap cur row,max row.
                if (k != k_)
                {
                    swap(p[k], p[k_]);
                    for (int i = 0; i < n; i++)
                        swap(A[k][i], A[k_][i]);                    
                }
                
                for (int i = k + 1; i < n; i++)
                {
                    // gen v
                    A[i][k] /= A[k][k];
                    for (int j = k + 1; j < n; j++)
                    {
                        A[i][j] -= A[i][k] * A[k][j];
                    }
                }
            }

            init_LU();
        } // !void lu_decomposition()

        void init_LU()
        {
            int n = A.size(); // get array size

            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    if (i > j)
                    {
                        L[i][j] = A[i][j];
                    } 
                    else
                    {
                        U[i][j] = A[i][j];
                        if (i == j)
                            L[i][i] = 1;
                    }
                }
            }

        } // !void init_LU()

        void lup_solve()
        {
            int n = A.size();
            int i = 0, j = 0;
            for (i = 0; i < n; i++)
            {
                x[i] = 0;
                y[i] = 0;
            }

            for (i = 0; i < n; i++)
            {
                y[i] = b[p[i]];
                for (j = 0; j < i; j++)
                    y[i] -= L[i][j] * y[j];
            }
            for (i = n - 1; i >= 0; i--)
            {
                double delt = y[i];
                for (j = n - 1; j > i; j--)
                    delt -= U[i][j] * x[j];
                x[i] = delt / U[i][j];
            }
        } // !void lup_solve()

        void show_y()
        {
            int n = A.size();
            std::cout << "###" << std::endl;
            for (int i = 0; i < n; i++)
            {
                std::cout << "y[" << i << "]=" << y[i] << std::endl;
            }
        } // !void show_y()

        void show_x()
        {
            int n = A.size();
            std::cout << "###" << std::endl;
            for (int i = 0; i < n; i++)
                std::cout << "x[" << i << "]=" << x[i] << std::endl;
        } // !void show_x()


        void save_solution()
        {
            int n = A.size();
            ofstream out("result.txt", ios::out);
            out << "-------------------------------------" << std::endl;
            std::cout << "-------------------------------------" << std::endl;
            for (int i = 0; i < n; i++)
            {                
                out << "x[" << i << "] = " << x[i]<< std::endl;
                std::cout << "x[" << i << "] = " << x[i]<< std::endl;
            }
            out << "-------------------------------------" << std::endl;
            std::cout << "-------------------------------------" << std::endl;
            out.close();
        }
    };
}

链接：http://pan.baidu.com/s/1hqJes4k 密码：elwz