logo

미분방정식으로 표현되는 시스템의 랴푸노프 스펙트럼과 그 수치적 계산법 📂동역학

미분방정식으로 표현되는 시스템의 랴푸노프 스펙트럼과 그 수치적 계산법

알고리즘

공간 X=RnX = \mathbb{R}^{n} 와 함수 f:X→Xf : X \to X 에 대해 다음과 같은 벡터필드가 미분 방정식으로 주어져 있다고 하자. x˙=f(x) \dot{x} = f(x) ff 의 자코비안 행렬 JJ 에 대해 다음과 같이 변분 방정식을 정의하자. Y˙=JY \dot{Y} = J Y 여기서 행렬함수 Y=Y(t)∈Rn×nY = Y(t) \in \mathbb{R}^{n \times n} 의 초기 조건은 항등행렬 Y(0)=InY(0) = I_{n} 로 둔다. 시스템 x˙=f(x)\dot{x} = f(x) 의 랴푸노프 스펙트럼 λ1,⋯ ,λn\lambda_{1} , \cdots , \lambda_{n} 은 x˙=f(x)\dot{x} = f(x) 의 수치해(이를테면 RK4 등으로 계산된) {xk}k=0K\left\{ x_{k} \right\}_{k=0}^{K} 에 대해 변분 방정식을 반복해서 풀며 근사적으로 계산된다. J(x)J(x) 는 원래 시스템의 트래젝터리 x(t)x(t) 에 종속적이다.

  • 여기서 변분 방정식 Y˙=JY\dot{Y} = J Y 를 푼다는 것은 퍼터베이션perterbation이라 불리는 행렬 혹은 벡터 YY 를 수치해를 얻을 때 사용했던 RK4와 같은 메서드로 푼다는 것이다. xkx_{k} 를 얻기 위해 사요했던 충분히 작은 시간 간격이 h≈0h \approx 0 이라면 변분 방정식도 같은 시간 간격 hh 를 사용해야한다.
  • 행렬 JJ 와 UU 에 대해 룽게-쿠타 메소드 RK4(J,U,h)RK4(J, U, h) 는 다음과 같이 정의된다. U(t+h)≈RK4(J,U,h)=h6(V1+2V2+2V3+V4)V1=JUV2=J(U+h2V1)V3=J(U+h2V2)V4=J(U+hV3) \begin{align*} U(t + h) \approx RK4(J, U, h) =& {\frac{ h }{ 6 }} \left( V_{1} + 2 V_{2} + 2 V_{3} + V_{4} \right)\\ V_{1} =& J U \\ V_{2} =& J \left( U + {\frac{ h }{ 2 }} V_{1} \right) \\ V_{3} =& J \left( U + {\frac{ h }{ 2 }} V_{2} \right) \\ V_{4} =& J \left( U + h V_{3} \right) \end{align*}

이제 기하적인 정의에 입각해서 그램-슈미트 직교화를 사용하는 전통적 방식과 QR 분해를 사용하는 방식을 소개한다.

그램-슈미트 직교화 기반 1

  • 행렬 UU 의 jj번째 열벡터를 U:jU_{:j} 와 같이 나타내자.
  • 행렬 UU 의 전치행렬을 UTU^{T} 와 같이 나타내자.
  • 행렬 UU 의 놈을 ∥U∥\left\| U \right\| 와 같이 나타내자.
알고리즘: 그램-슈미트 직교화 GS(A)GS(A)
In행렬 A∈Rn×nA \in \mathbb{R}^{n \times n}
1.U,V←A,AU, V \gets A, A
2.U:1←Vi1∥Vi1∥U_{:1} \gets {\frac{ V_{i1} }{ \left\| V_{i1} \right\| }}
3.for j=2,⋯ ,nj = 2, \cdots, n do
4.　　 for jp=1,⋯ ,(j−1)j_{p} = 1, \cdots, (j-1) do
5.　　　　 V:j←V:j−(U:jTU:jp)U:jpV_{:j} \gets V_{:j} - \left( U_{:j}^{T} U_{:j_{p}} \right) U_{:j_{p}}
6.　　 end for
7.　　 U:j←V:j∥V:j∥U_{:j} \gets {\frac{ V_{:j} }{ \left\| V_{:j} \right\| }}
8.end for
Out정규직교행렬 UU 와 직교행렬 VV 를 얻는다.

