📂동역학

선형 시스템의 랴푸노프 스펙트럼

정리

행렬 $A \in \mathbb{R}^{n \times n}$ 에 대해 다음과 같은 벡터필드가 선형 미분방정식으로 주어져 있다고 하자. $$\dot{\mathbf{x}} = A \mathbf{x} \qquad , \mathbf{x} \in \mathbb{R}^{n}$$

만약 $A$ 가 대칭행렬 이면, 이 시스템의 랴푸노프 스펙트럼은 ${\frac{ 1 }{ 2 }} \left( A + A^{T} \right)$ 의 고유값과 같다.

---

• 여기서 $A^{T}$ 는 행렬 $A$ 의 전치행렬이다.

설명

$$A = \begin{bmatrix} -0.2 & -0.3 \\ -0.2 & -0.7 \end{bmatrix}$$

julia> A = [-.2 -.3; -.3 -.7]
2×2 Matrix{Float64}:
 -0.2  -0.3
 -0.3  -0.7

julia> @showprogress for t in 0:h:tend
           # v, _ = RK4(T_A, v_[end])
           # push!(v_, v)

           U, V = gram_schmidt(U)
           λ += V |> eachcol .|> norm .|> log
           isinteger(t) && push!(λ_, λ / t)
           U = RK4(A, U, h)
       end
Progress: 100%|█████████████████████| Time: 0:00:05

구체적인 예로써 위와 같은 행렬 $A = A^{T}$ 로 주어진 선형 시스템에서 랴푸노프 스펙트럼을 계산해보자.

julia> λ / (tend)
2-element Vector{Float64}:
 -0.05958668545713238
 -0.8404133145577668

