1. 導入

状態がある決定論的法則に従って時間的に変化していくような系を一般的に決定論的力学と呼ぶ。力学系における運動状態を記述するためには、その運動状態を記述するのに必要な変数（状態変数）を考えることが必要である。この状態変数によって張られた空間を状態空間または相空間と呼ぶ。力学系の決定論的立場は、ある時刻における系の状態を\(x\)とすると、「状態\(x\)が次に状態空間中をどの方向にどれくらいの速さ（\(dx / dt\)）で変化するかが現在の状態\(x\)により唯一に確定している」、言い換えれば、「いかなる系も初期状態が決まればその未来は予測できる」という事実を前提にしている。ここで、位置\(x\)に依存するベクトル量が非線形関数で表される場合には、その挙動は非常に複雑かつ多様な様相を示すことがあり、現在では決定論的カオスと呼ばれている。なお今回、決定論的力学系において、数値計算によって観測される数値的カオスから生じるストレンジアトラクターを取り上げる。

ストレンジアトラクターとは、まず,

1.アトラクター：散逸系(エネルギーの出入りがある系、エネルギー非保存系)の運動は十分な時間が経つと特定の軌道や点に落ち着く。この運動の過渡状態後の安定した状態をアトラクターと言う。アトラクタは運動の複雑さや種類によって様々な形をとるが、基本的には平衡点アトラクター、周期アトラクター、準周期アトラクター、ストレンジアトラクターの４つに分類される。

2.カオス：決定論的な系において、初期条件の僅かな違いが瞬く間に増幅され、 予測不可能な挙動を示す現象をいう。

これらの用語を踏まえて、非整数次元のアトラクターやカオス理論でしか振る舞いを説明できない力学系のアトラクターをストレンジであると（非形式的に）いい、元はカオスアトラクターと呼ばれていたが、ダヴィッド・ルエールとFloris Takensが流体の力学系における一連の分岐の結果として生じるアトラクターを指してストレンジアトラクターという造語を使用した。また、ストレンジアトラクターという場合、集合（あるいは集合上での動き）を図示することが困難である場合であり、カントール集合と非可算無限集合の直積構造を持つことが多い。

ストレンジアトラクターの中でも３自由度常微分方程式で表される連続時間システムである以下の２つの方程式を具体例として扱う。

2. ローレンツ方程式

ストレンジアトラクターの比較的有名な例としてLorenz方程式がある。アメリカの気象学者のE. Lorenzが1963年に、大気の熱対流現象研究のために提案したモデルを Lorenz方程式と呼ぶ。本来は偏微分方程式で記述される熱対流現象を大幅な近似のもとで簡単化して導かれたものであり、３自由度常微分方程式で表される連続時間システムである。以下の式で表されるLorenz方程式によって得られるアトラクターをLorenzアトラクターと呼んでいる。これにより、E. Lorenzは気象の長期予報は不可能であることを証明した。

