1 ローレンツ方程式

ローレンツ方程式とは，カオス的振る舞いを示す非線形方程式の1つである． \[\begin{align*} \dfrac{dx}{dt}=&-px+py\\ \dfrac{dy}{dt}=&-xz+rx-y\\ \dfrac{dz}{dt}=&xy-bz \end{align*}\] ローレンツ(Lorentz)の1963年の論文の中で， \(p=10, r=28, b=8/3\)という設定が紹介されている．

2 シミュレーション

以下ではこのローレンツ方程式のシミュレーションを行う． 単純のためライブラリを使うが，レポートではできる限りライブラリを使わず，アルゴリズムの勉強をしてほしい．

必要なライブラリを読み込む．

library(deSolve)

## 
## Attaching package: 'deSolve'
## 
##  以下のオブジェクトは 'package:graphics' からマスクされています: 
## 
##      matplot

library(rgl)

パラメーターを設定する．

parameters <- c(s=10, b=8/3, r=28)
initial <- c(x=0, y=0.1, z=0)
times <- seq(0, 20, 0.1)

ローレンツ方程式を解いて，その結果をoutに出力する．

Lorentz <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dx <- s * (y-x)
    dy <- r * x - y - x * z
    dz <- -b * z + x * y
    list(c(dx,dy,dz))
  })
}
out <- ode(y = initial, times = times, func = Lorentz, parms = parameters)

その結果をcx,cy,czにコピーしながら動画として再生する．

imax <- nrow(out)
cx <- rep(NA,imax)
cy <- rep(NA,imax)
cz <- rep(NA,imax)

library(scatterplot3d)
xlim <- c(min(out[,"x"]),max(out[,"x"]))
ylim <- c(min(out[,"y"]),max(out[,"y"]))
zlim <- c(min(out[,"z"]),max(out[,"z"]))
for(i in 1:imax){
    cx[i] <- out[i,"x"]
    cy[i] <- out[i,"y"]
    cz[i] <- out[i,"z"]
    scatterplot3d(cx,cy,cz,xlim=xlim,ylim=ylim,zlim=zlim,type="l")
}