Google ColabのJulia言語で2次元複素TDGL方程式(Complex Time-Dependent Ginzburg Landau equation, ∂W/∂t = (1+iC₀)W + (1+iC₁)∇²W - (1+iC₂)|W|²W)を計算してスパイラルパターンをアニメーションにしてみる。

2次元複素TDGL方程式(Complex Time-Dependent Ginzburg Landau equation)は昔、Excel VBAで計算したり遊んでいた。

こんな形の方程式

∂W/∂t = (1+iC₀)W + (1+iC₁)∇²W - (1+iC₂)|W|²W

相転移に使われるGL方程式を複素数に拡張したもの。パラメータを選ぶとスパイラルパターン（らせん）が出ることはよく知られている。

今回はJuliaでやってみよう。VBAを移植するのはとても簡単。

コードはこんな感じで。

using Plots

using Printf

using Random

function main()

    #パラメータ設定

    n = 128

    m = 10000

    ndiv = 10

    dt = 0.1

    dx = 1.0

    C0 = 0.0

    C1 = 0.0

    C2 = 1.0

    #初期設定

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

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

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

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

    MT = MersenneTwister()

    Random.seed!(MT, 42)

    @inbounds for i in (Int(round(n * 0.4)) + 1):(Int(round(n * 0.6)) + 1)

        @simd for j in (Int(round(n * 0.4)) + 1):(Int(round(n * 0.6)) + 1)

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

            Y1[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]

            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

        # TDGLの計算

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

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

                R2 = X1[i, j] * X1[i, j] + Y1[i, j] * Y1[i, j]

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

                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] - C0 * Y1[i, j]) + (lapX - C1 * lapY) - R2 * (X1[i, j] - C2 * Y1[i, j]))                                    

                Y2[i, j] = Y1[i, j] + dt * ((C0 * X1[i, j] + Y1[i, j]) + (C1 * lapX + lapY) - R2 * (C2 * X1[i, j] + Y1[i, j]))                                                                                

            end

        end

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

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

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

                Y1[i, j] = Y2[i, j]

            end

        end

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

        if t % ndiv == 0

            push!(results, copy(X1))

        end

    end

    # 計算結果をアニメーションで表示、heatmapのスケールを-1から1に固定

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

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

    end

    gif(anim, "TDGLequation.gif", fps = 10)    

end

main()

とすれば数分で計算できる。

結果はこちら。

らせんを描いている！