波動方程式を有限差分法の陰解法で解く（Julia）

波動方程式を↓の文献にある式で数値計算する。

https://www.uni-muenster.de/imperia/md/content/physik_tp/lectures/ws2016-2017/num_methods_i/wave.pdf

https://www.uni-muenster.de/Physik.TP/archive/fileadmin/lehre/part2_hypebolic/node20.html

記号

uは波の変位
uの右上の添え字（j-1, j, j+1）は離散時間の番号
uの右下の添え字（i-1, i, i+1）は離散空間の番号
Δtは離散時間の単位時間（微小時間）
Δxは離散空間の単位長さ（微小長さ）
cは波の伝播速度

式の意味

文献の式の右辺を少し変形すると↓になる。

f:id:kazmus:20210821145200p:plain

左辺は $\frac{\partial^{2} u(x,t)}{\partial t^{2}}$ を現在時刻（j）で中心差分に置き換えたもの。

右辺は $\frac{\partial^{2} u(x,t)}{\partial x^{2}}$ を、1微小時刻前（j-1）における中心差分に置き換えたものと、1微小時刻後（j+1）における中心差分に置き換えたものの平均。

現在の値を、過去と未来の平均に置き換えられる理由は理解していない。もしかしたら冒頭の文献にあるかもしれない。

（陽解法の場合は右辺は現在時刻（j）だけで表し、 $l.h.s. = c^{2} \cdot\frac{u^{j} _{i+1} - 2u^{j} _{i} + u^{j} _{i-1}}{\Delta x^{2}}$ ）

行列形式にしやすいように変形

さらに変形して↓になる。

f:id:kazmus:20210821150204p:plain

ここで、定数部分を↓としてまとめて表示する。

f:id:kazmus:20210821150237p:plain

さらに、 $u^{j+1}_{??}$ の項（未来の項・未知数）を左辺に、それ以外を右辺に集めると↓。

今回の境界条件

$u^{j}_{i}$ の添え字の範囲は

時刻 $j=0 \cdots J$
空間 $i= 0 \cdots I$

境界条件は

初速度０： $u^{-1}_{i} = u^{0}_{i}$
両端固定： $u^{j}_{-1} = u^{j}_{I+1}=0$

行列表示

境界条件を適用して、行列を使って整理すると

f:id:kazmus:20210821225341p:plain これを使ってJuliaで実装する。

Juliaプログラムソースと実行結果

実行結果↓ 横軸が空間方向、縦軸が時間方向、色が変位の大きさ。

f:id:kazmus:20210821230417p:plain 左右端で位相が反転し、固定端反射している。

図の下方（初期状態付近）だと孤立波の形を保っているが、時間が経過すると波形が崩れる。

ソース↓

using Plots
using WAV
using LinearAlgebra

const t_max = 3000  #時間点数
const x_max = 200   #空間点数

const dt = 1.0 / 500    #微小時間
const dx = 1.0 / x_max  #微小長さ
const c = 0.5   #伝播速度
const r = c^2 * dt^2 / (2 * dx^2)

function main()
    global t_max
    global x_max
    global dt
    global dx
    global r

    u_init = zeros(x_max + 1)   #初期条件
    u_curr = zeros(x_max + 1)   #U(j)
    u_prev = zeros(x_max + 1)   #U(j-1)
    u_hist = zeros(t_max + 1, x_max + 1)    #記録用

    #Mを定義
    M = zeros(x_max + 1, x_max + 1)
    for x0 = 1:x_max+1
        for x1 = 1:x_max+1
            M[x0, x1] = iseven(x0 + x1) ? 2 * r + 1 : -r
        end
    end
    M = SymTridiagonal(M)

    #初期変位（孤立波）
    w_len = 10#孤立波の幅
    w_pos = 20#孤立波の位置
    for x = 1:w_len+1
        u_init[x+w_pos] = 1 - cos(2 * pi * (x - 1) / w_len)
    end

    for x = 1:x_max+1
        u_curr[x] = u_init[x]
        u_prev[x] = u_init[x]
    end
    @time for t = 1:t_max+1
        u_next = M \ (2 * u_curr - M * u_prev)  #次ステップを計算

        #変数を更新
        for x = 1:x_max+1
            u_hist[t, x] = u_curr[x]
            u_prev[x] = u_curr[x]
            u_curr[x] = u_next[x]
        end
    end

    heatmap(u_hist)
    savefig("waveeq_implicit.png")
    println("end")
end
main()