우선은 선형대수에서 널려진 그램-슈미트 직교화와 달리 위과 같이 축의 길이가 정규화된 UU 와 정규화되지 않고 직교하기만 한 행렬 VV 를 얻는 함수 GSGS 를 정의한다. 여기서 VV 의 각 열은 변분 방정식 Y˙=JY\dot{Y} = J Y 에서 퍼터베이션 YY 가 원래의 방정식 x˙=f(x)\dot{x} = f(x) 의 특정한 점 x(t)x(t) 에서 선형변환 J(x(t))J \left( x(t) \right) 에 의해 축의 길이가 얼마나 늘어나고 줄어드는지를 의미한다.

알고리즘: 랴푸노프 스펙트럼 계산
In트래젝터리 {xk}k=0K\left\{ x_{k} \right\}_{k=0}^{K}
자코비안 행렬 J(x)J(x)
1.λ←0\lambda \gets \mathbf{0}# 길이 nn 인 영벡터
2.U←InU \gets I_{n}# 초기화
3.for k=1,⋯ ,Kk = 1, \cdots, K do
4.　　 U,V←GS(U)U, V \gets GS(U)
5.　　 λ←λ+(log⁡∥V:1∥,⋯ ,log⁡∥V:n∥)\lambda \gets \lambda + \left( \log \left\| V_{:1} \right\| , \cdots , \log \left\| V_{:n} \right\| \right)
6.　　 U←RK4(J(xk),U,h)U \gets RK4 \left( J \left( x_{k} \right) , U , h \right)
7.end for
8.λ←λ/(hK)\lambda \gets \lambda / (hK)
Out랴푸노프 스펙트럼 λ\lambda

QR 분해 기반 2

행렬 AA 에 대해 QR(A)QR(A) 이 AA 의 QR 분해 결과인 QQ 와 RR 을 리턴하는 함수라 하자.

알고리즘: 랴푸노프 스펙트럼 계산
In트래젝터리 {xk}k=0K\left\{ x_{k} \right\}_{k=0}^{K}
자코비안 행렬 J(x)J(x)
1.λ←0\lambda \gets \mathbf{0}# 길이 nn 인 영벡터
2.J←J(x0)J \gets J\left( x_{0} \right)
2.Q,R←QR(J)Q, R \gets QR(J)# 초기화
3.for k=1,⋯ ,Kk = 1, \cdots, K do
4.　　 Q,R←QR(J)Q, R \gets QR(J)
5.　　 λ←λ+(log⁡∣R11∣,⋯ ,log⁡∣Rnn∣)\lambda \gets \lambda + \left( \log \left| R_{11} \right| , \cdots , \log \left| R_{nn} \right| \right)# RR 의 대각성분
6.　　 Q←RK4(J(xk),Q,h)Q \gets RK4 \left( J \left( x_{k} \right) , Q , h \right)
7.end for
8.λ←λ/(hK)\lambda \gets \lambda / (hK)
Out랴푸노프 스펙트럼 λ\lambda

설명

ODE로 표현되는 동역학계의 랴푸노프 스펙트럼을 계산하는 것은 바이퍼케이션 다이어그램을 계산하는 것에 비해 훨씬 어려운 일이다. 카오스라고 하는 개념 자체가 랴푸노프 스펙트럼으로 정의되는만큼 다이내믹스를 공부하는 입장에서는 대단히 중요한 기술이며, 계산량에 있어서 커다란 벽을 느끼게 되는 파트기도 하다. 당연하지만, 랴푸노프 지수에 대해 어렴풋이 알고 있는 것과 실제로 계산할 줄 아는 것에는 큰 차이가 있다.

