📂동역학

로렌츠 96 시스템

모델 ¹

로렌츠-96 시스템^{Lorenz-96 system}은 $N$ 개의 변수 $x_{1}, x_{2}, \ldots, x_{N}$ 와 다음과 같은 상미분방정식으로 표현된다. $$ {\frac{ d x_{k} }{ dt }} = x_{k-1} \left( x_{k+1} - x_{k-2} \right) - x_{k} + F(t) \qquad k = 1 , \cdots , N \pmod{N} $$

여기서 $F(t)$ 는 외력^{external forcing}을 반영하는 항이다.

설명

로렌츠-96과 대비되는 의미에서 원래 널리 알려진 로렌츠 어트랙터는 로렌츠-63이라고도 불린다². 로렌츠-96은 로렌츠-63과 달리 자연수 $N$ 에 대한 제한 없이 정의될 수 있으므로 예측 기법 벤치마크 등에서 고차원이면서 캐어릭한 시스템으로써 널리 사용되고 있다.

바이퍼케이션

다음은 $N = 6$ 이고 $F(t) = A \sin(2t) + 2$ 인 경우에 $A$ 의 값이 변화함에 따라 시스템이 어떻게 달라지는지를 보여주는 바이퍼케이션 다이어그램이다.

코드

function factory_lorenz96(A::Number; ic = [0,0,1,1,1,1,1], tspan = [0., 10.], dt = 1e-4)
    function sys(v::AbstractVector)
        T, x1, x2, x3, x4, x5, x6 = v

        dT = 1
        dx1 = x6*(x2 - x5) - x1
        dx2 = x1*(x3 - x6) - x2
        dx3 = x2*(x4 - x1) - x3
        dx4 = x3*(x5 - x2) - x4
        dx5 = x4*(x6 - x3) - x5
        dx6 = x5*(x1 - x4) - x6
        return [dT; [dx1, dx2, dx3, dx4, dx5, dx6] .+ (A*sin(2T) + 2)]
    end

        
    t_ = first(tspan):dt:last(tspan)
    len_t_ = length(t_)
    
    t, tk = .0, 0
    v = ic; DIM = length(v)
    traj = zeros(len_t_+2, 2DIM)
    while tk ≤ len_t_
        t,_,_ = v
        v, dv = RK4(sys, v, dt)

        if t ≥ first(t_)
            tk += 1
            traj[tk+1,         1:DIM ] =  v
            traj[tk  , DIM .+ (1:DIM)] = dv
        end
    end
    return traj[2:(end-2), :]
end
factory_lorenz96(T::Type, args...; kargs...) =
DataFrame(factory_lorenz96(args...; kargs...), ["t", "x1", "x2", "x3", "x4", "x5", "x6", "dt", "dx1", "dx2", "dx3", "dx4", "dx5", "dx6"])

A_ = 1:0.01:4
@showprogress @threads for k in eachindex(A_)
    A = A_[k]
    traj = factory_lorenz96(DataFrame, A, tspan = [0., 2500], dt = 1e-3)
    _traj = traj[150001:end, :]
    CSV.write("G:/lorenz96/lorenz96_$(lpad(k, 4, '0')).csv", _traj)
end

function arglmax(x)
    bits = circshift(x, 1) .< x .> circshift(x, -1)
    bits[1] = false
    bits[end] = false
    return findall(bits)
end
scatter(repeat([3.2], length(_traj.x1[almx1])), _traj.x1[almx1])

plt_bfcn = plot(legend = :none)
@showprogress for k in eachindex(A_)
    _traj = CSV.read("G:/lorenz96SINDy/lorenz96SINDy_$(lpad(k, 4, '0')).csv", DataFrame)
    x1 = _traj.x1[_traj.t .> 1000]
    almx1 = arglmax(x1)
    scatter!(plt_bfcn, repeat([α_[k]], length(x1[almx1])), x1[almx1], color = :red, ms = 1, ma = 0.05)
end
plot(plt_bfcn, xticks = [-3.5, 0, 1, 4], xflip = true)
png("bifurcation_plot.png")