第3-2回　TDMA法（Thomas’algorithm）[python]

• English (United States)
• 日本語

TDMA法とは何か？

係数行列が三重対角行列である方程式を解く

$$\left[\begin{array}{cccccc}b_0&c_0&&&& \\ a_0&b_0&c_1&&& \\ &\ddots& \ddots & \ddots && \\&& \ddots & \ddots & \ddots &\\ &&&a_{n-3}&b_{n-2}&c_{n-2} \\ &&&&a_{n-2} &b_{n-1} \end{array}\right] \left[\begin{array}{c} u_0 \\ u_1 \\ \vdots \\ \vdots \\u_{n-2} \\ u_{n-1} \end{array}\right] = \left[\begin{array}{c} d_0 \\ d_1 \\ \vdots \\ \vdots \\d_{n-2} \\ d_{n-1} \end{array}\right] $$

上記のようなn×n行列の三重対角行列を考えます。求めたいのはu行列です。
b,dはそれぞれn個あり、a,cはそれぞれn-1個あることに注意してください。

①　まず、行列の一行目の式は次のように書け、変形できます 。

$$b_{0}c_{0}+c_{0}u_{1}=d_1 ⇔ u_0+\frac{c_0}{b_0}u_1=\frac{d_1}{b_0}$$

ここで、\(e_0=\frac{c_0}{b_0} ,f_0=\frac{d_1}{b_0}\)と置き整理すると次の式を得ます。

$$u_0=f_0-e_{0}u_{1}$$

②　この一行前の結果を下記の二行目の行列に代入し整理します。

$$a_{0}u_{0}+b_{1}u_{1}+c_{1}u_{2}=d_1$$

$$⇔ a_{0}(f_{0}-e_{0}u_{1})+b_{1}u_{1}+c_{1}u_{2}=d_1$$

$$⇔ (b_{1}-a_{0}e_{0})u_{1}+c_{1}u_{2}=d_1-a_{0}f_{0}$$

$$⇔ u_{1}+\frac{c_{1}}{b_{1}-a_{0}e_{0}}u_{2}=\frac{d_{1}-a_{0}f_{0}}{b_{1}-a_{0}e_{0}}$$

ここで、\(e_{1}=\frac{c_{1}}{b_{1}-a_{0}e_{0}} ,f_1=\frac{d_{1}-a_{0}f_{0}}{b_{1}-a_{0}e_{0}}\)と置き整理すると次の式を得ます。

$$u_{1}=f_{1}-e_{2}u_{2}$$

③　②の動作を同様に繰り返すことで、行列の下から二段目では次の式が得られます。

$$u_{n-2}=f_{n-2}-e_{n-2}u_{n-1}$$

また、行列の一番下の段の式は次の通りです。

$$a_{n-2}u_{n-2}+b_{n-1}u_{n-1}=d_{n-1}$$

この二式を連立させ、次のように\(u_{n-1}\)を求めることができます。

$$a_{n-2}(f_{n-2}-e_{n-2}u_{n-1} )+b_{n-1}u_{n-1}=d_{n-1}$$

$$⇔ (b_{n-1}-a_{n-2}e_{n-2})u_{n-1}=d_{n-1}-a_{n-2}f_{n-2}$$

$$⇔ u_{n-1}=\frac{d_{n-1}-a_{n-2}f_{n-2}}{b_{n-1}-a_{n-2}e_{n-2}}$$

④　最後に②で作り出した下記の式に,n-2から0まで順番に代入を繰り返し行うことですべての\(u\)を求めることができます。

$$u_{i}=f_{i}-e_{i}u_{i+1} \hspace{12pt} (i=n-2,n-1,\cdots 0)$$

まとめ

結局のところ，この方法はガウスの消去法のように前進消去，後退代入を行っていますが，三重対角部分のみ計算を行うので計算量が少ないというわけです。

pythonを用いて、TDMA法のテストを行います。次のような行列を考えます。

$$\left[\begin{array}{ccc}3&1&0 \\ 1&4&2 \\ 0&2&5 \end{array}\right] \bf u = \left[\begin{array}{c} 5 \\ 15 \\ 19 \end{array}\right] $$

\(u= \left[\begin{array}{c} 1 \\ 2 \\ 3 \end{array}\right] \)となれば成功です。

この記事の最後に載せたpythonのコードで計算を行いました。結果は次の通りです。

[1.0, 2.0000000000000004, 2.9999999999999996]

上手くいきました。

解く順番は\(u_{n-1}\)から\(u_{0}\)となるので最後に\(u\)の順番を逆順にする『u.reverse()』を28行目に書いています。

次回はこの方法で定常移流拡散方程式を解いてみます。

by hide

import numpy

def TDMA(a, b, c, d):
    """
    a,cはn-1個の要素
    b,dはn個の要素
    """

    # ①
    e = /b[0]]
    f = [d[0]/b[0]]

    # ②
    for i in range(1, len(d)-1):
        e.append(c[i]/(b[i]-a[i-1]*e[i-1]))
        f.append((d[i]-a[i-1]*f[i-1])/(b[i]-a[0]*e[0]))

    # ③
    u = [(d[len(d)-1]-a[len(d)-2]*f[len(d)-2]) /
         (b[len(d)-1]-a[len(d)-2]*e[len(d)-2])]

    # ④
    for i in range(len(d)-2, -1, -1):
        u.append(f[i]-e[i]*u[len(d)-2-i])

    u.reverse()

    return u

#-----------------------------
#-- main                    --
#-----------------------------
a = [1.0, 2.0]
b = [3.0, 4.0, 5.0]
c = [1.0, 2.0]
d = [5.0, 15.0, 19.0]

kai = TDMA(a, b, c, d)
print(kai)