今年度は研究室のセミナーで，カオスの勉強をしています．今回は，ロジスティック写像の分岐図を描くRスクリプトを載せておきます．まずは，ロジスティック写像について簡単に説明します．
ロジスティック写像 (logistic map)
ロジスティック写像は，以下の差分方程式で定義されます．
ここで，，
は制御変数で，
の範囲で考えることが多いです．R.L. デバネーの「カオス力学系入門」では，
で，不安定なカオス軌道の存在を示したりしますが，通常は
を考えます．
ロジスティック写像の振る舞いについては，私が作った以下の動画を参考にしてください．検索すると，もっと詳しい説明を見つけられると思います．
youtu.be
Rで分岐図
以下のRスクリプトを実行すると，ロジスティック写像の分岐図を描くことができます．
###############
# 制御変数の範囲
r.range <- c(2.8, 4)
# 変化させるrの総数
n.r <- 500
###############
r.list <- seq(r.range[1],r.range[2],length=n.r)
###
# xの初期値
x0 <- 0.5
###
plot(c(),c(),xlim=r.range,ylim=c(0,1),xaxs="i",yaxs="i",xlab="r",ylab=expression({x}[n]))
###
# 過渡状態ステップ数
n.trans <- 500
# 周期の確認ステップ数
n.per <- 256
# 非周期の場合の追加ステップ数
n.plt <- 1000
for(r in r.list){
x.tmp <- x0
for(i in 1:n.trans){
x.tmp <- r*x.tmp*(1-x.tmp)
}
xn <- c()
xn[1] <- x.tmp
for(i in 2:n.per){
xn[i] <- r*xn[i-1]*(1-xn[i-1])
}
xn.u <- unique(xn)
n.xn <- length(xn.u)
if(n.xn <= 64){ # 64周期以下であることを判定
points(rep(r,n.xn),xn.u,"p",pch=16,cex=0.05,col=rgb(0,0.5,1,alpha=1))
}else{
for(i in (n.per+1):(n.per+n.plt)){
xn[i] <- r*xn[i-1]*(1-xn[i-1])
}
n.xn <- length(xn)
points(rep(r,n.xn),xn,"p",pch=16,cex=0.05,col=rgb(0,0.5,1,alpha=0.05))
}
# 最後の値を初期値に設定
x0 <- tail(xn,1)
}このスクリプトを，少し変えて描いた分岐図の例が以下です．
分岐図を描く領域を変えて，分岐図に見られる自己相似性を調べてみてください．
