雑感等

音楽,数学,語学,その他に関するメモを記す.

波動方程式を有限差分法の陰解法で解く(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}_{??}の項(未来の項・未知数)を左辺に、それ以外を右辺に集めると↓。

f:id:kazmus:20210821152657p:plain

今回の境界条件

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

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

境界条件

  • 初速度0: 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()

調べていたら「分数階波動方程式」というのもあるらしく、その計算法もあった。しかし読み解く気力がない。 https://www.researchgate.net/publication/259623599_Numerical_Solution_of_Fractional_Wave_Equation_using_Crank-Nicholson_Method