Thomas法で解く三重対角行列

　
この記事はアドカレに参加しています。

まえがき

　この記事は【S01】3次スプライン補間にて解説している内容をそのまま移植したものです。最後の方にあるコード例は追記したもので、二通りの方法で三重対角行列を解いているので、よければ。

三重対角行列を解くアルゴリズム

　Thomas法では以下のような行列を解くことができます。

\begin{bmatrix}
b_{0} & c_{0} & \cdots & \cdots & \cdots & \cdots & \cdots & 0 \\
a_{1} & b_{1} & c_{1} & 0 & \cdots & \cdots & \cdots & 0\\
0 & a_{2} & b_{2} & c_{2} & 0 & \cdots & \cdots & 0\\
\vdots & 0 & a_{3} & b_{3} & c_{3} & 0 & \cdots & 0\\
\vdots && \ddots & \ddots & \ddots & \ddots & \ddots \\
0 & \cdots & \cdots & 0 & a_{m-2} & b_{m-2} & c_{m-2} & 0\\
0 & \cdots & \cdots & \cdots & 0 & a_{m-1} & b_{m-1} & c_{m-1} \\
0 & \cdots & \cdots & \cdots & \cdots & \cdots & a_{m} & b_{m}
\end{bmatrix}

\begin{bmatrix}
x_{0} \\ x_{1} \\ x_{2} \\ x_{3} \\
\vdots \\ x_{m-2} \\ x_{m-1} \\ x_{m} 
\end{bmatrix}
=
\begin{bmatrix}
d_{0} \\ d_{1} \\ d_{2} \\ d_{3} \\
\vdots \\ d_{m-2} \\ d_{m-1} \\ d_{m}
\end{bmatrix}

$0～m$のすべての自然数$n$で、$a_{n}=c_{n}=0$かつ$b_{n}=1$であれば$x_{n}=d_{n}$となるので、行列を解いたことになります。

Thomas法での行列の解き方

まずは$a_{n}=0$、$b_{n}=1$となることを目指します。

一行目は$b_{0}$で割るだけでいいですね。演算後の値には$'$を付けています。

\begin{align}
c_{0}'&=c_{0}/b_{0} \\
d_{0}'&=d_{0}/b_{0}
\end{align}

一行目以降の$n$行目に対しては、上の行から順番に$a_{n}=0$となるようにします。さらに、$b_{n}=1$になるようにします。(上の行に$a_{n}$を掛けて引いた後に、$b_{n}$で割るという操作になる。)

\begin{align}
a_{n}'&=\frac{a_{n}-b_{n-1}'\times a_{n}}{b_{n}-c_{n-1}'\times a_{n}}
=\frac{a_{n}-1\times a_{n}}{b_{n}-c_{n-1}'\times a_{n}}=\frac{0}{b_{n}-c_{n-1}'\times a_{n}}=0\\
b_{n}'&=\frac{b_{n}-c_{n-1}'\times a_{n}}{b_{n}-c_{n-1}'\times a_{n}}=1 \\
c_{n}'&=\frac{c_{n}-0 \times a_{n}}{b_{n}-c_{n-1}'\times a_{n}}
=\frac{c_{n}}{b_{n}-c_{n-1}'\times a_{n}}\\
d_{n}'&=\frac{d_{n}-d_{n-1}'\times a_{n}}{b_{n}-c_{n-1}'\times a_{n}}
\end{align}

これで上三角行列となりました。

\begin{align}
\begin{bmatrix}
1 & c_{0}' & \cdots & \cdots & \cdots & \cdots & \cdots & 0 \\
0 & 1 & c_{1}' & 0 & \cdots & \cdots & \cdots & 0\\
0 & 0 & 1 & c_{2}' & 0 & \cdots & \cdots & 0\\
\vdots & 0 & 0 & 1 & c_{3}' & 0 & \cdots & 0\\
\vdots && \ddots & \ddots & \ddots & \ddots & \ddots \\
0 & \cdots & \cdots & 0 & 0 & 1 & c_{m-2}' & 0\\
0 & \cdots & \cdots & \cdots & 0 & 0 & 1 & c_{m-1}' \\
0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0 & 1
\end{bmatrix}
\begin{bmatrix}x_{0} \\ x_{1} \\ x_{2} \\ x_{3} \\ \vdots \\ x_{m-2} \\ x_{m-1} \\ x_{m} \end{bmatrix} &= 
\begin{bmatrix}d_{0}' \\ d_{1}' \\ d_{2}' \\ d_{3}' \\ \vdots \\ d_{m-2}' \\ d_{m-1}' \\ d_{m}' \end{bmatrix}
\end{align}

今度は単位行列を目指すために$c_{n}=0$となるようにしていきます。最後の行ではすでに$x_{m}=d_{m}''=d_{m}'$となっているので、下の行から処理していきます。(下の行に$c_{n}'$を掛けて引く操作。)

