Thomas法による数値解析

Thomas法が使える条件

係数行列が3重対角行列である連立一次方程式に対してこのThomas法を用いることで効率よく解を求めることができます．

流れ

1. 前進消去法で係数aの部分を0にする
2. 前進消去法で係数bの部分を1にする
3. 後退代入によって解を求める

前進消去

前進消去のざっくりとした手法のイメージとしては上(1行目)にある行から順番に行列を整理していくイメージです．今回であれば係数aを0、係数bを1にする作業を上にある行から行っていくイメージです．

まずはじめに係数a,bは0,1に変化させた時に対応する係数c,dの変数名をc',d'と置くことにします．
なので前進消去後の目標とする行列の形はこんな感じ．

1行目の処理

1行目で考えることとしては$b_1$を1にすることなので
$c_1' = c_1/b_1$
$d_1' = d_1/b_1$

k行目

とりあえず2行目で考えます．(すみませんちょっと数式markdownで書くの大変になってきたので手書きで．．)

こんな感じの状態なので、この状態から係数aを0、bを1にすることを考えればいいので $(② - ①*a_2)/(b_2-c_1'*a_2)$となります． そのため$c_2'=c_2/(b_2 - c_1' * a_2 )$、$ d_2' = (d_2 - d_1' * a_2)/(b_2 - c_1' * a_2) $となります．

一般化して書くと

c_k'=c_k / (b_k - c_{k-1}' * a_k ) \\
d_k' = (d_k - d_{k-1}' * a_k)/(b_k - c_{k-1}' * a_k)

後退代入

上の項までで前進代入だ終了してこのような形に行列がなりました．

後退代入では名前の通り行列の右下の成分から計算を行っていきます．ここで前進消去によって対角成分が1になっていることが効いてきます．しっくり来ない場合は前進消去後の行列と解Xの積を実際に手で計算してみるとわかるかと思います．

x_n = d_n' \\ 
x_k = d_k - c_kx_{k+1} \;\;\;\; (1<=k<n)

プログラム

既知である行列bと解の行列は次元数が同じなのでメモリ削減のために解は配列bに格納しています．また係数c', d'はp,qに置き換わっています．

.cpp

void Thomas(double *a, double *b, double *c, vector<double> d){
    vector<double> p, q;
    p.push_back(c[0] / b[0]);
    q.push_back(d[0] / b[0]);
    int n = d.size();
    
    
    for(int i = 1; i < n; i++){
        if(i < n-1){
            double _p = c[i] / (b[i] - c[i-1] * a[i]);
            p.push_back(_p);
        }
        double _q = (b[i] - a[i] * q[i-1]) / (b[i] - p[i-1] * a[i]);
        q.push_back(_q);
    }
    
    b[n-1] = q[n-1];
    for(int i = n-2; i >= 0; i--){
        b[i] = q[i] - p[i]*b[i+1];
    }
    
    for(int i = 0; i < n; i++){
        cout << d[i] << endl;
    }
}

(*実際にちゃんと動いているかどうかはこのサイトを使って確かめてみるといいもしれません)

備考

今回紹介したThomas法は拡散シミュレーションなどで用いられる陰解法の中で連立方程式を解くために用いられる．

参考

https://www1.gifu-u.ac.jp/~tanaka/numerical_analysis.pdf
http://www.mri-jma.go.jp/Publish/Technical/DATA/VOL_47/47_023.pdf
http://nas2008.akita-pu.ac.jp/~ozawa/Lecture/Grad_Schl/HandOut10.pdf