微分方程式で表されるシステムのリアプノフスペクトルとその数値計算法
アルゴリズム
空間 $X = \mathbb{R}^{n}$ と 関数 $f : X \to X$ に対して次のような ベクトル場 が 微分方程式 として与えられているとしよう。 $$ \dot{x} = f(x) $$ $f$ の ヤコビ行列 $J$ に対して次のように 変分方程式 を定義しよう。 $$ \dot{Y} = J Y $$ ここで 行列関数 $Y = Y(t) \in \mathbb{R}^{n \times n}$ の初期条件は 単位行列 $Y(0) = I_{n}$ にする。システム $\dot{x} = f(x)$ の リアプノフスペクトル $\lambda_{1} , \cdots , \lambda_{n}$ は $\dot{x} = f(x)$ の数値解（例えば RK4 などで計算されたもの） $\left\{ x_{k} \right\}_{k=0}^{K}$ に対して変分方程式を繰り返し解くことで近似的に計算される。$J(x)$ は元のシステムの トラジェクトリー $x(t)$ に依存する。
- ここで変分方程式 $\dot{Y} = J Y$ を解くとは、パートベーションperterbationと呼ばれる 行列 または ベクトル $Y$ を数値解を得る際に用いた RK4 などの手法で解くことを意味する。$x_{k}$ を得るために使用した十分に小さい時間間隔が $h \approx 0$ であるならば変分方程式も同じ時間間隔 $h$ を使う必要がある。
- 行列 $J$ と $U$ に対して ルンゲクッタ法 $RK4(J, U, h)$ は 次のように定義される。 $$ \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
- 行列 $U$ の $j$番目の列ベクトルを $U_{:j}$ のように表そう。
- 行列 $U$ の 転置行列を $U^{T}$ のように表そう。
- 行列 $U$ の ノルムを $\left\| U \right\|$ のように表そう。
| アルゴリズム: グラム-シュミット正規直交化 $GS(A)$ | ||
|---|---|---|
| In | 行列 $A \in \mathbb{R}^{n \times n}$ | |
| 1. | $U, V \gets A, A$ | |
| 2. | $U_{:1} \gets {\frac{ V_{i1} }{ \left\| V_{i1} \right\| }} $ | |
| 3. | for $j = 2, \cdots, n$ do | |
| 4. | for $j_{p} = 1, \cdots, (j-1)$ do | |
| 5. | $V_{:j} \gets V_{:j} - \left( U_{:j}^{T} U_{:j_{p}} \right) U_{:j_{p}} $ | |
| 6. | end for | |
| 7. | $U_{:j} \gets {\frac{ V_{:j} }{ \left\| V_{:j} \right\| }}$ | |
| 8. | end for | |
| Out | 正規直交行列 $U$ と直交行列 $V$ を得る。 |
まずは線形代数でよく知られたグラム-シュミット正規直交化とは異なり、このように軸の長さが正規化された $U$ と正規化されていないが直交だけしている行列 $V$ を得る関数 $GS$ を定義する。ここで $V$ の各列は変分方程式 $\dot{Y} = J Y$ でのパートベーション $Y$ が元の方程式 $\dot{x} = f(x)$ の特定の点 $x(t)$ での線形変換 $J \left( x(t) \right)$ により軸の長さがどれだけ長くなり、またどれだけ短くなるかを意味する。
| アルゴリズム: リアプノフスペクトル計算 | ||
|---|---|---|
| In | トラジェクトリー $\left\{ x_{k} \right\}_{k=0}^{K}$ ヤコビ行列 $J(x)$ | |
| 1. | $\lambda \gets \mathbf{0}$ | # 長さ $n$ のゼロベクトル |
| 2. | $U \gets I_{n}$ | # 初期化 |
| 3. | for $k = 1, \cdots, K$ do | |
| 4. | $U, V \gets GS(U)$ | |
| 5. | $\lambda \gets \lambda + \left( \log \left\| V_{:1} \right\| , \cdots , \log \left\| V_{:n} \right\| \right)$ | |
| 6. | $U \gets RK4 \left( J \left( x_{k} \right) , U , h \right)$ | |
| 7. | end for | |
| 8. | $\lambda \gets \lambda / (hK)$ | |
| Out | リアプノフスペクトル $\lambda$ |
QR分解に基づく 2
行列 $A$ に関して $QR(A)$ が $A$ の QR分解結果である $Q$ および $R$ を返す関数であるとしよう。
| アルゴリズム: リアプノフスペクトル計算 | ||
|---|---|---|
| In | トラジェクトリー $\left\{ x_{k} \right\}_{k=0}^{K}$ ヤコビ行列 $J(x)$ | |
| 1. | $\lambda \gets \mathbf{0}$ | # 長さ $n$ のゼロベクトル |
| 2. | $J \gets J\left( x_{0} \right)$ | |
| 2. | $Q, R \gets QR(J)$ | # 初期化 |
| 3. | for $k = 1, \cdots, K$ do | |
| 4. | $Q, R \gets QR(J)$ | |
| 5. | $\lambda \gets \lambda + \left( \log \left| R_{11} \right| , \cdots , \log \left| R_{nn} \right| \right)$ | # $R$ の対角成分 |
| 6. | $Q \gets RK4 \left( J \left( x_{k} \right) , Q , h \right)$ | |
| 7. | end for | |
| 8. | $\lambda \gets \lambda / (hK)$ | |
| Out | リアプノフスペクトル $\lambda$ |
説明
ODEで表される動的システムのリアプノフスペクトルを計算するのは、分岐図を計算することに比べてはるかに難しいことである。カオスという概念自体がリアプノフスペクトルで定義されるため、ダイナミクスを勉強する観点からは非常に重要な技術であり、計算量において大きな壁を感じるパートでもある。当たり前だが、リアプノフ指数についてぼんやりと知っていることと、実際に計算できることには大きな違いがある。
紹介された二つの方法は一見何の関連もないように見えるが、実はQR分解の存在証明を詳しく調べてみると本質的に非常に類似している。違いに注目すると次のように要約できる。
- GSに基づく: 幾何的に直感的で理解しやすい。アルゴリズムをよく読むとユニットボールが持つ $n$個の軸を列ベクトルとする行列 $Y$ を $J$ で変形した 楕円体の軸の長さを測るそのものだ。
- QRに基づく: 実装において比較的手間が少ない。QR分解後 $R$ の対角成分だけをしっかりと集めればよく、そもそもQR分解からパッケージを使用する可能性が高いため、プログラムの信頼性も高い。
結果の比較
$\rho \in [0, 100]$ で ローレンツアトラクタのリアプノフスペクトルを二つの方法で計算した結果は次のようになる。
- グラム-シュミットに基づく
- QR分解に基づく
- リファレンス 3
ご覧の通り、図としては区別がつかないほど似た結果を得ており、リファレンスと比較しても信頼に足るものを確認したので、具体的なデータとして比較してみよう。青色（左）はグラム-シュミットに基づくもので、オレンジ色（右）はQR分解に基づくもので計算したリアプノフスペクトルだ。
性能の比較
- グラム-シュミットに基づく（上）: 2:39:31
- QR分解に基づく（下）: 9:50:20
ところが、二つの手法には速度でかなりの差がある。どんな方法を使おうとも結果自体は出すことができるが、上記の数値を計算するためのベンチマークを見ると、3倍から4倍程度の速度差を確認できる。
次はより正確な比較のためのローレンツアトラクタを求める時間を測定したものである。約9分がかかったので、リアプノフスペクトルを計算すること自体には大きな影響を及ぼさなかったことが確認できる。
コード 4 5
全体コード
ローレンツアトラクタの $\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()
関連情報を見る
Parker, Thomas S., and Leon Chua. Practical numerical algorithms for chaotic systems. Springer Science & Business Media, 2012. ↩︎
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 ↩︎
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 ↩︎
Lutz Lehmann, Numerical computation of Lyapunov exponent, URL (version: 2023-12-07): https://scicomp.stackexchange.com/q/36023 ↩︎
https://github.com/JuliaDynamics/ChaosTools.jl/blob/4053d184aa91b62b550ecd722193225ea2b198cd/src/chaosdetection/lyapunovs/lyapunovspectrum.jl#L5-L51 ↩︎