로렌츠 어트랙터
Lorenz attractor
개요
로렌츠 방정식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해진다. 아주 사소한 초기조건의 변화에도 민감하게 반응하며 이를 예측하는 것은 아주 어렵다. 나비 한 마리의 날개짓이 폭풍이 될지도 모르는 것이다.
바이퍼케이션
위 그림은 로렌츠 어트랙터의 바이퍼케이션 다이어그램이다. 가로축은 $25$ 부터 $325$ 까지 조금씩 변하는 $\rho$ 의 값, 세로축은 충분히 긴 시간 후의 궤도에서 $z$ 방향의 극댓값을 찍어놓았다. $r= \rho$ 가 작을 땐 매우 불안한 움직임을 보이지만, 일정 수준 이상으로 커지자 안정된 궤도를 가짐을 알 수 있다. 다음의 설명을 읽길 권장한다:
위 그래프는 $\gamma=100$ 일 때의 궤적을 그린 것이다. 보통 $\gamma$ 가 어느정도 이상 크다면 초기에는 카오틱한 움직임을 보이는 듯 하지만 갈수록 피리어딕해질 수 있다.
위 그래프는 초기의 궤적을 제거하고 긴 시간이 지난 뒤 입자의 궤적만을 나타내는 것이다. 완전히 피리어딕하지는 않지만, 거의 행동을 예측할 수 있는 수준이 되었다.
궤적을 보기 좋게 $xz$ 그래프를 그려보면 위과 같은 모양을 한다. 분명 처음엔 카오틱해보였으나 시간이 갈수록 궤도가 안정되어가고, 극댓값은 세 개 정도로 뚜렷하게 나타난다. 실제로 바이퍼케이션 다이어그램을 보아도 똑같은 해석을 내놓을 수 있다.
한편 $\gamma = 28$ 일 때는 충분히 긴 시간이 지나도 다음과 같이 카오틱한 움직임을 보인다. 바이퍼케이션 다이어그램에서는 여러 점이 흩어져서 두꺼워진 모양새를 이룬다.
반면 $\gamma = 300$ 일 때는 다음과 같이 뚜렷하게 나타나는 두 개의 극댓값만을 가진다. 바이퍼케이션 다이어그램에서는 단 한 $\gamma$ 값에 대해 두 점이 대응됨으로써 두 개의 선이 나타나게 된다.
코드
매트랩
다음은 매트랩을 통해 로렌츠 어트랙터의 궤적을 그려주는 코드다.
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))