📂動力学

線形システムのリアプノフスペクトル

定理

行列 $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")