소개된 두 가지 방식은 일견 아무 상관도 없어보이지만, 사실 QR 분해의 존재성 증명을 자세히 살펴보면 본질적으로 매우 유사하다. 차이점에 주목해보자면 다음과 같이 요약할 수 있다.

  • GS 기반: 기하적으로 직관적이라서 이해하기 쉽다. 알고리즘을 잘 읽어보면 유닛 볼이 가진 nn개의 축을 열벡터로 가지는 행렬 YY 를 JJ 로 변형시킨 일립소이드 축의 길이를 재는 것 그 자체다.
  • QR 기반: 구현에 있어서 비교적 손이 덜 간다. QR 분해 후 RR 의 대각성분만 잘 모으면 되고, 애초에 QR 분해부터 패키지를 사용할 가능성이 높아 프로그램의 신뢰성도 높다.

결과 비교

ρ∈[0,100]\rho \in [0, 100] 에서 로렌츠 어트랙터의 랴푸노프 스펙트럼을 두가지 방법으로 계산한 결과는 다음과 같다.

  • 그램-슈미츠 기반

alt text

  • QR 분해 기반

alt text

  • 레퍼런스 3

alt text

보다시피 그림으로 그려서는 구분할 수 없을 정도로 비슷한 결과를 얻었고, 레퍼런스와 비교해도 신뢰할만한 것을 확인했으니 구체적인 데이터로써 비교해보자. 파란색(왼쪽)은 그램-슈미츠 기반, 주황색(오른쪽)은 QR 분해 기반으로 계산한 랴푸노프 스펙트럼이다.

alt text

성능 비교

alt text

  • 그램-슈미츠 기반(위): 2:39:31
  • QR 분해 기반(아래): 9:50:20

그런데 두 기법은 속도에서 상당한 차이가 난다. 어떤 방법을 사용하든 결과 자체는 낼 수 있지만, 위 수치들을 계산하기 위한 벤치마크를 보면 3배에서 4배 정도의 속도 차이를 확인할 수 있다.

다음은 보다 정확한 비교를 위해 로렌츠 어트랙터를 구하는 시간을 측정한 것이다. 약 9분이 소요되었으니 랴푸노프 스펙트럼을 계산하는 것 자체에는 큰 영향을 미치지 못한 것을 확인할 수 있다.

alt text

코드 4 5

전체코드

로렌츠 어트랙터의 ρ∈[0,100]\rho \in [0, 100] 에서 랴푸노프 스펙트럼을 계산하는 줄리아 코드다.

using CSV, DataFrames, LinearAlgebra, ProgressMeter, Plots, Dates; @info now()
using Base.Threads: @threads # Base.Threads.nthreads()

function RK4(f::Function, v::AbstractVector, h=1e-3)
    V1 = f(v)
    V2 = f(v + (h/2)*V1)
    V3 = f(v + (h/2)*V2)
    V4 = f(v + h*V3)
    return v + (h/6)*(V1 + 2V2 + 2V3 + V4), V1
end
function RK4(J::AbstractMatrix, U::AbstractMatrix, dt=1e-3)
    V1 = J*U
    V2 = J*(U + (dt/2)*V1)
    V3 = J*(U + (dt/2)*V2)
    V4 = J*(U + dt*V3)
    return U + (dt/6)*(V1 + 2V2 + 2V3 + V4)
end
function factory_lorenz(ρ::Number; ic = [10.,10.,10.], tspan = [0., 1000.], dt = 1e-3)
    σ = 10
    β = 8/3
    function sys(v::AbstractVector)
        x, y, z = v
        dx = σ*(y - x)
        dy = x*(ρ - z) - y
        dz = x*y - β*z
        return [dx, dy, dz]
    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_
        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_lorenz(T::Type, args...; kargs...) =
DataFrame(factory_lorenz(args...; kargs...), ["x", "y", "z", "dx", "dy", "dz"])

