目次

連立方程式 (Simultaneous Equation)

連立方程式 (Simultaneous Equation)

数値計算ではそのモデルを表す方程式を離散化する.
このとき, 連立一次方程式の問題に帰着されることが多い.
連立一次方程式の解法としては, Gaussの消去法やGauss Seidel法などがある.

Gaussの消去法 (Gaussian Elimination)

Gaussの消去法は係数行列を掃き出し法により上三角行列に変形し, 後退代入により解を求めるものである.
このとき, 部分pivot選択をおこない, 丸め誤差の影響を小さくする.

#include <cmath>    // for std::fabs, std::pow
#include <iostream> // for std::cerr, std::cin, std::cout
#include <optional> // for std::optional
#include <string>   // for std::string
#include <vector>   // for std::vector
#include <utility>  // for std::move
#include <boost/format.hpp>          // for boost::format
#include <boost/multi_array.hpp>     // for boost::multi_array
#include <boost/range/algorithm.hpp> // for boost::range::copy

namespace{
    using dmatrix = boost::multi_array<double, 2>;
    using dvector = std::vector<double>;

    static auto constexpr MAXBUFSIZE = 32;
    static auto const eps = pow(2, -50);

    template <typename T, bool check_positive_num = true>
    T input_parameter(std::string const &str){
        T param;
        while (true){
		    std::cout << str;
            std::cin >> param;

		    std::cin.clear();
		    std::cin.ignore(MAXBUFSIZE, '\n');

            // 初期条件として<=0の場合に再入力が要らない時は<type, false>
            if constexpr (check_positive_num){
                if (!std::cin.fail() && param > static_cast<T>(0)){ 
			        break;
		        }
            } 
            else{
                if (!std::cin.fail()){
                    break;
                }
            }
	    }

        return param;
    }

    void display_result(dvector const &x);

    // 後退代入
    void Gauss_backward(dmatrix &A);

    // Gaussの消去法
    std::optional<dvector> Gauss_elimination(dmatrix &A, double eps);

    // 掃き出し法
    bool Gauss_forward(dmatrix &A, double eps);

    // pivot選択
    void Gauss_pivot(dmatrix &A, int k);

    dvector input_row(int size);

    dvector make_result(dmatrix const &A);    
}

int main(){
    auto const size = input_parameter<int> ("Enter the size of a expansion coefficient matrix n-rows and (n+1)-colowns\n");

    dmatrix A(boost::extents[size][size + 1]);
    
    std::cout << "Enter elements of the expansion coefficient matrix" << std::endl;
    for (int i = 0; i < size; i++){
        auto const row = input_row(size + 1);
        boost::range::copy(row, A[i].begin());
    }

    if (auto const res = Gauss_elimination(A, eps); res){
        auto const x(*res);
        display_result(x);
    }
    else{
        std::cerr << "Error! Gaussian Elimination failed" << std::endl;

        return -1;
    }
}

namespace{
    void display_result(dvector const &x){
        auto const size = static_cast<int>(x.size());
        for (int i = 0; i < size; i++){
            std::cout << "x[" << i << boost::format("] = %.7f\n") % x[i];
        }
    }

    void Gauss_backward(dmatrix &A){
        auto const size = static_cast<int>(A.size());

        A[size - 1][size] /= A[size - 1][size - 1];
        for (int i = size - 2; i >= 0; i--){
            auto ax = 0.0;
            for (int j = i + 1; j < size; j++){
                ax += A[i][j] * A[j][size];
            }

            A[i][size] = (A[i][size] - ax) / A[i][i];
        }
    }

    std::optional<dvector> Gauss_elimination(dmatrix &A, double eps){
        if (!Gauss_forward(A, eps)){
            return std::nullopt;
        }
        
        Gauss_backward(A);

        return std::make_optional(make_result(A));
    }

    bool Gauss_forward(dmatrix &A, double eps){
        auto const size = static_cast<int>(A.size());

        for (int k = 0; k < size - 1; k++) {
            Gauss_pivot(A, k);

            for (int i = k + 1; i < size; i++){
                for (int j = k; j < size + 1; j++){
                    A[i][j] -= A[i][k] * (A[k][j] / A[k][k]);
                    auto tmp = A[i][j];
                    auto a = 0;
                }
            }
        }

        if (std::fabs(A[size - 1][size - 1]) < eps){
            return false;
        }

        return true;
    }

    void Gauss_pivot(dmatrix &A, int k){
        auto const size = static_cast<int>(A.size());

        auto pivot_ren = k;
        auto pivot = std::fabs(A[k][k]);
        for (int i = k + 1; i < size; i++){
            if (std::fabs(A[i][k]) > pivot){
                pivot_ren = i;
                pivot = std::fabs(A[i][k]);
            }
        }

        if (k != pivot_ren){
            dvector tmp(A[k].begin(), A[k].end());
            boost::range::copy(A[pivot_ren], A[k].begin());
            boost::range::copy(tmp, A[pivot_ren].begin());
        }
    }

