ローレンツ・アトラクター 📂動力学

ローレンツ・アトラクター

概要

ローレンツ方程式^{lorenz equation}とは、大気の対流を連立常微分方程式によって表現する数学的モデルである。

システム

$A\_Trajectory\_Through\_Phase\_Space\_in\_a\_Lorenz\_Attractor.gif$

$$ \begin{align*} {{dx} \over {dt}} =& - \sigma x + \sigma y \\ {{dy} \over {dt}} =& - xz + \rho x - y \\ {{dz} \over {dt}} =& xy - \beta z \end{align*} $$

変数

$x(t)$: $t$ 時点での粒子の $x$ 座標を示す。
$y(t)$: $t$ 時点での粒子の $y$ 座標を示す。
$z(t)$: $t$ 時点での粒子の $z$ 座標を示す。

パラメータ

$\sigma$: 粘性と熱伝導率に関するパラメータである プラントル数^{prandtl number}。
$\rho$: 流体の熱伝達方法に関するパラメータである レイリー数^{rayleigh number}。
$\beta$: ウィキペディア¹によれば次元に関連するパラメータらしい。教科書でもこれが何かを正確に説明していない場合が多い。

初めてシミュレーションを行う際、これらのパラメータは慣例的にローレンツに従って $\sigma = 10$, $\rho = 28$, $\displaystyle \beta = {{8} \over {3}}$ と設定する場合が多い。レイリー数が $\rho \approx 24.74$ を超えると、ローレンツアトラクタが描く粒子の動きはカオス^chaoticになる。非常にわずかな初期条件の変化にも敏感に反応し、それを予測することは極めて難しい。蝶の羽ばたきが嵐になるかもしれないのだ。

説明

ビフルケーション

$20190207\_143724.png$ alt text

上の図はローレンツアトラクタのビフルケーションダイアグラムである²。横軸は $25$ から $325$ まで少しずつ変わる $\rho$ の値、縦軸は十分に長い時間後の軌道での $z$ 方向の極大値を示している。$r= \rho$ が小さい時は非常に不安定な動きを見せるが、ある一定以上に大きくなると安定した軌道を持つことが分かる。以下の説明を読むことをお勧めする：

$20190207\_142108.png$

上のグラフは $\gamma=100$ の時の軌跡を描いたものである。通常$\gamma$がある程度以上大きいと、初めはカオスな動きをするように見えるが、次第に周期的になることがある。

$20190207\_142944.png$

上のグラフは初期の軌跡を削除し、長時間後の粒子の軌跡だけを表している。完全に周期的ではないが、ほぼ予測可能な行動を示すようになった。

$20190207\_160923.png$

軌跡を見やすくするために $xz$ グラフを描くと上記のような形になる。確かに最初はカオスに見えたが、時間が経つにつれて軌道が安定し、極大値は3つ程度で明瞭に現れる。実際にビフルケーションダイアグラムを見ても同じような解釈ができる。

一方、$\gamma = 28$ の場合は十分に長い時間が経過しても次のようにカオスな動きを示す。ビフルケーションダイアグラムでは多くの点が散らばり、太い形状を形成する。

$20190207\_160633.png$

反対に、$\gamma = 300$ の場合は次のように明瞭に現れる2つの極大値だけが存在する。ビフルケーションダイアグラムでは、1つの $\gamma$ 値に対応する2つの点が現れることで2つの線が示される。

$20190207\_161132.png$

リアプノフスペクトル

以下は $\rho \in [0, 100]$ におけるリアプノフスペクトルを示す³。

alt text

コード

マトラブ

以下はマトラブを使ってローレンツアトラクタの軌跡を描画するコードである。

sigma = 10;
beta = 8/3;
rho = 28;
f = @(t,v) [-sigma*v(1) + sigma*v(2); rho*v(1) - v(2) - v(1)*v(3); -beta*v(3) + v(1)*v(2)];
[t,v] = ode45(f,[0 10000],[1 1 1]);     % Runge-Kutta 4th/5th order ODE solver
x=v(:,1); x=x(50000:end);
y=v(:,2); y=y(50000:end);
z=v(:,3); z=z(50000:end);
plot3(x,y,z)
%plot(x,z)

フォートラン

以下はフォートランを使ってビフルケーションダイアグラムのソースを作成し、OUTPUT.csvファイルに出力するコードである。ODEを解くために使用されたメソッドはRK4である。

program LorenzBifurcation
  implicit none
  double precision :: sigma=10.d0, rho=28.d0, beta=8.d0/3.d0
  Integer, parameter :: D = 3
  double precision :: v(D)=1.d0, h=0.001d0, z1, z2, z3
  integer :: i, j, k, preitr=10**4-1000, endstep=10**4
  character :: answer
  abstract interface
     function func (v)
       Integer, parameter :: D = 3
       double precision :: v(D), func(D)
     end function func
  end interface
  open(UNIT=34,FILE="OUTPUT.csv",ACTION="WRITE")
  rho=5.d0
  h=0.01d0
  do while(rho<325)
    print '(A,f6.4,A)', "Now", (rho-5.d0)/320.d0, "%"
    do i=1, preitr
       v = RK4(v,f)
    end do
    z1 = preitr
    z2 = preitr
    z3 = preitr
    do i=preitr+1, endstep
       v = RK4(v,f)
       z1 = z2
       z2 = z3
       z3 = v(3)
       if(z1<z2 .and. z2>z3) then
         write (34,*) rho, ",", z2
       end if
    end do
    v=1.d0
    rho=rho+0.1d0
  end do
