数値流体力学9～1次元拡散方程式～

##はじめに
　今までは、スカラー移流方程式を扱ってきました。今回からは、二階の偏微分方程式を扱っていきます。今回はその中でも基本的な、1次元拡散方程式について扱っていきます。

##シリーズ構成
####1. スカラー移流方程式
　1.1 [輸送速度が正の場合](https://qiita.com/vardia/items/3bc2b21e75a4b7557897)
　1.2 [輸送速度の符号が不明の時の線形問題](https://qiita.com/vardia/items/0c20be629343492d0cb8)
　1.3 [数値流束を利用した方法](https://qiita.com/vardia/items/a48b7dc9a6431e0847ed)
　1.4 [輸送速度が未知量であるとき(Burgers方程式)](https://qiita.com/vardia/items/2c2c6471deea37e11a98)
　1.5 [多次元への拡張(2次元スカラー移流方程式)](https://qiita.com/vardia/items/473192fadac80baefa89)
　1.6 [実践的な計算法(TVD方程式)](https://qiita.com/vardia/items/26d01a78adcdf2df125a)
####2. 移流方程式の時間積分法
　2.1 [時間陽解法について](https://qiita.com/vardia/items/7a72953183f065a7dd8c)
　2.2 [時間陰解法について(Crank-Nicholson法、近似LDU分解)](https://qiita.com/vardia/items/0986ac24d6c0748e92b9)
####3. 2階偏微分方程式
　3.1 1次元熱伝導方程式の数値計算法
　3.2 楕円方程式の数値計算法
####4. 圧縮性流れの数値計算法
　4.1 オイラー方程式の数値計算法
　4.2 MacCormack法による数値計算法
　4.3 TVD法による数値計算法

##陽的解法
　1次元拡散方程式は以下の式で表せます。

```math
\frac{\partial{q}}{\partial{t}}-\alpha\frac{\partial^2{q}}{\partial{x^2}}=0 \tag{1}
```
ここで、$\alpha$は拡散係数です。時間微分項にオイラーの陽的差分を、空間微分項に中心差分を適用すると、

```math
\frac{q^{n+1}-q^{n}}{\Delta{t}}=\frac{q_{j+1}-2q_{j}-q_{j-1}}{\Delta{x^2}}\tag{2}
```
これは、次のような簡単な前進スキームになります。

```math
q^{n+1}-q^{n}=q^{n}+\mu(q_{j+1}-2q_{j}-q_{j-1})\tag{3}
```
ここで、
$$\mu=\alpha\frac{\Delta{t}}{\Delta{x^2}}\tag{4}$$
である。移流方程式と同様に、von Nuemannの安定性解析を行うことができます。
$$g=1+\mu(e^{-i\theta}-2+e^{i\theta})\tag{5}$$
または、
$$g=1-2\mu(1-cos{\theta})\tag{6}$$
これから、位相誤差は生じないことがわかります。$|g|<1$であるためには、
$$-1\leqq1-2\mu(1-cos{\theta})\leqq\tag{7}1$$
$$\Leftrightarrow 0\leqq 2\mu sin^2\frac{\theta}{2} \leqq 1\tag{8}$$
任意の$\theta$について安定であるためには、
$$0 \leqq \mu \leqq \frac{1}{2} \tag{9}$$

 では、実際に1次元拡散方程式を数値計算してみましょう。今回は、初期条件・境界条件を以下のように定めます。今回は、拡散現象の中でも温度拡散について扱っていこうと思います。

```math
\left \{
\begin{array}{l}
境界条件:q(0,t)=q(1,t)=0, \\
初期条件:q(x)=sin(\pi x)
\end{array}
\right.
```
初期温度$sin(\pi x)$の棒の両端を$0℃$で冷やした時の熱拡散を計算することになります。

コードは以下の通りです。

```python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

def init(jmax):
    x = np.linspace(0, 1, jmax)
    q = np.sin(np.pi * x)
    return (x, q)

def do_computing(x, q, a, dt, dx, nmax):
    fig=plt.figure(figsize=(7,7), dpi=100)
    plt.rcParams["font.size"] = 22

    for n in range(1, nmax+1):
        qold = q.copy()

        for j in range(1, jmax-1):
            dq = a * dt * (qold[j+1] - 2.0 * qold[j] + qold[j-1]) / (dx ** 2)
            q[j] = qold[j] + dq
        
        q[0] = 0
        q[-1] = 0
        
        im = plt.plot(x, q, color='red', lw=2)
        ims.append(im)
    
    ani = animation.ArtistAnimation(fig,ims)
    plt.grid(color='black', linestyle='dashed', linewidth=0.5)
    plt.xlim([0,1])
    plt.xlabel('x')
    plt.ylabel('q')
    plt.show()

jmax = 21
nmax = 500

a = 1
dt = 0.01
dx = 0.2 
x, q = init(jmax)
ims = []

do_computing(x, q, a, dt, dx, nmax)
```
実行すると、熱拡散の様子が見られます。アニメーションの時間が長いので一部だけ掲載したいと思います。
![熱拡散_Trim.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2329221/08195791-6dd8-c7a2-4cac-0baf90d9e38e.gif)

##陰的解法(Crank-Nicholson法)
　陽的解法では、時間が1次精度となってしまいます。そこで、時間微分項にも中心差分を適用することにより時間も2次精度にすることを考えましょう。そこで、第8回で扱ったCrank-Nicholson法を1次元拡散方程式に適用しましょう。すると、

```math
\frac{q^{n+1}-q^{n}}{\Delta{t}}=\frac{\alpha}{2}\Bigl\{\frac{q^{n+1}_{j+1}-2q^{n+1}_{j}+q^{n+1}_{j-1}}{\Delta{x^2}}-\frac{q^{n}_{j+1}-2q^{n}_{j}+q^{n}_{j-1}}{\Delta{x^2}}\Bigr\}\tag{10}
```
式(10)を変形をすると、


```math
(\frac{2}{\nu}+2)q^{n+1}_{j}-(q^{n+1}_{j+1}+q^{n+1}_{j-1})=(\frac{2}{\nu}-2)q^{n}_{j}+(q^{n}_{j+1}+q^{n}_{j-1})\tag{4}
```

もちろん、これを解くときは3項方程式の形をとります。

```math
\left[
\begin{array}{ccc}
\frac{2}{\nu}+2 & 1 & \cdots & \cdots & 0 \\
1 & \ddots & & & \vdots \\
\vdots & &\ddots & &\vdots \\
\vdots &  &  &\ddots & 1\\
0 & \cdots &\cdots & 1 & \frac{2}{\nu}+2
\end{array}
\right]
\left[
\begin{array}{ccc}
\vdots \\
q^{n+1}_{j-1} \\
q^{n+1}_{j} \\
q^{n+1}_{j+1} \\
\vdots
\end{array}
\right]
=
\left[
\begin{array}{ccc}
 \\
 \\
 \\
RHS \\
 \\
 \\
 \\
\end{array}
\right]
\tag{11}
```
RHSはRight Side Handの略で式(11)の右辺の係数ベクトルを表します。
　von Neumannの安定性解析もしておきましょう。

````math
g=\frac{1-\mu(1-cos\theta)}{1+\mu(1-cos\theta)}\tag{12}
```
$\theta=0$のとき$g=1$となります。もし。$\mu$が0よりも大きい場合、$g$が$\theta$とともに減少して、$\pi$のときに最小値をとります。したがって、安定限界は、

```math
\frac{1-2\mu}{1+2\mu}>-1\tag{13}
```
ですが、これはすべての$\mu$に対して成り立ちます。つまり、無条件で安定ということになります。ですが、計算が安定に進むということであり、その計算結果が正しいということを保証するものではありません。

　では、実際に計算してみましょう。

```python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

Nx = 100
Nt = 5000
delta_x = 1/(Nx-1)
delta_t = 0.05
count_T = []
count_x = []
ims = []

alpha = 0.1

TT = np.zeros([Nx, Nt])

#初期条件の設定
for i in range(Nx):
    TT[:,0] = np.sin(np.pi*i)

#境界条件の設定
for i in range(Nt):
    TT[0, i] = 0
    TT[-1, i] = 0

#熱伝達率の行列を作成    
coe = np.ones([Nx, Nt])
for i in range(Nx):
    coe[i,:] = alpha

#安定性を決める数　alpha*dt/dx^2　をもつ行列の作成
q = np.zeros([Nx, Nt])
vT = np.ones([Nx, Nt])
for i in range(Nx):
    vT[i,:] = delta_t * coe[i,:] / (delta_x ** 2)
    
fig = plt.figure(figsize=(7,7), dpi=100)
plt.rcParams["font.size"] = 18  
plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['ytick.direction'] = 'in'
plt.rcParams['legend.frameon'] = True
plt.rcParams['legend.edgecolor'] = 'black'

count_x.append(0)
count_T.append(0)
plt.plot(count_x, count_T, c='black', marker='o', lw=2, label='x=0.8[mm]')

#クランクニコルソン法の適用
for t in range(Nt-1):
    #連立方程式の係数行列Aの作成
    A = np.zeros([Nx-2,Nx-2])
    for i in range(Nx-2):
        A[i,i] = 2/vT[i,t] +2
        if i >=1 :
            A[i-1,i] = -1
        if i <= Nx-4 :
            A[i+1,i] = -1
    
    #係数行列bの作成
    b = np.zeros([Nx-2])
    for i in range(Nx-2):
        b[i] = TT[i,t]+(2/vT[i+1,t]-2) * TT[i+1,t] + TT[i+2,t] + q[i+1,t]
    b[0] += TT[0,t+1]
    b[Nx-3] += TT[-1,t+1]
    
    #Aq=bをqについて解く
    T = np.linalg.solve(A,b)
    for i in range(Nx-2):
        TT[i+1,t+1] = T[i]
        count_T.append(T[i])
        count_x.append(i)
        
    im = plt.plot(count_x, count_T, c='black', marker='o', lw=2)
    ims.append(im)

ani = animation.ArtistAnimation(fig,ims)     
plt.gca().spines['right'].set_visible(False)
plt.gca().spines['top'].set_visible(False)
plt.xlabel('x s')
plt.ylabel('T ℃')
plt.legend()
plt.show()
```

##まとめ
　今回は、二階偏微分方程式の数値計算初回として、1次元拡散方程式を扱いました。次回は、ポアソン方程式で有名な楕円方程式を扱います。

##参考文献
藤井孝蔵、立川智章、Pythonで学ぶ流体力学の数値計算法、2020/10/20