今日の講義資料です．散逸力学系のカオスアトラクター（下図）をＲを使って描画してみます．ただし，ここで説明する例は，ただの3次元プロットです．考える微分方程式はレスラー方程式です．パラメタは，カオスになる代表的な値を使います．方程式については，以下のWikipediaを参考にしてくださ（詳しく書いていないですが）．
レスラー方程式 - Wikipedia
Ｒを使って数値解を3次元プロット
以下のＲスクリプトの
plt1 <- persp(x,y,z,theta=30,phi=10,.....)
に含まれる，thetaとphiの値を変えれば，見る角度を変えられます．thetaで水平面内の角度，phiで上下方向の角度を変えられます．
# 数値計算にdeSolveパッケージを使いました
#install.packages("deSolve")
require(deSolve)
### レスラーモデルの定義 ###################
Lorenz.eq <- function(t,state,parameters) {
with(as.list(c(state, parameters)), {
dxdt <- -y - z
dydt <- x + a * y
dzdt <- b + x * z - c * z
list(c(dxdt, dydt, dzdt))
})
}
# よく使われている制御変数の値を使用
control.parms <- c(a=0.2,b=0.2,c=5.7)
# 初期条件はトランジェントを捨てずにすませるためアトラクタに近い点をとった．
ini.cond <- c(x=-4.2842519245, y=-2.438439520, z=0.01965899)
# 時間ステップの設定
time.step <- seq(0, 200, 0.01)
#######################################
# deSolveを使って数値計算
PROF <- ode(y=ini.cond,times=time.step,func=Lorenz.eq,parms=control.parms)
###
# グラフの描画
x.range <- 12
x <- y <- c(-x.range, x.range)
z0 <- 0
z <- matrix(rep(z0,4), nrow=2)
###
# thetaとphiの値を変えれば見る角度が変る
par(mfrow=c(2,2))
plt1 <- persp(x,y,z,theta=30,phi=10,xlim=c(-x.range,x.range),ylim=c(-x.range,x.range),zlim=c(z0,2*x.range),expand=1,cex.lab=1.4)
lines(trans3d(PROF[,"x"],PROF[,"y"],PROF[,"z"], pmat=plt1), col = 2)
#
plt1 <- persp(x,y,z,theta=80,phi=10,xlim=c(-x.range,x.range),ylim=c(-x.range,x.range),zlim=c(z0,2*x.range),expand=1,cex.lab=1.4)
lines(trans3d(PROF[,"x"],PROF[,"y"],PROF[,"z"], pmat=plt1), col = 2)
#
plt1 <- persp(x,y,z,theta=20, phi=10,xlim=c(-x.range,x.range),ylim=c(-x.range,x.range),zlim=c(z0,2*x.range),expand=1,cex.lab=1.4)
lines(trans3d(PROF[,"x"],PROF[,"y"],PROF[,"z"], pmat=plt1), col = 2)
#
plt1 <- persp(x,y,z,theta=-45,xlim=c(-x.range,x.range),ylim=c(-x.range,x.range),zlim=c(z0,2*x.range), phi=30,expand=1,cex.lab=1.4)
lines(trans3d(PROF[,"x"],PROF[,"y"],PROF[,"z"], pmat=plt1), col = 2)このスクリプトの実行結果が以下です．
