Google ColabのJulia言語で1次元のGray-Scottモデル(∂u/∂t=u²v-(F＋k)u＋Du∂²u/∂x²,∂v/∂t=-u²v＋F(1-v)＋Dv∂²v/∂x²)を計算してパルスが次々分裂する様子を見る。空間6次の差分、時間8次のルンゲクッタ法で計算。

今回はGray-Scott方程式を計算。2次元だと動物の模様のようなものができて面白いが今回は1次元でやってみよう。

∂u/∂t=u²v-(F＋k)u＋Du∂²u/∂x²
∂v/∂t=-u²v＋F(1-v)＋Dv∂²v/∂x²

コードはこんな感じ。

using Plots

using Printf

function main()

    n = 200 #空間分割数

    m = 2000000 #時間分割数

    #u, vの初期化

    u = zeros(Float64, n + 6)

    v = zeros(Float64, n + 6)

    ud = zeros(Float64, n + 6)

    vd = zeros(Float64, n + 6)

    #物理パラメータ

    xmax = 1.6

    xmin = 0.0

    t = 0.0

    du = 0.00001

    dv = 0.00002

    Fc = 0.04

    kc = 0.06075

    dt = 0.01

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

    x = zeros(Float64, n + 6)

    #表示用

    results = []

    ndiv = 2000

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

    a = zeros(Float64, 13)

    b = zeros(Float64, 13, 13)

    c = zeros(Float64, 13)

    ku = zeros(Float64, n, 13)

    kv = 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] = exp(-(x[i] - xmax / 2.0) ^ 2 * 1000.0)

        v[i] = 1.0 - u[i]

    end

    for i in 1:m

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

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

            ku[ii - 3, 1] = f(ii, t + dt, u, v, du, dv, Fc, kc, dx, n) * dt

            kv[ii - 3, 1] = g(ii, t + dt, u, v, du, dv, Fc, kc, dx, n) * dt

        end

        @inbounds for j in 2:13

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

                ud[ii] = u[ii]

                vd[ii] = v[ii]

            end

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

                @inbounds @simd for ii in 2:(n + 1)

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

                    vd[ii] = vd[ii] + b[j, k] * kv[ii - 3, k]

                end

            end

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

                ku[ii - 3, j] = f(ii, t + a[j] * dt, ud, vd, du, dv, Fc, kc, dx, n) * dt

                kv[ii - 3, j] = g(ii, t + a[j] * dt, ud, vd, du, dv, Fc, kc, dx, n) * dt

            end

        end

        @inbounds @simd for j in 1:13

            for ii in 2:(n + 1)

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

                v[ii] = v[ii] + c[j] * kv[ii - 3, j]

            end

        end

        t = t + dt

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

        if i % ndiv == 0

            push!(results, copy(u))

        end

    end

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

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

        plot(x, results[i], ylim = (0.0, 0.5), linewidth = 4,title="Gray-Scott model", label="", xlabel="x", ylabel="u", size=(900,500))

    end

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

end

function f(i, t, u, v, du, dv, Fc, kc, 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]

    d2 = (2 * 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)

    return du * d2 + u[i] * u[i] * v[i] - (Fc + kc) * u[i]

end

function g(i, t, u, v, du, dv, Fc, kc, dx, n)

    v[1] = v[n + 1]

    v[2] = v[n + 2]

    v[3] = v[n + 3]

    v[n + 4] = v[4]

    v[n + 5] = v[5]

    v[n + 6] = v[6]

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

    return dv * d2 - u[i] * u[i] * v[i] + Fc * (1.0 - v[i])

end

実行結果。パルスが次々分裂していくのが見える。