今回は，Ｒを使った微分方程式の数値計算についての説明というよりも，数値計算結果を3次元プロットする方法の紹介です．私が以下で使っている方法は，ベストではないかもしれませんが，とりあえず，3次元空間内の軌道を立体的に描くことができます．結果はこんな感じです．
使ったスクリプトは以下です．事前にdeSolveパッケージをRにインストールしてください．
# 数値計算にdeSolveパッケージを使いました
#install.packages("deSolve")
require(deSolve)
### ローレンツモデルの定義 ###################
Lorenz.eq <- function(t,state,parameters) {
with(as.list(c(state, parameters)), {
dxdt <- sigma*(y-x)
dydt <- rho*x-y-x*z
dzdt <- -beta*z+x*y
list(c(dxdt, dydt, dzdt))
})
}
# よく使われている制御変数の値を使用
control.parms <- c(rho=28,sigma=10,beta=8/3)
# 初期条件はトランジェントを捨てずにすませるためアトラクタに近い点をとった．
ini.cond <- c(x=-2.292091, y=0.09829925, z=24.50526)
# 時間ステップの設定
time.step <- seq(0, 100, 0.01)
#######################################
# deSolveを使って数値計算
PROF <- ode(y=ini.cond,times=time.step,func=Lorenz.eq,parms=control.parms)
###
# グラフの描画
x.range <- 25
x <- y <- c(-x.range, x.range)
z <- matrix(rep(0,4), nrow=2)
###
par(mfrow=c(3,1))
plt1 <- persp(x,y,z,theta=30,xlim=c(-x.range,x.range),ylim=c(-x.range,x.range),zlim=c(0,2*x.range), phi=10,expand=1,cex.lab=1.4)
lines(trans3d(PROF[,"x"],PROF[,"y"],rep(0,length(PROF[,"z"])), pmat=plt1), col = 8)
lines(trans3d(PROF[,"x"],PROF[,"y"],PROF[,"z"], pmat=plt1), col = 2)
#
plt1 <- persp(x,y,z,theta=25,xlim=c(-x.range,x.range),ylim=c(-x.range,x.range),zlim=c(0,2*x.range), phi=45,expand=1,cex.lab=1.4)
lines(trans3d(PROF[,"x"],PROF[,"y"],rep(0,length(PROF[,"z"])), pmat=plt1), col = 8)
lines(trans3d(PROF[,"x"],PROF[,"y"],PROF[,"z"], pmat=plt1), col = 2)
#
plt1 <- persp(x,y,z,theta=20,xlim=c(-x.range,x.range),ylim=c(-x.range,x.range),zlim=c(0,2*x.range), phi=90,expand=1,cex.lab=1.4)
lines(trans3d(PROF[,"x"],PROF[,"y"],rep(0,length(PROF[,"z"])), pmat=plt1), col = 8)
lines(trans3d(PROF[,"x"],PROF[,"y"],PROF[,"z"], pmat=plt1), col = 2)
