1次元波動方程式を差分法で解く（python)

シンプルな波動方程式

 \frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2} \\
0<x<l,0<t

これを解きたい

境界条件

u(0)=0,\frac{\partial u}{\partial x}\Big|_{x=l} = 0

左端は境界での$u$の値を指定するディリクレ境界条件で、右端は境界値の微分値を指定するノイマン境界条件

左は固定端、右は自由端という境界条件になっている。

初期条件

u(x,0) = f(x)\\
\frac{\partial u}{\partial t}=0

弦の振動や縦棒の振動であれば、初速は0ということ。

離散式

x_i = i\Delta x\\
t^n = n\Delta t\\
u_{i}^n = u(x_i,t^n)\\
(i = 1\cdots m, \Delta x = \frac{l}{m})

と置く。ただし$m$は空間分割数とする。

初期条件　境界条件

まず初期解を設定する。

u_{i}^0 = f(x_i)\\
u_{i}^1 = u_{i}^0

時間の一階微分が0なので$u_{i}^1 = u_{i}^0$となっている。

次に境界条件を設定する。
しかし、困ったことに右端のノイマン境界条件は空間微分で表されるため、表現するためには2要素使わなければならない。そのため、右端の境界にもう一つ偽の要素を加える。つまり、

u_{m+1}^n = u_{m}^n

とすれば、一階の空間微分が0であるという、右端のノイマン境界条件が満たされる。
そして、左端のディリクレ境界条件は、以下のように表せる。

u_1^n = 0

いよいよ方程式そのものを離散化していく。

空間微分の離散化

\frac{\partial^2 u}{\partial x^2}\Big|_{x=x_i}=\frac{u_{i+1}-2u_{i}+u_{i-1}}{\Delta x^2}

$u_{i+1}$と$u_{i-1}$を$x=x_i$でテーラー展開して求めるとこれを導けます。

時間発展の離散化

空間微分と同様に考えると，

\frac{\partial^2 u}{\partial t^2} =\frac{u^{n+1}-2u^{n}+u^{n-1}}{\Delta t^2}

であるから，

\begin{align}
u^{n+1}_i &= 2u^{n}_i-u^{n-1}_i + \Delta t^2 \frac{\partial^2 u}{\partial t^2}\Big|_{x=x_i}\\
&= 2u^{n}_i-u^{n-1}_i + c^2 \Delta t^2 \frac{u_{i+1}-2u_{i}+u_{i-1}}{\Delta x^2}
\end{align}

コードを書こう

サンプルクソコード

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from copy import deepcopy


fig, ax = plt.subplots()

c = 1.
x_left = 0.
x_right = 1.
t = 0.
# 101個の点で離散化を行うという意味
num_points = 101
dx = (x_right-x_left)/(num_points-1)
# ノイマン境界のために余分に一つ確保
x = np.linspace(x_left, x_right+dx, num_points+1)
# u = x/1000 で初期化
u = 1e-3*x
# ノイマン境界条件
u[-1] = u[-2]

u_before = deepcopy(u)
u_next = deepcopy(u)

# 一定時間ごとのuの値をここに積んでいく
u_at_t = [deepcopy(u)]


dt = 0.1 * dx/c

count = 0


while t < 0.5:
    print(np.diff(u,2))
    u_next[1:num_points] = 2.*u[1:num_points] - u_before[1:num_points] + c*c*dt*dt/(dx*dx) * np.diff(u, 2)

    # ノイマン境界条件
    u_next[-1] = u_next[-2]
    u_next[0] = 0.
    u_before = deepcopy(u)
    u = deepcopy(u_next)
    u_next = np.zeros(u.shape)
    t += dt
    count += 1
    # 16ステップに一回
    if count == 16:
        u_at_t.append(deepcopy(u))
        count = 0


def update(u, x):
    plt.clf()
    plt.xlim(x_left, x_right)
    plt.ylim(-1e-3, 1e-3)
    plt.plot(x, u)


# アニメーション表示するための関数
anim = FuncAnimation(fig, update, fargs=(
    x,), frames=u_at_t, interval=10)

plt.show()

# anim.save("foo.mp4", writer="ffmpeg",fps=10)

動画を保存するのは重たい作業になるので注意

続く

参考

動画作成のためのテクニック
https://qiita.com/msrks/items/e264872efa062c7d6955
差分式
https://qiita.com/Ushio/items/0249fd7a5363ccd914dd