# License

https://raw.githubusercontent.com/computational-sediment-hyd/tridiagonal-algorithm-matrix-of-periodic-boundary-conditions/master/LICENSE

# tridiagonal algorithm matrix of periodic boundary conditions

## はじめに

- 周期境界の三重対角は代表的な問題．
- その解法は，Sherman-Morrison formulaを使う．
- 参考(一部間違いあり)
:https://www.cfd-online.com/Wiki/Tridiagonal_matrix_algorithm_-_TDMA_%28Thomas_algorithm%29 


## Sherman-Morrison formulaとは
参考：https://mathtrain.jp/woodbury

- 逆行列の補助定理（Woodburyの恒等式）の特別な場合を示す．

逆行列の補助定理とは下式である．

$$
(A+BDC)^{-1}=A^{-1}-A^{-1}B(D^{-1}+CA^{-1}B)^{-1}CA^{-1}
$$

証明は省略する．

$
A:n\times n \,,\,
B:n\times k \,,\,
C:k\times n \,,\,
D:k\times k
$
が条件．

$k=1\,,\, D=1$のときは，

$$
(A+BC)^{-1}=A^{-1}-A^{-1}B(1+CA^{-1}B)^{-1}CA^{-1}
$$

であり，本式がシャーマンモリソンの公式と呼ばれる．

この式の良いところは，$A$の逆行列が計算できると$A+BC$の逆行列が計算できる点です．
今回の問題を例にとると，$A$にTDMAができれば，$A+BC$の逆行列が簡単に計算できるということ．





## 周期境界三重対角への適用

実際に計算してみます．

$$
A'x=d
$$

$$
A' =
\left(
\begin{array}{cccc}
b_1 & c_1 & 0 & a_1 \\
a_2 & b_2 & c_2 & 0 \\
0 & a_3 & b_3 & c_3 \\
c_4 & 0 & a_4 & b_4 \\
\end{array}
\right)
$$
とすると，

$$
A' = A + B \times C
$$

$$
A =
\left(
\begin{array}{cccc}
2 b_1 & c_1 & 0 & 0 \\
a_2 & b_2 & c_2 & 0 \\
0 & a_3 & b_3 & c_3 \\
0 & 0 & a_4 & b_4 + \dfrac{a_1 c_4}{b_1} \\
\end{array}
\right)
$$

$$
B =
\left(
\begin{array}{c}
-b_1  \\
0  \\
0  \\
c_4  \\
\end{array}
\right)
$$

$$
C =
\left(
\begin{array}{cccc}
1 & 0 & 0 & -\dfrac{a_1}{b_1}
\end{array}
\right)
$$

となるので，$A$にTDMAが適用できます

あとは，公式にあてはめて展開する.

$$
A y = d \,,\, A z = B \\
y = A^{-1}  d \,,\, z = A^{-1}  B
$$
とおくと

$$
x =  y - \dfrac{C \cdot y}{1 + C \cdot z}z
$$

となり，右辺第二項は$z$のスカラー倍のため，2回のTDMAで容易に計算可能である．

## 補足：TDMA（導出は省略）
 - 一般的にはメモリ節約版を使いますが（パタンカーの影響？），分かりにくいので通常版を使います．

三重対角行列
$$
a_i x_{i-1} + b_i x_i + c_i x_{i+1} = d_i
$$
について，
 
$$
\begin{align}
P_i &= -\frac{c_i}{a_i P_{i-1} + b_i} \\
Q_i &= \frac{-a_iQ_{i-1}+d_i}{a_i P_{i-1} + b_i} \\
P_1 &= -\frac{c_1}{b_1} \\
Q_1 &= \frac{d_1}{b_1}
\end{align}
$$

とすると，

$$
\begin{align}
x_n &= Q_n \\
x_i &= P_ix_{i+1} + Q_i
\end{align}
$$

## 実装

### Periodic boundary TDMA source code


```python
import numpy as np

def pTDMA(a,b,c,d):
    
    def TDMA(a,b,c,d):
        x, P, Q = np.empty_like(a), np.empty_like(a), np.empty_like(a)
        n = len(x)
        
        P[0], Q[0] = -c[0]/b[0], d[0]/b[0]
        for i in range(1, n):
            P[i] = -  c[i]                /(a[i]*P[i-1] + b[i])
            Q[i] =  (-a[i]*Q[i-1] + d[i])/(a[i]*P[i-1] + b[i])
            
        x[-1] = Q[-1]
        for i in range(n-2, -1, -1):
            x[i] = P[i]*x[i+1] + Q[i]
            
        return x
    
    a, b, c, d = np.array(a,dtype=float), np.array(b,dtype=float), np.array(c,dtype=float), np.array(d,dtype=float)
    a1, b1, c1, d1 = np.copy(a), np.copy(b), np.copy(c), np.copy(d) #deep copy
    a1[0], c1[-1] = 0.0, 0.0
    b1[0] = 2.0*b[0]
    b1[-1] = b[-1] + a[0]*c[-1]/b[0]
    d1[:] = 0.0
    d1[0] = -b[0]
    d1[-1] = c[-1]
    
    y = TDMA(a1,b1,c1,d)
    z = TDMA(a1,b1,c1,d1)
    
    w = ( y[0] - a[0]/b[0]*y[-1] ) / ( 1 + (z[0]- a[0]/b[0]*z[-1]) )
    x = y - w * z
    
    return x
```

### test


```python
n = 10
a = np.ones(n)
b = np.ones(n)
c = np.ones(n)
d = np.arange(1,n+1,1,dtype=float)

a *= -0.2
c *= 0.2

x = pTDMA(a,b,c,d)
print(x)
```

    [ 2.81456954  2.02649007  2.68211921  3.61589404  4.60264901  5.60264901
      6.58940397  7.65562914  8.31125828 11.09933775]


 - check print d


```python
print( b[0]*x[0] + c[0]*x[1] + a[0]*x[-1] )

for i in range(1, n-1):
    print( a[i]*x[i-1] + b[i]*x[i] + c[i]*x[i+1] )
    
print( a[-1]*x[-2] + b[-1]*x[-1] + c[-1]*x[0]  )
```

    1.0000000000000004
    2.0
    3.0
    3.9999999999999996
    5.0
    6.0
    7.0
    8.0
    9.0
    9.999999999999998


5c87ad0ed967bc69fe10ff8bee3e787ecd1c5949
ipynb
Jupyter Notebook
tridiagonal algorithm matrix of periodic boundary conditions.ipynb
computational-sediment-hyd/tridiagonal-algorithm-matrix-of-periodic-boundary-conditions
b3b756895b5a3ed13340c1d0eb307dced70613c4
MIT
tridiagonal algorithm matrix of periodic boundary conditions.ipynb
computational-sediment-hyd/tridiagonal-algorithm-matrix-of-periodic-boundary-conditions
b3b756895b5a3ed13340c1d0eb307dced70613c4
MIT
tridiagonal algorithm matrix of periodic boundary conditions.ipynb
computational-sediment-hyd/tridiagonal-algorithm-matrix-of-periodic-boundary-conditions
b3b756895b5a3ed13340c1d0eb307dced70613c4
MIT
Qwen/Qwen-72B
1. YES
2. YES