Google ColabのJulia言語でFPUT問題(Fermi–Pasta–Ulam–Tsingou、非線形結合した振動子が最初に与えたモードに戻る再帰現象)をDifferentialEquations.jlの2階の常微分方程式ソルバーDynamicalODEProblemでシンプレクティック8次のKahanLi8で計算、振動子の動きも動画にしてみる。

FPU問題というのがある。最近ではFPUT問題（Fermi–Pasta–Ulam–Tsingou）ともいわれる（TはManiacでコードを書いた女性プログラマ）。

https://en.wikipedia.org/wiki/Fermi%E2%80%93Pasta%E2%80%93Ulam%E2%80%93Tsingou_problem

非線形結合した振動子が最初に与えたモードに戻る再帰現象を発見した。

これを確かめてみよう。

まず、必要なパッケージをインストールする。が45分かかった…ここだけがGoogle ColabでJuliaをやるときのネック（Pythonは一瞬）。

using Pkg

Pkg.add("DifferentialEquations")

Pkg.add("RecursiveArrayTools")

DifferentialEquations.jlには2階の常微分方程式専用のソルバーがある。

Dynamical, Hamiltonian, and 2nd Order ODE Solvers

DynamicalODEProblem{isinplace}(f1, f2, v0, u0, tspan, p = NullParameters(); kwargs...)
SecondOrderODEProblem{isinplace}(f, du0, u0, tspan, p = NullParameters(); kwargs...)
HamiltonianProblem{T}(H, p0, q0, tspan, p = NullParameters(); kwargs...)

おお、HamiltonianProblem使えば一瞬じゃないか、と思ったら…そんな関数はないとエラー、じゃあSecondOrderODEProblemは…これもエラー、DynamicalODEProblemは…これは使える。

まあシンプレクティック積分を使うだろうな、と思ったらいろいろある。その中で

”For symplectic methods, higher order algorithms are the most efficient when higher accuracy is needed, and when less accuracy is needed lower order methods do better. Optimized efficiency methods take more steps and thus have more force calculations for the same order, but have smaller error. Thus, the “optimized efficiency” algorithms are recommended if your force calculation is not too sufficiency large, while the other methods are recommended when force calculations are really large (for example, like in MD simulations VelocityVerlet is very popular since it only requires one force calculation per timestep). A good go-to method would be McAte5, and a good high order choice is KahanLi8.
”

ということなのでシンプレクティック8次のKahanLi8を使う。

コードはこんな感じで。RecursiveArrayToolsを使うとsolverから望むデータが取り出せる（とGeminiに教えてもらった）。

using DifferentialEquations

using Plots

using RecursiveArrayTools

using Printf

n = 32 #振動子の数

tmax = 12000.0 #最大時間

dt = 0.2 #時間刻み幅

function Velocity!(dx, x, v, p, t)

    dx .= v

end

#　非線形な結合を設定

function Force!(dv, x, v, p, t)

    dv[1] = 0.0

    dv[n + 2] = 0.0

    for i in 2:(n+1)

        dv[i] = (x[i + 1] + x[i - 1] - 2.0 * x[i]) + 0.25 * ((x[i + 1] - x[i])^2 - (x[i] - x[i - 1])^2)

    end

end

# 初期条件

x0 = zeros(Float64, n + 2)

v0 = zeros(Float64, n + 2)

for i in 2:(n+1)

    x0[i] = sin(pi * (i - 1) / (n + 1))

end

# 2階常微分方程式をシンプレクティック8次で計算

prob = DynamicalODEProblem(Velocity!, Force!, x0, v0, (0.0, tmax))

sol = solve(prob, KahanLi8(), dt = dt)

#モードごとのエネルギーを計算

m = length(sol.t)

a = zeros(Float64, n)

ad = zeros(Float64, n)

E = zeros(Float64, n, m)

for i in 1:m

    for k in 1:n

        a[k] = 0.0

        ad[k] = 0.0

        for j in 2:(n+1)

            a[k] += sol[i].x[1][j] * sin((j - 1) * k * pi / (n + 1))

            ad[k] += sol[i].x[2][j] * sin((j - 1) * k * pi / (n + 1))

        end

        E[k, i] = sqrt(2 / (n + 1)) * ((ad[k] ^ 2) / 2 + 2 * (a[k] ^ 2) * (sin(pi * k / 2 / (n + 1)) ^ 2))

    end

end

#5つのモードのエネルギーをプロット

plot(dt * (1:m), E[1,:], label = "mode 1", xlabel = "time", ylabel = "Energy", legend = :topright, legendfontsize=6, ylim=(0.0, 0.4))

plot!(dt * (1:m),E[2,:], label = "mode 2")

plot!(dt * (1:m),E[3,:], label = "mode 3")

plot!(dt * (1:m),E[4,:], label = "mode 4")

plot!(dt * (1:m),E[5,:], label = "mode 5")

ではプロットした結果。

確かにt=10000くらいで最初のモードに戻ってる！

この振動子の動きを動画にしてみた。

@gif for i in 1:10:m

    plot(sol[i].x[1], xlim=(1, n), ylim=(-1.2, 1.2), title = "FPU problem", label="t=$(@sprintf("%.0f", i*dt))",legend=:bottomleft, xlabel="Position", ylabel="Amplitude", size=(800,500))

end

で動画にできる。

結果はこちら。確かにt=10000くらいで最初の振動に戻っている。