로렌츠 어트랙터
개요
로렌츠 방정식lorenz equation이란 대기의 대류를 연립 상미분 방정식으로써 표현하는 수학적 모델이다.
시스템
$$ \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$: 위키피디아1에 따르면 차원에 관련된 파라미터인 것 같다. 교재에서도 이것이 무엇인지 정확하게 설명하지 않는 경우가 많다.
처음 시뮬레이션을 할 때 이 파라미터들은 관례적으로 로렌츠가 했던대로 $\sigma = 10$, $\rho = 28$, $\displaystyle \beta = {{8} \over {3}}$ 와 같이 두는 경우가 많다. 레일리 수가 $\rho \approx 24.74$ 을 넘어가면 로렌츠 어트랙터가 묘사하는 입자의 움직임은 캐어릭chaotic해진다. 아주 사소한 초기조건의 변화에도 민감하게 반응하며 이를 예측하는 것은 아주 어렵다. 나비 한 마리의 날개짓이 폭풍이 될지도 모르는 것이다.
설명
바이퍼케이션
위 그림은 로렌츠 어트랙터의 바이퍼케이션 다이어그램이다2. 가로축은 $25$ 부터 $325$ 까지 조금씩 변하는 $\rho$ 의 값, 세로축은 충분히 긴 시간 후의 궤도에서 $z$ 방향의 극댓값을 찍어놓았다. $r= \rho$ 가 작을 땐 매우 불안한 움직임을 보이지만, 일정 수준 이상으로 커지자 안정된 궤도를 가짐을 알 수 있다. 다음의 설명을 읽길 권장한다:
위 그래프는 $\gamma=100$ 일 때의 궤적을 그린 것이다. 보통 $\gamma$ 가 어느정도 이상 크다면 초기에는 카오틱한 움직임을 보이는 듯 하지만 갈수록 피리어딕해질 수 있다.
위 그래프는 초기의 궤적을 제거하고 긴 시간이 지난 뒤 입자의 궤적만을 나타내는 것이다. 완전히 피리어딕하지는 않지만, 거의 행동을 예측할 수 있는 수준이 되었다.
궤적을 보기 좋게 $xz$ 그래프를 그려보면 위과 같은 모양을 한다. 분명 처음엔 카오틱해보였으나 시간이 갈수록 궤도가 안정되어가고, 극댓값은 세 개 정도로 뚜렷하게 나타난다. 실제로 바이퍼케이션 다이어그램을 보아도 똑같은 해석을 내놓을 수 있다.
한편 $\gamma = 28$ 일 때는 충분히 긴 시간이 지나도 다음과 같이 카오틱한 움직임을 보인다. 바이퍼케이션 다이어그램에서는 여러 점이 흩어져서 두꺼워진 모양새를 이룬다.
반면 $\gamma = 300$ 일 때는 다음과 같이 뚜렷하게 나타나는 두 개의 극댓값만을 가진다. 바이퍼케이션 다이어그램에서는 단 한 $\gamma$ 값에 대해 두 점이 대응됨으로써 두 개의 선이 나타나게 된다.
랴푸노프 스펙트럼
다음은 $\rho \in [0, 100]$ 에서의 랴푸노프 스펙트럼을 나타낸다3.
코드
매트랩
다음은 매트랩을 통해 로렌츠 어트랙터의 궤적을 그려주는 코드다.
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 DataFrames, CSV, Plots, LinearAlgebra, ProgressMeter
using Base.Threads: @threads
function lorenz(v::AbstractVector; ρ = 28)
x, y, z = v
dx = 10*(y - x)
dy = x*(ρ - z) - y
dz = x*y - (8/3)*z
return [dx, dy, dz]
end
J_lorenz(x,y,z,ρ) = [
-10 10 0
ρ-z -1 -x
y x -8/3
]
function lyapunov_lorenz()
r_range = 0:0.1:100
schedules = DataFrame(idx = eachindex(r_range), r = r_range)
result = DataFrame(ρ = Float64[], λ1 = Float64[], λ2 = Float64[], λ3 = Float64[])
h = 1e-4; tend = 10_000_000
@showprogress @threads for dr = eachrow(schedules)[1:1:end]
Q = I(3); v = ones(3); λ = zeros(3);
for i in 1:tend
v = v + h*lorenz(v, ρ = dr.r); # push!(v_, v)
J = J_lorenz(v..., dr.r)
Q = (I(3) + h*J)*Q
Q, R = qr(Q)
λ += R .|> abs |> diag .|> log
end
λ ./= h*tend
push!(result, [dr.r, λ...])
end
sort!(result, :ρ)
plot(xticks = 0:20:100, legend = :none, size = [600, 300])
plot!(result.ρ, result.λ1, lw = 2, color = 1)
plot!(result.ρ, result.λ2, lw = 2, color = 2)
plot!(result.ρ, result.λ3, lw = 2, color = 3)
png("lyapunov")
end
lyapunov_lorenz()
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 ↩︎