Google ColabのJulia言語でソリトンを生み出すKdV方程式を計算してGIFアニメにしてみる。ZabuskyとKruskalが1965年に使った計算法（Leap frog法で非線形のところは周りとの平均値）を使ったが、当時は数日かかっていただろう計算が20秒ほどでできる。

今回はKdV方程式　∂u/∂t+u*∂u/∂x+δ²*∂³u/∂x³ = 0。

ZabuskyとKruskalのこの論文の方法を使う。Leap frog法で、非線形の部分に隣との平均を取るのがポイント。

Interaction of "Solitons" in a Collisionless Plasma and the Recurrence of Initial States

コードはこんな感じ。

using Plots

n = 202

m = 32000

c = 0.000484

dt = 0.0001

dx = 2.0 / (n - 2.0)

u1 = zeros(Float64, n + 2)

u2 = zeros(Float64, n + 2)

u3 = zeros(Float64, n + 2)

ϕ = zeros(Float64, n - 2, m) #表示用

for i in 3:n

    x = (i - 3.0) * dx

    u1[i] = cos(π * x)

end

u1[1] = u1[n - 1]

u1[2] = u1[n]

u1[n + 1] = u1[3]

u1[n + 2] = u1[4]

for i in 3:n

    u2[i] = u1[i] - (dt / dx / 2.0) * u1[i] * (u1[i + 1] - u1[i - 1]) - c * (dt / (dx * dx * dx) / 2.0) * (u1[i + 2] - 2 * u1[i + 1] + 2 * u1[i - 1] - u1[i - 2])

end

for i in 1:(n-2)

    ϕ[i, 1] = u1[i + 2]

    ϕ[i, 2] = u2[i + 2]

end

for j in 3:m

    for i in 3:n

        u3[i] = u1[i] - (dt / dx) * ((u2[i + 1] + u2[i] + u2[i - 1]) / 3.0) * (u2[i + 1] - u2[i - 1]) - c * (dt / (dx * dx * dx)) * (

            u2[i + 2] - 2 * u2[i + 1] + 2 * u2[i - 1] - u2[i - 2])

    end

    for i in 3:n

        u1[i] = u2[i]

        u2[i] = u3[i]

    end

    u1[1] = u1[n - 1]

    u1[2] = u1[n]

    u1[n + 1] = u1[3]

    u1[n + 2] = u1[4]

    u2[1] = u2[n - 1]

    u2[2] = u2[n]

    u2[n + 1] = u2[3]

    u2[n + 2] = u2[4]

    for i in 1:(n-2)

        ϕ[i, j] = u3[i + 2]

    end

end

x = zeros(Float64, n - 2)

for i in 1:(n-2)

    x[i] = dx * i

end

GIFアニメにするには、 

@gif for i in 1:50:m

    plot(x, ϕ[:, i], xlim=(0, 2), ylim=(-1, 3), title = "KdV Equation", label="", xlabel="x", ylabel="phi", size=(800,500))

end

でOK。結果はこちら。

過去のもの：

Google ColabのJulia言語でDifferentialEquationsパッケージを使ってピタゴラスの三体問題を35段14次ルンゲクッタFeagen法で計算し、Colab上でGIFアニメにしてみた。dtminを小さくしないと途中で止まってしまうことにハマった…  

 Google ColabのJulia言語でDifferentialEquationsパッケージを使って35段14次ルンゲクッタFeagen法（BigFloat使用）、オイラー法、4段4次のルンゲクッタ法、Tsit5、Dormand＆Princeの8次(全部Float64)でローレンツ方程式をアダプティブは切って固定刻み幅で計算して比較。

Julia言語でタッパーの自己言及式（不等式を計算して図示するとまた不等式になる）を描いてみる。543桁の数を含む計算が必要だが、デフォルトで任意精度演算が可能なので容易にできた。