Google ColabのJulia言語で時空カオスを生む蔵本・シバシンスキー方程式(Kuramoto–Sivashinsky equation)∂u/∂t＋u∂u/∂x＋∂²u/∂²＋∂⁴/∂x⁴=0を空間6次の差分、時間は8次のルンゲクッタ法で計算して動画にしてみる。

今回は蔵本・シバシンスキー方程式。時間・空間的にカオスを生じる最も簡単な偏微分方程式と言われている。金平糖の方程式と呼ばれているとかいないとか。

∂u/∂t＋u∂u/∂x＋∂²u/∂²＋∂⁴/∂x⁴=0

これは何回もExcel VBAで計算しているのでそれを移植した。

空間方向は6次の差分（4階微分が出てくるので他の微分もそれに合わせた）でこんな感じで書ける。

1階：　(-u(x-3dx)/9 + u(x-2dx)-5u(x-dx)+5u(x+dx) -u(x+2dx)+u(x+3dx)/9)/(20dx/3)
2階：　(2u(x+3dx)-27u(x+2dx)+270u(x+dx)-490u(x)+270u(x-dx)-27u(x-2dx)+2u(x-3dx))/(180dx^2)
3階：　(-u(x+3dx)+8u(x+2dx)-13u(x+dx)+13u(x-dx)-8u(x-2dx)+u(x-3dx))/(8dx^3)
4階：　(-u(x+3dx)+12u(x+2dx)-39u(x+dx)+56u(x)-39u(x-dx)+12u(x-2dx)-u(x-3dx))/(6dx^4)

時間方向はRunge Kutta Felhbergの8次を使う。

コードはこんな感じ。

using Plots

using Random

using Printf

