📂줄리아

줄리아에서 랴푸노프 스펙트럼 계산하는 법

코드

줄리아에서 랴푸노프 스펙트럼을 계산하는 방법에 대해 알아보려고 한다. 줄리아의 기본적인 미분방정식 패키지인 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: 랴푸노프 스펙트럼을 계산하는 함수다. 반복횟수를 두 번째 인자로 받아주고, 충분히 큰 횟수가 되어야 정확한 결과를 얻을 수 있다.