ローレンツ方程式の実装とアトラクタの描画

lorenz-attractor.R

library(pforeach)
library(rgl) #3Dプロットに必要
library(ggplot2)
library(gridExtra)

x <- 0.1 #初期値
y <- 0.1 #初期値
z <- 0.1 #初期値

h <- 0.01 #増分
maxit <- 200 #hの最大値

p <- 10 #初期値
r <- 28 #初期値
b <- 8/3 #初期値

iter <- seq(0.01, maxit, by = h)

dx.dt <- function(xx, yy, zz, t) { #dx/dt の定義
    return(-p*xx + p*yy)
}

dy.dt <- function(xx, yy, zz, t) { #dy/dt の定義
    return(-xx*zz + r*xx - yy)
}

dz.dt <- function(xx, yy, zz, t) { #dz/dt の定義
    return(xx*yy - b*zz)
}

npforeach (t = iter, .c = rbind)({ #x, y, zそれぞれについてルンゲクッタ法を適用すればよい．
    k1.x <- dx.dt(x, y, z, t)
    k1.y <- dy.dt(x, y, z, t)
    k1.z <- dz.dt(x, y, z, t)
    
    k2.x <- dx.dt(x + h/2 * k1.x,
                  y + h/2 * k1.x,
                  z + h/2 * k1.x,
                  t + h/2)
    k2.y <- dy.dt(x + h/2 * k1.y,
                  y + h/2 * k1.y,
                  z + h/2 * k1.y,
                  t + h/2)
    k2.z <- dz.dt(x + h/2 * k1.z,
                  y + h/2 * k1.z,
                  z + h/2 * k1.z,
                  t + h/2)

    k3.x <- dx.dt(x + h/2 * k2.x,
                  y + h/2 * k2.x,
                  z + h/2 * k2.x,
                  t + h/2)
    k3.y <- dy.dt(x + h/2 * k2.y,
                  y + h/2 * k2.y,
                  z + h/2 * k2.y,
                  t + h/2)
    k3.z <- dz.dt(x + h/2 * k2.z,
                  y + h/2 * k2.z,
                  z + h/2 * k2.z,
                  t + h/2)

    k4.x <- dx.dt(x + h * k3.x,
                  y + h * k3.x,
                  z + h * k3.x,
                  t + h)
    k4.y <- dy.dt(x + h * k3.y,
                  y + h * k3.y,
                  z + h * k3.y,
                  t + h)
    k4.z <- dz.dt(x + h * k3.z,
                  y + h * k3.z,
                  z + h * k3.z,
                  t + h)
    
    x <- x + h/6 * (k1.x + 2*k2.x + 2*k3.x + k4.x)
    y <- y + h/6 * (k1.y + 2*k2.y + 2*k3.y + k4.y)
    z <- z + h/6 * (k1.z + 2*k2.z + 2*k3.z + k4.z)
    
    data.frame(x = x, y = y, z = z, t = t)
}) -> result

plot3d(result$x, result$y, result$z, type = 'l', lwd = 1, col = "orange")

ggplot(result, aes(x, z, color=y)) + geom_path(size = 0.5)