function main()

    n = 200 #空間分割数

    m = 400000 #時間分割数

    u = zeros(Float64, n + 6)

    ud = zeros(Float64, n + 6)

    xmax = 50.0

    xmin = 0.0

    x = zeros(Float64, n + 6)

    D = 0.15

    t = 0.0

    dt = 0.0005

    dx = (xmax - xmin) / (n - 3)

    MT = MersenneTwister()

    Random.seed!(MT, 42)

    #表示用

    results = []

    ndiv = 500

    #8次のルンゲクッタ法の係数

    a = zeros(Float64, 13)

    b = zeros(Float64, 13, 13)

    c = zeros(Float64, 13)

    ku = zeros(Float64, n, 13)

    a[1] = 0.0

    a[2] = 2.0 / 27.0

    a[3] = 1.0 / 9.0

    a[4] = 1.0 / 6.0

    a[5] = 5.0 / 12.0

    a[6] = 0.5

    a[7] = 5.0 / 6.0

    a[8] = 1.0 / 6.0

    a[9] = 2.0 / 3.0

    a[10] = 1.0 / 3.0

    a[11] = 1.0

    a[12] = 0.0

    a[13] = 1

    b[2, 1] = 2.0 / 27.0

    b[3, 1] = 1.0 / 36

    b[3, 2] = 1.0 / 12.0

    b[4, 1] = 1.0 / 24.0

    b[4, 3] = 1.0 / 8.0

    b[5, 1] = 5.0 / 12.0

    b[5, 3] = -25.0 / 16.0

    b[5, 4] = 25.0 / 16.0

    b[6, 1] = 1.0 / 20.0

    b[6, 4] = 1.0 / 4.0

    b[6, 5] = 1.0 / 5.0

    b[7, 1] = -25.0 / 108.0

    b[7, 4] = 125.0 / 108.0

    b[7, 5] = -65.0 / 27.0

    b[7, 6] = 125.0 / 54.0

    b[8, 1] = 31.0 / 300.0

    b[8, 5] = 61.0 / 225.0

    b[8, 6] = -2.0 / 9.0

    b[8, 7] = 13.0 / 900.0

    b[9, 1] = 2.0

    b[9, 4] = -53.0 / 6.0

    b[9, 5] = 704.0 / 45.0

    b[9, 6] = -107.0 / 9.0

    b[9, 7] = 67.0 / 90.0

    b[9, 8] = 3.0

    b[10, 1] = -91.0 / 108.0

    b[10, 4] = 23.0 / 108.0

    b[10, 5] = -976.0 / 135.0

    b[10, 6] = 311.0 / 54.0

    b[10, 7] = -19.0 / 60.0

    b[10, 8] = 17.0 / 6.0

    b[10, 9] = -1.0 / 12.0

    b[11, 1] = 2383.0 / 4100.0

    b[11, 4] = -341.0 / 164.0

    b[11, 5] = 4496.0 / 1025.0

    b[11, 6] = -301.0 / 82.0

    b[11, 7] = 2133.0 / 4100.0

    b[11, 8] = 45.0 / 82.0

    b[11, 9] = 45.0 / 164.0

    b[11, 10] = 18.0 / 41.0

    b[12, 1] = 3.0 / 205.0

    b[12, 6] = -6.0 / 41.0

    b[12, 7] = -3.0 / 205.0

    b[12, 8] = -3.0 / 41.0

    b[12, 9] = 3.0 / 41.0

    b[12, 10] = 6.0 / 41.0

    b[13, 1] = -1777.0 / 4100.0

    b[13, 4] = -341.0 / 164.0

    b[13, 5] = 4496.0 / 1025.0

    b[13, 6] = -289.0 / 82.0

    b[13, 7] = 2193.0 / 4100.0

    b[13, 8] = 51.0 / 82.0

    b[13, 9] = 33.0 / 164.0

    b[13, 10] = 12.0 / 41.0

    b[13, 12] = 1.0

    c[6] = 34.0 / 105.0

    c[7] = 9.0 / 35.0

    c[8] = 9.0 / 35.0

    c[9] = 9.0 / 280.0

    c[10] = 9.0 / 280.0

    c[12] = 41.0 / 840.0

    c[13] = 41.0 / 840.0

    #初期値

    for i in 1:(n + 6)

        x[i] = xmin + (i - 3) * dx

    end

    for i in 4:(n + 3)

        u[i] = 0.1 * (2.0 * rand(MT) - 1.0)

    end

    u[1] = u[n + 1]

    u[2] = u[n + 2]

    u[3] = u[n + 3]

    u[n + 4] = u[4]

    u[n + 5] = u[5]

    u[n + 6] = u[6]

    for i in 1:m

        # 8次のルンゲクッタ法計算

         @inbounds for ii in 4:(n + 3)

            ku[ii - 3, 1] = fu(ii, t + dt, u, dx, n) * dt

          end

        @inbounds for j in 2:13

            @inbounds @simd for ii in 4:(n + 3)

                ud[ii] = u[ii]

            end

            @inbounds for k in 1:(j - 1)

                @inbounds @simd for ii in 4:(n + 3)

                    ud[ii] = ud[ii] + b[j, k] * ku[ii - 3, k]

                end

            end

            @inbounds for ii in 4:(n + 3)

                ku[ii - 3, j] = fu(ii, t + a[j] * dt, ud, dx, n) * dt

            end

        end

        @inbounds @simd for j in 1:13

            for ii in 4:(n + 3)

                u[ii] = u[ii] + c[j] * ku[ii - 3, j]

            end

        end

        t = t + dt

        # ndivステップごとに結果を配列に格納

        if i % ndiv == 0

            push!(results, copy(u))

        end

    end

    # 計算結果をアニメーションで表示

    anim = @animate for i in 1:length(results)

        plot(x, results[i], ylim = (-3, 3), linewidth = 3,title="Kuramoto Sivashinsky equation", label="", xlabel="x", ylabel="u", size=(900,500))

    end

    gif(anim, "Kuramoto_Sivashinsky.gif", fps = 30)

end

 function fu(i, t, u, dx, n)

    u[1] = u[n + 1]

    u[2] = u[n + 2]

    u[3] = u[n + 3]

    u[n + 4] = u[4]

    u[n + 5] = u[5]

    u[n + 6] = u[6]

    d1 = (u[i + 3] / 9.0 - u[i + 2] + 5.0 * u[i + 1] - 5.0 * u[i - 1] + u[i - 2] - u[i - 3] / 9.0) / (20.0 * dx / 3.0)

    d2 = (2.0 * u[i + 3] - 27.0 * u[i + 2] + 270.0 * u[i + 1] - 490.0 * u[i] + 270.0 * u[i - 1] - 27.0 * u[i - 2] + 2 * u[i - 3]) / (180.0 * dx * dx)

    d4 = (-u[i + 3] + 12.0 * u[i + 2] - 39.0 * u[i + 1] + 56.0 * u[i] - 39.0 * u[i - 1] + 12.0 * u[i - 2] - u[i - 3]) / (6.0 * dx * dx * dx * dx)

    return -u[i] * d1 - d2 - d4

 end

main()で実行させるとこうなる。

全く予想できない動きをしている。