\begin{align}
c_{n}''&=c_{n}'-b_{n+1}'\times c_{n}'=c_{n}'-1\times c_{n}'=0 \\
d_{n}''&=d_{n}'-d_{n+1}''\times c_{n}'
\end{align}

これで単位行列となったので、$0～m$の自然数$n$で$x_{n}=d_{n}''$となりました。

\begin{align}
\begin{bmatrix}
1 & 0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0 \\
0 & 1 & 0 & 0 & \cdots & \cdots & \cdots & 0\\
0 & 0 & 1 & 0 & 0 & \cdots & \cdots & 0\\
\vdots & 0 & 0 & 1 & 0 & 0 & \cdots & 0\\
\vdots && \ddots & \ddots & \ddots & \ddots & \ddots \\
0 & \cdots & \cdots & 0 & 0 & 1 & 0 & 0\\
0 & \cdots & \cdots & \cdots & 0 & 0 & 1 & 0 \\
0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0 & 1
\end{bmatrix}
\begin{bmatrix}x_{0} \\ x_{1} \\ x_{2} \\ x_{3} \\ \vdots \\ x_{m-2} \\ x_{m-1} \\ x_{m} \end{bmatrix} &= 
\begin{bmatrix}d_{0}'' \\ d_{1}'' \\ d_{2}'' \\ d_{3}'' \\ \vdots \\ d_{m-2}'' \\ d_{m-1}'' \\ d_{m}'' \end{bmatrix}
\end{align}

Thomas法まとめ

要点をまとめます。

まず、1行目の処理をします。

\begin{align}
c_{0}'&=c_{0}/b_{0} \\
d_{0}'&=d_{0}/b_{0}
\end{align}

次に、上から順に$n$行目$(0<n\leqq m)$に対して以下の処理をします。