julia> eigen((A + A')/2).values
2-element Vector{Float64}:
 -0.8405124837953326
 -0.0594875162046673

수치 계산 결과를 비교해보면 이론적인 예측치와 상당히 유사하게 계산되는 것을 확인할 수 있다.

이렇게 쉬운 시스템의 랴푸노프 스펙트럼을 굳이 계산하는 의미는 코드의 검증에 있다. 보통의 복잡한 시스템과 달리 선형 시스템은 이론적으로 랴푸노프 스펙트럼을 계산할 수 있어서, RK4 및 그램 슈미트 직교화 혹은 QR 분해와 그 순서가 정확한지 확인할 수 있는 예시가 되는 것이다.

증명 ¹

랴푸노프 스펙트럼의 정의:

변분 방정식: $f$ 의 자코비안 행렬 $J$ 에 대해 다음을 변분 방정식^{variational equation}이라 한다. $$\dot{Y} = J Y$$ 여기서 행렬함수 $Y = Y(t) \in \mathbb{R}^{n \times n}$ 의 초기 조건은 항등행렬 $Y(0) = I$ 로 둔다. … 기하적으로, $Y$ 는 원래 시스템의 $x(0)$ 에서 조금 움직인 $x(t)$ 로 변하는동안 그 탄젠트 벡터 자체가 어떻게 작용하는지를 보여준다고 생각할 수 있다.

$$\lambda_{k} := \lim_{t \to \infty} \log \left[ \left( {\frac{ \left\| Y(t, v) v \right\| }{ \left\| v \right\| }} \right)^{1/t} \right]$$ 위와 같이 정의된 $\left\{ \lambda_{1} , \cdots , \lambda_{n} \right\}$ 를 랴푸노프 스펙트럼^{Lyapunov spectrum}이라 하거나, $$\Lambda_{v} := \lim_{t \to \infty} \left[ Y(t)^{\ast} Y(t) \right]^{1/2t}$$ 위와 같이 정의된 행렬 $\Lambda_{v}$ 의 고유값 $\mu_{1} , \cdots , \mu_{n}$ 에 로그를 취한 $\lambda_{k} := \log \mu_{k}$ 를 랴푸노프 스펙트럼이라 한다.

$\dot{\mathbf{x}} = A \mathbf{x}$ 은 선형 미분방정식이므로 그 자코비안 $J$ 역시 $J = A$ 가 되어 원래 시스템의 변분 방정식은 항등행렬인 초기값 $Y = I_{n}$ 에 대해 $\dot{Y} = A Y$ 이고 그 정확한 해는 $Y(t) = e^{A t}$ 가 된다. 랴푸노프 스펙트럼의 정의에 따라 다음을 얻는다. $$\begin{align*} & \log \lim_{t \to \infty} \left[ Y(t)^{T} Y(t) \right]^{1/2t} \\ =& \lim_{t \to \infty} \log \left[ Y(t)^{T} Y(t) \right]^{1/2t} \\ =& \lim_{t \to \infty} {\frac{ 1 }{ 2t }} \log \left[ Y(t)^{T} Y(t) \right] \\ =& \lim_{t \to \infty} {\frac{ 1 }{ 2t }} \log \left[ e^{A^{T} t} e^{A t} \right] \end{align*}$$

행렬의 지수와 로그: 만약 두 행렬 $A, B \in \mathbb{C}^{n \times n}$ 가 $AB = BA$ 를 만족하면, 다음이 성립한다. $$\log AB = \log A + \log B - 2 \pi i \mathcal{U} \left( \log A + \log B \right)$$ 특히 $A$ 의 고유값 $\mu_{k}$ 와 $B$ 의 고유값 $\nu_{k}$ 에 대해 다음의 따름정리를 얻을 수 있다. $$\log AB = \log A + \log B \iff \forall k = 1 , \cdots , n : \arg \mu_{k} + \arg \nu_{k} \in ( - \pi , \pi ]$$

$A \in \mathbb{R}^{n \times n}$ 는 실수 엔트리만을 가지므로 에르미트 행렬이고, 그 고유값은 실수이므로 $e^{A^{T} t}$ 과 $e^{A t}$ 의 모든 고유값도 실수다. 고유값이 실수라는 것은 이들의 편각이 $0$ 이라는 것이므로 행렬 로그 $\log A B = \log A + \log B$ 가 성리해서 다음을 얻는다. $$\begin{align*} & \lim_{t \to \infty} {\frac{ 1 }{ 2t }} \log \left[ e^{A^{T} t} e^{A t} \right] \\ =& \lim_{t \to \infty} {\frac{ 1 }{ 2t }} \left[ \log e^{A^{T} t} + \log e^{A t} \right] \\ =& \lim_{t \to \infty} {\frac{ 1 }{ 2t }} \left[ A^{T} t + A t \right] \\ =& {\frac{ 1 }{ 2 }} \left[ A^{T} + A \right] \end{align*}$$

■

코드

다음은 줄리아를 통해 실제로 랴푸노프 스펙트럼을 계산한 결과다.

using ProgressMeter, LinearAlgebra

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
function RK4(J::AbstractMatrix, U::AbstractMatrix, dt=1e-2)
    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

A = [-.2 -.3; -.3 -.7]
# T_A(v) = A*v

U = I(2)
λ = zeros(2); λ_ = [λ]
# v = [10, 10.0]; v_ = [v]
h = 1e-3; tend = 1000
@showprogress for t in 0:h:tend
    # v, _ = RK4(T_A, v_[end])
    # push!(v_, v)

    U, V = gram_schmidt(U)
    λ += V |> eachcol .|> norm .|> log
    isinteger(t) && push!(λ_, λ / t)
    U = RK4(A, U, h)
end
λ / (tend)
eigen((A + A')/2).values

using Plots
plot(xlims = [1, 100], xlabel = "t", legend = :right)
scatter!(0:10:1001, first.(λ_[1:10:end]), label = "λ1", msw = 0)
scatter!(0:10:1001,  last.(λ_[1:10:end]), label = "λ2", msw = 0)
hline!(eigen((A + A')/2).values, color = :black, label = "Theoretical")