ローレンツ方程式は、以下の非線形常微分方程式で表される。 \[\begin{equation} \left\{ \begin{aligned} \frac{dx}{dt} &= -px + -py \\ \frac{dy}{dt} &= -xz + rx - y \\ \frac{dz}{dt} &= xy + bz \end{aligned} \right. \end{equation}\]

ただし \(p,r,b\) は定数。

まずは、ローレンツ方程式を解析をするためにルンゲクッタ法のコードを準備する。

rk4 <- function(dx.dt=function(t,x){}, h=1E-7,
                x0=1, start=0, end=1){
  
  
  t <- seq(start,end,by=h)
  nsteps <- length(t)-1
  nvar <- length(x0)
  
  m <- data.frame(matrix(0, nrow = nsteps+1, ncol = nvar+1))
  m[1,] <- c(start,x0)
  m[,1] <- t
  
  for(i in 1:nsteps){ #ルンゲクッタ法
    
    s1 <- dx.dt(t[i],c(t(m[i,-1])))
    s2 <- dx.dt(t[i+h/2],c(t(m[i,-1]))+h*s1/2)
    s3 <- dx.dt(t[i+h/2],c(t(m[i,-1]))+h*s2/2)
    s4 <- dx.dt(t[i+h],c(t(m[i,-1]))+h*s3)
    
    m[i+1,-1] <- m[i,-1] + h/6*(s1+2*s2+2*s3+s4)
  }
  return(m)
  
}

この下で、解析を進めていく。

変数 \(x,y,z\) の初期値を

x_0 <- 0
y_0 <- 1
z_0 <- 2

とし、定数 \(p,q,r\) を

p <- 10
r <- 28
b <- 8/3

と定める。このことにより、ストレンジアトラクターが生じる。

さて、３変数ベクトル値関数を\(f\) として

f <- function(t,a){
  x <- a[1]
  y <- a[2]
  z <- a[3]
  
  dx <- -p*x + p*y
  dy <- -x*z + r*x - y
  dz <- x*y - b*z
  
  return(c(dx, dy, dz))
}

と定め、時間 \(t\) を \(0.02\) 毎に区切って解くと、

h <- 0.02
x0 <- c(x_0, y_0, z_0)

interval <- 0.02
skip <- round(interval/h)

start <- 0
#end <- 1 
end <- 80

df <- rk4(dx.dt=f, h=h, x0=x0, start=start, end=end)

head(df)

##     X1        X2       X3       X4
## 1 0.00 0.0000000 1.000000 2.000000
## 2 0.02 0.1824975 1.028832 1.897992
## 3 0.04 0.3455760 1.145772 1.805039
## 4 0.06 0.5079994 1.344684 1.721650
## 5 0.08 0.6848621 1.628108 1.649483
## 6 0.10 0.8896762 2.005817 1.591629

となる。また、グラフは、

library(scatterplot3d)

を読み込んで、

scatterplot3d(df[,2],df[,3],df[,4],type="l",xlab = "Lorenz方程式(p = 10, r = 28, b = 8/3)",ylab = "",zlab = "")

と表せる。

3. ローレンツ方程式：シュミレーション

ローレンツ方程式を解いた結果を動画にする。

len <- length(df[,1])
dur <- 30
 for(i in 1:(len-dur)){
    if(i %% skip == 0){
      scatterplot3d(df[1:(i+dur),2],df[1:(i+dur),3],df[1:(i + dur),4], 
                    xlim = c(-20, 20), ylim = c(-20, 30), zlim = c(0,50),
                    type="l")
    }}

4. ラングフォード方程式

数学者Langfordが、不変トーラスからカオスへの分岐を研究するために用いた方程式であり、以下の式のように３変数の微分方程式として表される。 \[\begin{equation} \left\{ \begin{aligned} \frac{dx_{2}}{dt} &= (z_{2} - u)x_{2} - wy_{2} \\ \frac{dy_{2}}{dt} &= wx_{2} + (z_{2} - u)y_{2} \\ \frac{dz_{2}}{dt} &= n + dz_{2} - \frac{z_{2}^2}{3} - (x_2^2 + y_{2}^2)(1 + lz_{2}) + ez_{2}x_{2}^3 \end{aligned} \right. \end{equation}\]

ただし, \(d,u,n,w,l,e\) は定数。

同様の操作を施していく。

変数 \(x_{2},y_{2},z_{2}\) の初期値を

x2_0 <- 1
y2_0 <- 1
z2_0 <- 1

とし、定数\(d,u,n,w,l,e\)を

d <- 1
u <- 0.7
n <- 0.6
w <- 3.5
l <- 0.25
e <- 0

と定める。

f <- function(t,r){
  x2 <- r[1]
  y2 <- r[2]
  z2 <- r[3]
  
  dx2 <- z2*x2 - u*x2 - w*y2
  dy2 <- w*x2 + z2*y2 - u*y2
  dz2 <- n + d*z2 - z2^3/3 - x2^2 - y2^2 - l*z2*x2^2  - l*z2*y2^2+ e*z2*x2^3
  
  return(c(dx2, dy2, dz2))
}

h <- 0.02
x0 <- c(x2_0, y2_0, z2_0)

interval <- 0.02
skip <- round(interval/h)

start <- 0
#end <- 1 
end <- 80

df <- rk4(dx.dt=f, h=h, x0=x0, start=start, end=end)

head(df)

##     X1        X2       X3        X4
## 1 0.00 1.0000000 1.000000 1.0000000
## 2 0.02 0.9329593 1.073652 0.9751613
## 3 0.04 0.8600856 1.142260 0.9499972
## 4 0.06 0.7817878 1.205326 0.9245366
## 5 0.08 0.6985205 1.262388 0.8988094
## 6 0.10 0.6107808 1.313024 0.8728466

scatterplot3d(df[,2],df[,3],df[,4],type="l",xlab = "Langford方程式(d =  1, u = 0.7, n = 0.6, w = 3.5, l = 0.25, e = 0)", ylab = "",zlab = "")

と表せる。

5. ラングフォード方程式：シュミレーション

ラングフォード方程式を解いた結果を動画にする。

len <- length(df[,1])
dur <- 30
for(i in 1:(len-dur)){
    if(i %% skip == 0){
      scatterplot3d(df[1:(i+dur),2],df[1:(i+dur),3],df[1:(i + dur),4], 
                    xlim = c(-2, 2), ylim = c(-2, 2), zlim = c(-1,2),
                    type="l")
    }
  }