Google ColabのJulia言語で2次元Swift-Hohenberg方程式（∂φ/∂t=φ - φ³ - (∇²+ k₀²)²φ、熱対流などを表す）を差分法で計算してパターンを動画にしてみる。

前回は複素TDGL方程式をやったので、今回は類似のSwift-Hohenberg方程式

∂φ/∂t = φ - φ³ - (∇²+ k₀²)²φ

をやってみよう。空間の4階微分が出てくる。コードはこんな感じで。

TDGLのプログラム流用なので無駄があり変数名がおかしいがまあ気にしないで。

using Plots

using Printf

using Random

function main()

    #パラメータ設定

    n = 128

    m = 1000

    ndiv = 2

    dt = 0.02

    dx = 1.0

    C1 = 1.0

    #初期設定

    X1 = zeros(n + 2, n + 2)

    Y1 = zeros(n + 2, n + 2)

    X2 = zeros(n + 2, n + 2)

    MT = MersenneTwister()

    Random.seed!(MT, 42)

    @inbounds for j in 2:(n + 1)

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

            X1[i, j] = 0.1 * (2.0 * rand(MT) - 1.0)

        end

    end

    # 結果を格納する配列

    results = []

    for t in 1:m

        #境界条件

        @inbounds for i in 2:(n + 1)

            X1[1, i] = X1[2, i]

            X1[n + 2, i] = X1[n + 1, i]

            X1[i, 1] = X1[i, 2]

            X1[i, n + 2] = X1[i, n + 1]

        end

        #一つ目のラプラシアン

        @inbounds for j in 2:(n + 1)

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

                lapX = (X1[i + 1, j] + X1[i - 1, j] + X1[i, j + 1] + X1[i, j - 1] - 4.0 * X1[i, j]) / (dx * dx)

                Y1[i, j] = lapX + C1 * C1 * X1[i, j]                                                                              

            end

        end

        #境界条件

        @inbounds for i in 2:(n + 1)    

            Y1[1, i] = Y1[2, i]

            Y1[n + 2, i] = Y1[n + 1, i]

            Y1[i, 1] = Y1[i, 2]

            Y1[i, n + 2] = Y1[i, n + 1]

        end

        #Swift-Hohenberg計算

        @inbounds for j in 2:(n + 1)

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

                lapY = (Y1[i + 1, j] + Y1[i - 1, j] + Y1[i, j + 1] + Y1[i, j - 1] - 4.0 * Y1[i, j]) / (dx * dx)

                X2[i, j] = X1[i, j] + dt * (X1[i, j] * (1.0 - X1[i, j] ^ 2) - (lapY + C1 * C1 * Y1[i, j]))                                                                            

            end

        end

        @inbounds for j in 2:(n + 1)

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

                X1[i, j] = X2[i, j]

            end

        end

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

        if t % ndiv == 0

            push!(results, copy(X1))

        end

    end

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

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

        heatmap(results[i], title="t = $(@sprintf("%.3f", i * ndiv * dt))", clim=(-1.2, 1.2), size = (950, 800))

    end

    gif(anim, "Swift-Hohenberg.gif", fps = 10)    

end

これでmain()とすると動画ができる。