    dvector input_row(int size){
        dvector row(size);

        auto success = true;
        do{
            auto i = 0;
            for (; i < size; i++) {
                std::cin >> row[i];

                if (std::cin.fail()) {
                    success = false;
                    break;
                }
            }
            if (i == size) {
                success = true;
            }

            std::cin.clear();
            std::cin.ignore(MAXBUFSIZE, '\n');
        } while (!success);

        return row;
    }

    dvector make_result(dmatrix const &A){
        auto const size = static_cast<int>(A.size());

        dvector x(size);

        for (int i = 0; i < size; i++){
            x[i] = A[i][size];
        }

        return x;
    }
}

Gauss Seidel法 (Gauss Seidel Method)

Gauss Seidel法は連立一次方程式の解を反復法によって求めるものである.

Gauss Seidel法は係数行列が正定値対称であるとき収束する.
また, 係数行列が狭義対角優位であるとき収束する.

#include <cassert>  // for assert
#include <cmath>    // for std::fabs, std::pow, std::sqrt
#include <iostream> // for std::cerr, std::cin, std::cout
#include <optional> // for std::optional
#include <string>   // for std::string
#include <utility>  // for std::move
#include <vector>   // for std::vector

namespace{
    static auto constexpr MAXBUFSIZE = 32;
    static auto constexpr MAXLOOP = 100000;
    static auto const eps = std::pow(2.0, -50);

    template <typename T, bool check_positive_num = true>
    T input_parameter(std::string const & str){
        T param;
        while (true){
		    std::cout << str;
            std::cin >> param;

		    std::cin.clear();
		    std::cin.ignore(MAXBUFSIZE, '\n');

            // 初期条件として<=0の場合に再入力が要らない時は<type, false>
            if constexpr (check_positive_num){
                if (!std::cin.fail() && param > static_cast<T>(0)){ 
			        break;
		        }
            } 
            else{
                if (!std::cin.fail()){
                    break;
                }
            }
	    }

        return param;
    }

    double calc_norm(std::vector<double> const &afterx, std::vector<double> const &beforex);

    std::optional<std::vector<double>> Gauss_Seidel(std::vector<std::vector<double>> const &A, std::vector<double> const &b, double eps);

    double inner_product(std::vector<double> const &a, std::vector<double> const &b);
}

int main(){
    auto const size = input_parameter<int> ("Enter the size of a coefficient matrix\n");

    // Ax = b 
    std::vector<std::vector<double>> A(size, (std::vector<double>(size)));
    std::vector<double> b(size);
    
    std::cout << "Enter elements of the coefficient matrix" << std::endl;
    for (int i = 0; i < size; i++){
        for (int j = 0; j < size; j++){
            A[i][j] = input_parameter<double, false>("");
        }
    }

    std::cout << "Enter elements of a right side matrix" << std::endl;
    for (int i = 0; i < size; i++){
        b[i] = input_parameter<double, false>("");
    }

    // Gauss Seidel methodが成功して解が収束していれば, resに解が入る
    if (auto const res = Gauss_Seidel(A, b, eps); res){
        // resから中身を取り出す
        auto const x = *res;

        for (auto i = 0; i < size; i++){
            std::cout << x[i] << std::endl;
        }
    }
    else{  
        // resに中身が入っていなかったら、解が収束していないのでエラーメッセージ
        std::cerr << "Error! Gauss_Seidel method failed!" << std::endl;

        return -1; 
    }
}

namespace{
    double calc_norm(std::vector<double> const &afterx, std::vector<double> const &beforex){
        auto const size = afterx.size();
        
        // afterxとbeforexのサイズが違っていたら強制終了
        assert(size == beforex.size());

        auto norm = 0.0;

        for (auto i =0ULL; i < size; i++){
            norm += (afterx[i] - beforex[i]) * (afterx[i] - beforex[i]);
        }

        return std::sqrt(norm);
    }

    std::optional<std::vector<double>> Gauss_Seidel(std::vector<std::vector<double>> const &A, std::vector<double> const &b, double eps){
        auto const size = A.size();

        // vect Aとvect bのサイズが違っていたら強制終了
        assert(size == b.size());

        std::vector<double> x(size, 0.0);

        for (auto i = 0; i < MAXLOOP; i++){
            auto const beforex(x);

            for (auto j =0ULL; j < size; j++){
                x[j] = (b[j] - (inner_product(A[j], x) - A[j][j] * x[j])) / A[j][j];
            }

            // 残差vectorのnormがEPS未満になったら終了
            if (calc_norm(x, beforex) < eps){
                // 収束したので、解ベクトルxをstd::moveで右辺値参照にキャストし、std::optionalでラップして返す
                return std::make_optional(std::move(x));
            }
        }

        return std::nullopt;
    }

