はじめに

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

シリーズ構成

1. スカラー移流方程式

　1.1 輸送速度が正の場合
　1.2 輸送速度の符号が不明の時の線形問題
　1.3 数値流束を利用した方法
　1.4 輸送速度が未知量であるとき(Burgers方程式)
　1.5 多次元への拡張(2次元スカラー移流方程式)
　1.6 実践的な計算法(TVD方程式)

2. 移流方程式の時間積分法

　2.1 時間陽解法について
　2.2 時間陰解法について(Crank-Nicholson法、近似LDU分解)

3. 2階偏微分方程式

　3.1 1次元熱伝導方程式の数値計算法
　3.2 楕円方程式の数値計算法

4. 圧縮性流れの数値計算法

　4.1 オイラー方程式の数値計算法
　4.2 MacCormack法による数値計算法
　4.3 TVD法による数値計算法

陽的解法

　1次元拡散方程式は以下の式で表せます。

\frac{\partial{q}}{\partial{t}}-\alpha\frac{\partial^2{q}}{\partial{x^2}}=0 \tag{1}

ここで、$\alpha$は拡散係数です。時間微分項にオイラーの陽的差分を、空間微分項に中心差分を適用すると、

\frac{q^{n+1}-q^{n}}{\Delta{t}}=\frac{q_{j+1}-2q_{j}-q_{j-1}}{\Delta{x^2}}\tag{2}

これは、次のような簡単な前進スキームになります。

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

\left \{
\begin{array}{l}
境界条件:q(0,t)=q(1,t)=0, \\
初期条件:q(x)=sin(\pi x)
\end{array}
\right.

初期温度$sin(\pi x)$の棒の両端を$0℃$で冷やした時の熱拡散を計算することになります。

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

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)

実行すると、熱拡散の様子が見られます。アニメーションの時間が長いので一部だけ掲載したいと思います。

陰的解法(Crank-Nicholson法)

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

\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)を変形をすると、

(\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項方程式の形をとります。

\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の安定性解析もしておきましょう。

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