contains
  function RK4(v,f)
    procedure (func) :: f
    double precision :: v(D), V1(D), V2(D), V3(D), V4(D), RK4(D)
    V1 = f(v)
    V2 = f(v + (h/2)*V1)
    V3 = f(v + (h/2)*V2)
    V4 = f(v + h*V3)
    RK4 = v + (h/6)*(V1 + 2*V2 + 2*V3 + V4)
  end function RK4
  
  function f(v)
    double precision :: v(D), f(D)
    f(1) = sigma*(v(2)-v(1))
    f(2) = v(1)*(rho-v(3))-v(2)
    f(3) = v(1)*v(2) - beta*v(3)
  end function f
end program LorenzBifurcation

以下はマトラブを使って上記で得たOUTPUT.csvファイルを使用し、ビフルケーションダイアグラムを描画するコードである。

[file,path] = uigetfile('*.csv');
data = csvread([path,file]);
 
sz = 0.6;
scatter(data(:,1),data(:,2),sz,'.');

ジュリア

以下はジュリアを使ってローレンツアトラクタを見つけ、図を描くコードである。

using DifferentialEquations
using Plots

function parameterized_lorenz!(du,u,p,t)
    x,y,z = u
    σ,ρ,β = p
    du[1] = dx = σ*(y-x)
    du[2] = dy = x*(ρ-z) - y
    du[3] = dz = x*y - β*z
end
u0 = [1.0,0.0,0.0]
tspan = (0.0,100.0)
p = [10.0,28.0,8/3]
prob = ODEProblem(parameterized_lorenz!,u0,tspan,p)
solution = solve(prob)
plot(solution, vars = (1,2,3))

以下はビフルケーションダイアグラムを描くコードである。

using DataFrames, CSV, Plots

function RK4(f::Function, v::AbstractVector, h=10^(-2))
    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 factory_lorenz(idx::Int64, ρ::Number; ic = [10.,10.,10.], tspan = [0., 10.])
    σ = 10
    β = 8/3
    function lorenz(v::AbstractVector)
        x, y, z = v
        dx = σ*(y - x)
        dy = x*(ρ - z) - y
        dz = x*y - β*z
        return [dx, dy, dz]
    end

    dt = 10^(-3)
    t_ = first(tspan):dt:last(tspan)
    
    ndatapoints = count(first(tspan) .< t_ .≤ last(tspan))
    len_t_ = length(t_)

    v = ic; DIM = length(v)
    traj = zeros(2DIM, len_t_+1)
    traj[1:DIM, 1] = v

    for tk in eachindex(t_)
        v, dv = RK4(lorenz, v, dt)
        if tk+1 ≥ (len_t_ - ndatapoints)
            traj[        1:DIM , tk+1] =  v
            traj[DIM .+ (1:DIM), tk  ] = dv
        end
    end
    traj = traj[:, 1:(end-1)]'
    traj = traj[(end-ndatapoints):end, :]

    return traj
end
factory_lorenz(T::Type, args...; ic = [10,10,10.], tspan = [0., 100.]) =
DataFrame(factory_lorenz(args...; ic = ic, tspan), ["x", "y", "z", "dx", "dy", "dz"])

r_range = 25:0.1:325
schedules = DataFrame(idx = eachindex(r_range), r = r_range)
CSV.write("G:/DDM/lyapunov/lorenz_schedules.csv", schedules)
# schedules = CSV.read("G:/DDM/lyapunov/lorenz_schedules.csv", DataFrame)
@showprogress @threads for dr in eachrow(schedules)
    data = factory_lorenz(DataFrame, dr.idx, dr.r, tspan = [0, 100])[90000:end, :]
    CSV.write("G:/DDM/lyapunov/lorenz/$(lpad(dr.idx, 5, '0')).csv", data)
end

idcs = Int64[]; vrtc = Float64[]; hrzn = Float64[]
@showprogress for dr = eachrow(schedules)
    filename = "G:/DDM/lyapunov/lorenz/$(lpad(dr.idx, 5, '0')).csv"
    !isfile(filename) && continue

    data = CSV.read(filename, DataFrame)
    idx_sampled = findall((circshift(data.z, 1) .< data.z) .&& (circshift(data.z, -1) .< data.z))[2:(end-1)]
    sampledz = data.z[idx_sampled]
    append!(idcs, fill(dr.idx, length(sampledz)))
    append!(hrzn, fill(dr.r, length(sampledz)))
    append!(vrtc, sampledz)
end
CSV.write("G:/DDM/lyapunov/lorenz_bifurcation.csv", DataFrame(; vrtc, hrzn))
scatter(hrzn, vrtc, color = :black, ms = 1, legend = :none, msw = 0, ma = 0.1, xticks = [25, 325], size = [800, 800], xlabel = L"\rho");
png("G:/DDM/lyapunov/lorenz_bifurcation.png")

以下はリアプノフ指数を計算するコードである。

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

lyapunov_lorenz_GS()

https://en.wikipedia.org/wiki/Lorenz_system#Overview ↩︎
Yorke. (1996). CHAOS: An Introduction to Dynamical Systems: p368. ↩︎
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 ↩︎