function gram_schmidt(J)
    N = size(J, 1)
    U, V = deepcopy(J), deepcopy(J)
    U[:,1] = V[:,1] / norm(V[:,1])
    for j in 2:N
        for jp in 1:(j-1)
            V[:,j] -= (J[:,j]'U[:,jp])*U[:,jp]
        end
        U[:,j] = V[:,j] / norm(V[:,j])
    end
    return U, V
end
J_lorenz(x,y,z,ρ) = [
     -10   10  0
      ρ-z -1  -x
        y  x  -8/3
]

const dt = 1e-3
ρ_range = 0:.1:100
cd("G:/DDM")

function lyapunov_lorenz_GS()
    schedules = DataFrame(idx = eachindex(ρ_range), ρ = ρ_range)
    schedules[!, :λ1] .= .0; schedules[!, :λ2] .= .0; schedules[!, :λ3] .= .0; 
    @showprogress @threads for dr = eachrow(schedules)
        data = factory_lorenz(DataFrame, dr.ρ; dt)
        
        λ = zeros(3);
        J = J_lorenz(collect(data[1, 1:3])..., dr.ρ)
        U = I(3)
        for i in 2:nrow(data)
            U, V = gram_schmidt(U)
            λ += V |> eachcol .|> norm .|> log          

            U = RK4(J, U, dt)
            J = J_lorenz(collect(data[i, 1:3])..., dr.ρ)
        end

        λ ./= dt*nrow(data)
        dr[[:λ1, :λ2, :λ3]] .= sort(λ, rev = true)
    end
    sort!(schedules, :ρ)
    CSV.write("lyapunov/lorenzGS.csv", schedules, bom = true)
    plot(xticks = 0:20:100, legend = :none, xlabel = "ρ", ylabel = "λ", size = [600, 300])
    plot!(schedules.ρ, schedules.λ1, lw = 2, color = 1)
    plot!(schedules.ρ, schedules.λ2, lw = 2, color = 2)
    plot!(schedules.ρ, schedules.λ3, lw = 2, color = 3)
    png("lyapunov/lorenzGS.png")
end

function lyapunov_lorenz_QR()
    schedules = DataFrame(idx = eachindex(ρ_range), ρ = ρ_range)
    schedules[!, :λ1] .= .0; schedules[!, :λ2] .= .0; schedules[!, :λ3] .= .0; 
    @showprogress @threads for dr = eachrow(schedules)
        data = factory_lorenz(DataFrame, dr.ρ; dt)
        
        λ = zeros(3);
        J = J_lorenz(collect(data[1, 1:3])..., dr.ρ)
        Q, _ = qr(J)
        for i in 2:nrow(data)
            Q, R = qr(Q)
            λ += R .|> abs |> diag .|> log   

            Q = RK4(J, Matrix(Q), dt)
            J = J_lorenz(collect(data[i, 1:3])..., dr.ρ)
        end

        λ ./= dt*nrow(data)
        dr[[:λ1, :λ2, :λ3]] .= sort(λ, rev = true)
    end
    sort!(schedules, :ρ)
    CSV.write("lyapunov/lorenzQR.csv", schedules, bom = true)
    plot(xticks = 0:20:100, legend = :none, xlabel = "ρ", ylabel = "λ", size = [600, 300])
    plot!(schedules.ρ, schedules.λ1, lw = 2, color = 1)
    plot!(schedules.ρ, schedules.λ2, lw = 2, color = 2)
    plot!(schedules.ρ, schedules.λ3, lw = 2, color = 3)
    png("lyapunov/lorenzQR.png")
end

lyapunov_lorenz_GS()
lyapunov_lorenz_QR()

같이보기

  1. Parker, Thomas S., and Leon Chua. Practical numerical algorithms for chaotic systems. Springer Science & Business Media, 2012. ↩︎

  2. von Bremen, Hubertus F., Firdaus E. Udwadia, and Wlodek Proskurowski. “An efficient QR based method for the computation of Lyapunov exponents.” Physica D: Nonlinear Phenomena 101.1-2 (1997): 1-16. https://doi.org/10.1016/S0167-2789(96)00216-3 ↩︎

  3. Balcerzak, M., Pikunov, D. & Dabrowski, A. The fastest, simplified method of Lyapunov exponents spectrum estimation for continuous-time dynamical systems. Nonlinear Dyn 94, 3053–3065 (2018). https://doi.org/10.1007/s11071-018-4544-z ↩︎

  4. Lutz Lehmann, Numerical computation of Lyapunov exponent, URL (version: 2023-12-07): https://scicomp.stackexchange.com/q/36023 ↩︎

  5. https://github.com/JuliaDynamics/ChaosTools.jl/blob/4053d184aa91b62b550ecd722193225ea2b198cd/src/chaosdetection/lyapunovs/lyapunovspectrum.jl#L5-L51 ↩︎