Rで力学系の解軌道を描画する

Rでは簡単に力学系の軌道を描くことができます。

rglでの軌道を描画に関しては@sfchaosさんのブログレプリカ交換モンテカルロ法を用いた力学系の軌道・パラメータ探索ですでに紹介されており、
またR-bloggerではggplot2をつかってきれいな軌道を書く例が紹介されています。
Coding a dynamic systems and controlling it via a graphical user interface | (R news & tutorials)

しかしコードも短く気軽にかけるので以下に張っておきます。

例としてLorenzアトラクタを用います。
loz.f.R

library(deSolve)
loz.params <- list(dl = 10,
               bt = 8/3,
               ro = 28)
loz <- function(t, y, p){
    dl <- p[['dl']]
    ro <- p[['ro']]
    bt <- p[['bt']]
    list(c(dl * (y[2] - y[1]),
           y[1] * (ro - y[3]) - y[2],
           y[1] * y[2] - bt * y[3] ))
}

これをanimationパッケージを使って解軌道の発展を可視化したものが以下のコードになります。
一定期間内の解をlineでつないだ静止画像の開始時間をずらしたものをまとめて動画にしています。

library(animation)
source("loz.f.R")

f <- lotz
lim <- 35
size <- 30
len <- 4

par <- lotz.params
stp <- 100

main <- function(){
  iniseq <-seq(-10,10,length=len)
  inigrid <- expand.grid(iniseq,iniseq,iniseq)
  yi <- c(0.1,0.1,0.1)
  
#  for (tmax in seq(100,3000,by=stp)){
  tmax <- 1000
  for (it in 1:100){
    times <- c(0, 0.01*(1:tmax))

    y <- lsoda(yi, times, f, par)

    scatterplot3d(y[,2:4],xlim=c(-lim,lim),
                  ylim=c(-lim,lim),
                  zlim=c(0,2*lim),
                  type='l',color="blue")
    yi <- y[stp,2:4]
#    yi <- y[nrow(y),2:4]
  }
  
}

saveMovie({
   main()
 }, movie.name = "loz.gif", interval = 0.1, nmax = 30, ani.width = 600, ani.height = 600)

(図をクリックして開かないとgitアニメにならないみたいです)
またrglのwebGL機能を使うとマウスでぐりぐり回転させることができます。

library(rgl)
source("loz.f.R")

times <- c(0, 0.01*(1:2000))
f <- lotz

lim <- 35
size <- 30
len <- 4

iniseq <-seq(-10,10,length=len)
inigrid <- expand.grid(iniseq,iniseq,iniseq)

par <- lotz.params

col <- c()
y <- c()
for (i in 1:nrow(y0s)){
                yt <- lsoda(y0s[i,], times, f, par)
                y <- rbind(y,yt)
                col <- c(col,rep(i,length(times)))
}

#open3d()                                                                                                                                                                                                                    
plot3d(y[,2:4],type='l',col=col,size=c(600,1000))
writeWebGL(width=700, height=550)

実行するとwebGLというフォルダにhtmlファイルができます。
現状中心を移動させる方法がわからないのが難点です。

感想
lorenz系はキャッチーでこういった例にはうってつけなのですがその性質は意外と複雑ですね。