ローレンツ・アトラクター
概要
ローレンツ方程式は、大気の対流を連立常微分方程式で表した数学的モデルである。
システム
$$ \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$: 粘性と熱伝導率に関するパラメータであるプラントル数だ。
- $\rho$: 流体の熱伝達方法に関するパラメータであり、レイリー数だ。
- $\beta$: ウィキペディア1によると、次元に関連するパラメータらしい。教科書でもこれが何であるかを正確に説明しないことが多い。
初めのシミュレーションでは、これらのパラメータはローレンツが行ったように$\sigma = 10$、$\rho = 28$、$\displaystyle \beta = {{8} \over {3}}$ とすることが多い。レイリー数が $\rho \approx 24.74$ を超えると、ローレンツアトラクタが描く粒子の運動はカオスになる。非常に小さい初期条件の変化にも敏感に反応し、その予測は非常に困難だ。蝶が一回羽ばたくことで嵐が起こりうるということだ。
バイファーケーション
上の図はローレンツアトラクターのバイファーケーション図である。横軸は$25$から$325$まで少しずつ変わる$\rho$の値で、縦軸は十分に長い時間後の軌道で$z$方向の極大値をプロットしたものだ。$r= \rho$が小さいときは非常に不安定な動きを見せるが、一定レベル以上に大きくなると安定した軌道を持つようになることが分かる。次の説明を読むことをお勧めする:
上のグラフは、$\gamma=100$ の時の軌跡を描いたものである。通常、$\gamma$がある程度以上大きい場合、初めはカオス的な動きを見せるように思われるが、時間が経つにつれ周期的になることができる。
上のグラフは、最初の軌跡を取り除き、長い時間が経過した後の粒子の軌跡のみを示している。完全に周期的ではないが、ほぼ行動を予測できるレベルになった。
$xz$ のグラフを描くと、上のような形をしている。最初はカオス的に見えたが、時間が経つにつれて軌道が安定していき、極大値は3つほどにはっきりと現れる。実際にバイファーケーション図を見ても同じ解釈を出すことができる。
$\gamma = 28$ のときは、十分に長い時間が経っても次のようにカオス的な動きを見せる。バイファーケーション図では、多くの点が散らばって厚くなった形状を成す。
一方、$\gamma = 300$のときは、次のように明確に現れる2つの極大値を持つ。バイファーケーション図では、一つの$\gamma$値に対して2点が対応することで、2つの線が現れるようになる。
コード
MATLAB
以下はMATLABを使ってローレンツアトラクターの軌跡を描くコードである。
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)
Fortran
以下はFortranを使ってバイファーケーション図のソースを作り、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
ファイルを使ってバイファーケーション図を描くMATLABのコードだ。
[file,path] = uigetfile('*.csv');
data = csvread([path,file]);
sz = 0.6;
scatter(data(:,1),data(:,2),sz,'.');
Julia
以下はJuliaでローレンツアトラクターを見つけて描くコードである。
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))