\begin{align}
c_{n}'&=\frac{c_{n}}{b_{n}-c_{n-1}'\times a_{n}} \\
d_{n}'&=\frac{d_{n}-d_{n-1}'\times a_{n}}{b_{n}-c_{n-1}'\times a_{n}}
\end{align}

$m$行目の値は既に分かりました。

x_{m}=d_{m}''=d_{m}'

下から順に$n$行目$(m>n\geqq 1)$に対して以下の処理をします。

x_{n}=d_{n}''=d_{n}'-d_{n+1}''\times c_{n}'

Thomas法での行列の解き方(省メモリver.)

係数の分の配列を使用できるのであれば、こちらの方法の方が少ないメモリで処理することができます。

先ほどの方法との違いはいつ$b_{n}=1$とするかです。
先程の方法では、$a_{n}=0$、$b_{n}=1$、$c_{n}=0$の順で解いていきましたが、こちらの方法では$a_{n}=0$、$c_{n}=0$、$b_{n}=1$の順で解いていきます。

まずは一行目以外の$n$行目に対して、上の行から順番に$a_{n}=0$になるようにします。$w_{n}=\frac{a_{n}}{b_{n-1}'}$として、

\begin{align}
a_{n}'&=a_{n}-b_{n-1}'\times w_{n}=a_{n}-b_{n-1}' \times \frac{a_{n}}{b_{n-1}'}=0 \\
b_{n}'&=b_{n}-c_{n-1}\times w_{n} \\
c_{n}'&=c_{n}-0\times w_{n}=c_{n} \\
d_{n}'&=d_{n}-d_{n-1}'\times w_{n}
\end{align}

$a_{n}=0$となったので、上三角行列となりました。

\begin{align}
\begin{bmatrix}
b_{0} & c_{0} & \cdots & \cdots & \cdots & \cdots & \cdots & 0 \\
0 & b_{1}' & c_{1} & 0 & \cdots & \cdots & \cdots & 0\\
0 & 0 & b_{2}' & c_{2} & 0 & \cdots & \cdots & 0\\
\vdots & 0 & 0 & b_{3}' & c_{3} & 0 & \cdots & 0\\
\vdots && \ddots & \ddots & \ddots & \ddots & \ddots \\
0 & \cdots & \cdots & 0 & 0 & b_{m-2}' & c_{m-2} & 0\\
0 & \cdots & \cdots & \cdots & 0 & 0 & b_{m-1}' & c_{m-1} \\
0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0 & b_{m}'
\end{bmatrix}
\begin{bmatrix}x_{0} \\ x_{1} \\ x_{2} \\ x_{3} \\ \vdots \\ x_{m-2} \\ x_{m-1} \\ x_{m} \end{bmatrix} &= 
\begin{bmatrix}d_{0} \\ d_{1}' \\ d_{2}' \\ d_{3}' \\ \vdots \\ d_{m-2}' \\ d_{m-1}' \\ d_{m}' \end{bmatrix}
\end{align}

このとき、最終行では$b_{m}'$で割ってあげると$x_{m}$が求まります。

x_{m}=d_{m}''=d_{m}'/b_{m}'

$x_{m}$が求まったので、下から順に$c_{n}=0$となるようにしていきます。

\begin{align}
b_{n}''&=b_{n}'-0\times \frac{c_{n}'}{b_{n+1}'}=b_{n}' \\
c_{n}''&=c_{n}-b_{n+1}''\times \frac{c_{n}}{b_{n+1}''}=0 \\
d_{n}''&=d_{n}'-d_{n+1}''\times \frac{c_{n}'}{b_{n+1}'}
\end{align}

同時に$b_{n}=1$となるように割ってあげます。

\begin{align}
b_{n}''&=b_{n}'/b_{n}'=1 \\
d_{n}''&=\frac{d_{n}'-d_{n+1}''\times \frac{c_{n}}{b_{n+1}'}}{b_{n}'}
\end{align}

一行目も同様に$c_{n}=0$、$b_{n}=1$となるようにします。

\begin{align}
b_{0}'&=b_{0}/b_{0}=1 \\
d_{0}'&=\frac{d_{0}-d_{1}''\times \frac{c_{0}}{b_{1}'}}{b_{0}}
\end{align}

これで単位行列となったので、$1～m$の自然数$n$で$x_{n}=d_{n}''$、一行目で$x_{0}=d_{n}'$となりました。

\begin{align}
\begin{bmatrix}
1 & 0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0 \\
0 & 1 & 0 & 0 & \cdots & \cdots & \cdots & 0\\
0 & 0 & 1 & 0 & 0 & \cdots & \cdots & 0\\
\vdots & 0 & 0 & 1 & 0 & 0 & \cdots & 0\\
\vdots && \ddots & \ddots & \ddots & \ddots & \ddots \\
0 & \cdots & \cdots & 0 & 0 & 1 & 0 & 0\\
0 & \cdots & \cdots & \cdots & 0 & 0 & 1 & 0 \\
0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0 & 1
\end{bmatrix}
\begin{bmatrix}x_{0} \\ x_{1} \\ x_{2} \\ x_{3} \\ \vdots \\ x_{m-2} \\ x_{m-1} \\ x_{m} \end{bmatrix} &= 
\begin{bmatrix}d_{0}' \\ d_{1}'' \\ d_{2}'' \\ d_{3}'' \\ \vdots \\ d_{m-2}'' \\ d_{m-1}'' \\ d_{m}'' \end{bmatrix}
\end{align}

Thomas法まとめ(省メモリver.)

要点をまとめます。

まず、上から順に$n$行目$(0<n\leqq m)$に対して以下の処理をします。

