Thomas法で解く三重対角行列

# Thomas法で解く三重対角行列
　
この記事は[アドカレ](https://qiita.com/advent-calendar/2023/hachinohe_kosen)に参加しています。

# まえがき
　この記事は[【S01】3次スプライン補間](https://qiita.com/metaphysical_bard/items/e558b005255d1767fc70)にて解説している内容をそのまま移植したものです。最後の方にあるコード例は追記したもので、二通りの方法で三重対角行列を解いているので、よければ。

# 三重対角行列を解くアルゴリズム
　Thomas法では以下のような行列を解くことができます。
```math
\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}$で割るだけでいいですね。演算後の値には$'$を付けています。
```math
\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}$で割るという操作になる。)
```math
\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}
```
これで上三角行列となりました。
```math
\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}'$を掛けて引く操作。)
```math
\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}''$となりました。
```math
\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行目の処理をします。
```math
\begin{align}
c_{0}'&=c_{0}/b_{0} \\
d_{0}'&=d_{0}/b_{0}
\end{align}
```
次に、上から順に$n$行目$(0<n\leqq m)$に対して以下の処理をします。
```math
\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$行目の値は既に分かりました。
```math
x_{m}=d_{m}''=d_{m}'
```
下から順に$n$行目$(m>n\geqq 1)$に対して以下の処理をします。
```math
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}'}$として、
```math
\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$となったので、上三角行列となりました。
```math
\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}$が求まります。
```math
x_{m}=d_{m}''=d_{m}'/b_{m}'
```
$x_{m}$が求まったので、下から順に$c_{n}=0$となるようにしていきます。
```math
\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$となるように割ってあげます。
```math
\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$となるようにします。
```math
\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}'$となりました。
```math
\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)$に対して以下の処理をします。
```math
\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}$が求まります。
```math
x_{m}=d_{m}''=d_{m}'/b_{m}'
```
下から順に$n$行目$(m>n\geqq 1)$に対して以下の処理をします。
```math
x_{n}=d_{n}''=\frac{d_{n}'-d_{n+1}''\times c_{n}}{b_{n}'}
```

# プログラム例
　プログラム例では以下のような行列を解いています。
```math
\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法による数値解析](https://qiita.com/keito_takaishi/items/18f6943eb857b3360472)
・[第3-2回　TDMA法（Thomas’algorithm）[python]](https://teru-hide.com/cfd-python-tdma/)
・[Tridiagonal matrix algorithm](https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm)

# むすび
　Thomas法の紹介でした。