📂Julia

Juliaでリャプノフスペクトルを計算する方法

コード

Juliaでリャプノフスペクトルを計算する方法について見ていく。Juliaの基本的な微分方程式パッケージであるDifferentialEquations.jlと互換性のあるChaosTools.jlを使うので、扱いやすい¹。

上の図は、実際にローレンツアトラクターにおいて$\sigma \in [6, 15]$、$\rho \in [120, 150]$、$\beta \in [3, 5]$の範囲でリャプノフスペクトルを計算した結果である。

using DifferentialEquations, ChaosTools, ProgressMeter, Base.Threads, Plots

function sys(du, u, p, t)
    x, y, z = u; σ, ρ, β = p
    
    du[1] = σ*(y - x)
    du[2] = x*(ρ - z) - y
    du[3] = x*y - β*z
    return du
end

# sol = factory_lorenz63(DataFrame, [10, 28, 8/3])
# plot(sol.x, sol.y, sol.z, alpha = .5)

pm, pM = -2, 3; p0, p1 = 0, 1;
p_ = range(pm, pM, length = 2001)
σ_ = range(6, 15, length = 2001)
ρ_ = range(120, 150, length = 2001)
b_ = range(3, 5, length = 2001)
lpnv = callbfcn()
@showprogress @threads for k in eachindex(p_)
    lpnv[p_[k]] = lyapunovspectrum(CoupledODEs(sys, [100.0, 100, 100], [σ_[k], ρ_[k], b_[k]]), 100000)
end

• sys: 基本的にはDifferentialEquations.jlで使われる形式で微分方程式を定義すればよい。公式ドキュメントではSVectorを使うなど、やや不便な部分があるため、一般的な配列を使う方法に変えた。方程式が複雑になるほどこの方式のほうが便利であろうし、速度が非常に重要になる場合は、もともと提供されているガイドを参考にするのがよい。
• CoupledODEs: 方程式の次に、初期条件、パラメータが与えられる。
• lyapunovspectrum: リャプノフスペクトルを計算する関数である。反復回数を第二引数として受け取り、十分に大きな回数でなければ正確な結果は得られない。