로렌츠 어트랙터 📂동역학

로렌츠 어트랙터

Lorenz attractor

개요

로렌츠 방정식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$: 위키피디아1에 따르면 차원에 관련된 파라미터인 것 같다. 교재에서도 이것이 무엇인지 정확하게 설명하지 않는 경우가 많다.

처음 시뮬레이션을 할 때 이 파라미터들은 관례적으로 로렌츠가 했던대로 $\sigma = 10$, $\rho = 28$, $\displaystyle \beta = {{8} \over {3}}$ 와 같이 두는 경우가 많다. 레일리 수가 $\rho \approx 24.74$ 을 넘어가면 로렌츠 어트랙터가 묘사하는 입자의 움직임은 캐어릭Chaotic해진다. 아주 사소한 초기조건의 변화에도 민감하게 반응하며 이를 예측하는 것은 아주 어렵다. 나비 한 마리의 날개짓이 폭풍이 될지도 모르는 것이다.

바이퍼케이션

20190207\_143724.png

위 그림은 로렌츠 어트랙터의 바이퍼케이션 다이어그램이다. 가로축은 $25$ 부터 $325$ 까지 조금씩 변하는 $\rho$ 의 값, 세로축은 충분히 긴 시간 후의 궤도에서 $z$ 방향의 극댓값을 찍어놓았다. $r= \rho$ 가 작을 땐 매우 불안한 움직임을 보이지만, 일정 수준 이상으로 커지자 안정된 궤도를 가짐을 알 수 있다. 다음의 설명을 읽길 권장한다:

20190207\_142108.png

위 그래프는 $\gamma=100$ 일 때의 궤적을 그린 것이다. 보통 $\gamma$ 가 어느정도 이상 크다면 초기에는 카오틱한 움직임을 보이는 듯 하지만 갈수록 피리어딕해질 수 있다.

20190207\_142944.png

위 그래프는 초기의 궤적을 제거하고 긴 시간이 지난 뒤 입자의 궤적만을 나타내는 것이다. 완전히 피리어딕하지는 않지만, 거의 행동을 예측할 수 있는 수준이 되었다.

20190207\_160923.png

궤적을 보기 좋게 $xz$ 그래프를 그려보면 위과 같은 모양을 한다. 분명 처음엔 카오틱해보였으나 시간이 갈수록 궤도가 안정되어가고, 극댓값은 세 개 정도로 뚜렷하게 나타난다. 실제로 바이퍼케이션 다이어그램을 보아도 똑같은 해석을 내놓을 수 있다.

한편 $\gamma = 28$ 일 때는 충분히 긴 시간이 지나도 다음과 같이 카오틱한 움직임을 보인다. 바이퍼케이션 다이어그램에서는 여러 점이 흩어져서 두꺼워진 모양새를 이룬다.

20190207\_160633.png

반면 $\gamma = 300$ 일 때는 다음과 같이 뚜렷하게 나타나는 두 개의 극댓값만을 가진다. 바이퍼케이션 다이어그램에서는 단 한 $\gamma$ 값에 대해 두 점이 대응됨으로써 두 개의 선이 나타나게 된다.

20190207\_161132.png

코드

매트랩

다음은 매트랩을 통해 로렌츠 어트랙터의 궤적을 그려주는 코드다.

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))

  1. https://en.wikipedia.org/wiki/Lorenz_system#Overview ↩︎

댓글