    double inner_product(std::vector<double> const &a, std::vector<double> const &b){
        auto value = 0.0;

        auto const size = a.size();
        
        // vect aとvect bのサイズが違っていたら強制終了
        assert(size == b.size());

        for (auto i = 0ULL; i < size; i++){
            value += a[i] * b[i];
        }

    return value;
    }
}

TDMA法 (Tri-Diagonal Matrix Algorithm method)

執筆中

その他の解法

連立一次方程式を解く方法としては他にはCramerの公式(Cramer rule)やJacobi法(Jacobi rule)などがある.
Gauss Seidel法における計算で緩和係数を考えて調整するものをSOR法(Successive Over-Relaxation)という.
実際の数値計算ではCG法(conjugate gradient method)やCholesky分解(Cholesky decomposition), LU分解(LU decomposition)が主に用いられる.

Eigen

C++には行列計算ライブラリとしてEigenがある.
行列に関する計算は基本的にEigenを使う.

#include <iostream>    // for std::cout, std::endl;
#include <Eigen/Dense> // for Eigen::MatrixXd, Eigen::VectorXd, Eigen::ColPivHouseholderQR

int main(){
    Eigen::Matrix3d A;
    Eigen::Vector3d b;
    A << 2,3,1, 4,-5,-1, 6,-7,1;
    b << 7, 3, 5;

    Eigen::ColPivHouseholderQR<Eigen::Matrix3d> dec(A);
    Eigen::Vector3d x = dec.solve(b);

    std::cout << x << std::endl;
}

補遺 Eigenについて

EigenはC++のテンプレートを使用されて作られている線形代数ライブラリである.
Eigenの基本的なコマンドを紹介する.
ここで紹介していない便利なコマンドもたくさんあるので, 参考文献を適宜参考にしてほしい.

#include <iostream>    // for std::cout, std::endl
#include <Eigen/Dense> // for Eigen::

int main(){
    // 3*3 matrix A, Bと3*1 vector cの宣言と初期化
    Eigen::Matrix3d A;
    Eigen::Matrix3d B;
    Eigen::Vector3d c;
    Eigen::Vector3d d;
    Eigen::Matrix3d id = Eigen::Matrix3d::Identity(3, 3); // 単位行列で初期化
    A << 2,3,1, 4,-5,-1, 6,-7,1;
    B << 1,3,4, -5,-6,3, 6,-1,0;
    c << 7, 3, 5;
    d << -3, 1, 8;

    // matrix A, vecotr c, matrix idの出力
    std::cout << A << std::endl;
    std::cout << c << std::endl;
    std::cout << id << std::endl;

    // 要素へのアクセス(0行0列から始まることに注意)
    std::cout << A(2, 2) << std::endl;        // matrix Aの(2,2)成分
    std::cout << A.row(0) << std::endl;       // matrix Aの0行目
    std::cout << B.col(1) << std::endl;       // matrix Bの1列目
    std::cout << c.segment(0,2) << std::endl; //vector cの0要素目から1要素まで
    
    // 演算
    std::cout << A + B << std::endl; // matrix Aとmatrix Bの和  
    std::cout << A - B << std::endl; // matrix Aとmatrix Bの差
    std::cout << 2 * A << std::endl; // matrix Aの各要素を2倍
    std::cout << A * B << std::endl; // matrix Aとmatrix Bの積

    // 転置, 随伴, 逆行列
    std::cout << A.transpose() << std::endl; // matrix Aの転置行列
    std::cout << B.adjoint() << std::endl;   // matrix Bの随伴行列
    std::cout << A.inverse() << std::endl;   // matrix Aの逆行列

    // vectorの内積と外積
    std::cout << c.dot(d) << std::endl;   // vector cとvector dの内積
    std::cout << c.cross(d) << std::endl; //vector cとvector dの外積

    // Ax = cを解く
    Eigen::ColPivHouseholderQR <Eigen::Matrix3d> dec(A);
    Eigen::Vector3d x = dec.solve(c);
    std::cout << x << std::endl;

    // LU分解でAx = cを解く
    Eigen::FullPivLU <Eigen::Matrix3d> lu(A); //Eigen::PartialPivLUでも良い
    Eigen::Vector3d y = lu.solve(c);
    std::cout << y << std::endl;

    // 固有値分解
    Eigen::EigenSolver <Eigen::Matrix3d> eigensolver(A);
    std::cout << eigensolver.eigenvalues() << std::endl;  // 固有値
    std::cout << eigensolver.eigenvectors() << std::endl; // 固有ベクトル
}