\begin{align}
w_{n}&=\frac{a_{n}}{b_{n-1}'} \\
b_{n}'&=b_{n}-c_{n-1}\times w_{n} \\
d_{n}'&=d_{n}-d_{n-1}'\times w_{n}
\end{align}

最終行だけ$b_{m}'$で割ることで、$x_{m}$が求まります。

x_{m}=d_{m}''=d_{m}'/b_{m}'

下から順に$n$行目$(m>n\geqq 1)$に対して以下の処理をします。

x_{n}=d_{n}''=\frac{d_{n}'-d_{n+1}''\times c_{n}}{b_{n}'}

プログラム例

　プログラム例では以下のような行列を解いています。

\begin{align}
\begin{bmatrix}
1 & 2 & 0 & 0 & 0 & 0  \\
3 & 4 & 5 & 0 & 0 & 0 \\
0 & 6 & 7 & 8 & 0 & 0 \\
0 & 0 & 9 & 10 & 11 & 0 \\
0 & 0 & 0 & 12 & 13 & 14 \\
0 & 0 & 0 & 0 & 15 & 16 
\end{bmatrix}
\begin{bmatrix}x_{0} \\ x_{1} \\ x_{2} \\ x_{3} \\ x_{4} \\ x_{5} \\ x_{6} \end{bmatrix} &= 
\begin{bmatrix} 5 \\ 26 \\ 65 \\ 122 \\ 197 \\ 171 \end{bmatrix}
\end{align}

$x_{n}$の解は$n+1$となります。
TDMA1関数は係数の値を書き換えずに解を求めます。TDMA2関数は係数の値を書き換えますが、省メモリです。

TDMA_M.hpp

/*

TDMA_M.hpp

*/

#include <vector>

void TDMA1(int n, std::vector<double>& a, std::vector<double>& b, std::vector<double>& c, std::vector<double>& d) {
	std::vector<double> p(n), q(n);
	p[0] = c[0] / b[0];
	q[0] = d[0] / b[0];

	for (int i = 1; i < n; i++) {
		double w = b[i] - p[i - 1] * a[i];
		p[i] = c[i] / w;
		q[i] = (d[i] - q[i - 1] * a[i]) / w;
	}

	d[n - 1] = q[n - 1];
	for (int i = n - 2; i >= 0; i--)
		d[i] = q[i] - d[i + 1] * p[i];
}

void TDMA2(int n, std::vector<double>& a, std::vector<double>& b, std::vector<double>& c, std::vector<double>& d) {
	for (int i = 1; i < n; i++) {
		double w = a[i] / b[i - 1];
		b[i] = b[i] - c[i - 1] * w;
		d[i] = d[i] - d[i - 1] * w;
	}

	d[n - 1] = d[n - 1] / b[n - 1];
	for (int i = n - 2; i >= 0; i--)
		d[i] = (d[i] - d[i + 1] * c[i]) / b[i];
}

main.cpp

/*

main.cpp

*/

#include <iostream>
#include <vector>
#include <windows.h>
#include "TDMA_M.hpp"

int main() {
	std::cout << "start" << std::endl;

	std::vector<double> a = { 0, 3, 6, 9, 12, 15 };
	std::vector<double> b = { 1, 4, 7, 10, 13, 16 };
	std::vector<double> c = { 2, 5, 8, 11, 14, 0 };
	std::vector<double> d = { 5, 26, 65, 122, 197, 171 };

	std::cout << "TDMA1" << std::endl;
	TDMA1(d.size(), a, b, c, d);
	for (int i = 0; i < d.size(); i++)
		std::cout << "d[" << i << "] = " << d[i] << std::endl;

	d = { 5, 26, 65, 122, 197, 171 };
	std::cout << "TDMA2" << std::endl;
	TDMA2(d.size(), a, b, c, d);
	for (int i = 0; i < d.size(); i++)
		std::cout << "d[" << i << "] = " << d[i] << std::endl;

	system("pause");
	return 0;
}

参考文献

・Thomas法による数値解析
・第3-2回　TDMA法（Thomas’algorithm）[python]
・Tridiagonal matrix algorithm

むすび

　Thomas法の紹介でした。