ナンクル力学系

学んだ事を書き連ねていこう。

Rの練習のために力学系を描いてみた

with one comment

Rを使う練習。使い方間違ってる気がしてならない。

odesolveとscatterplot3dが要ります。

まずは、NN屋さんっぽくCTRNN

3つのニューロン使えば、カオスになるやつがあるよ、という例。

library(odesolve)
library(scatterplot3d)

sigmoid <- function(x){
  (tanh(x/2) + 1) / 2
}

dydt <- function(t, y, p){
  wt <- p[['wt']]
  bs <- p[['bs']]
  tu <- p[['tu']]

  list(( - y + wt %*% sigmoid( y +  bs ) ) / tu)
}

params <- list(wt = matrix(
                 c( 5.422 , -0.24 , 0.535 ,
                   -0.018 , 4.59  , -2.25 ,
                   2.75   , 1.21  , 3.885 ), 3,3),
               bs = c(-4.108, -2.787, -1.114),
               tu = c(1,2.5,1))
times <- c(0, 0.1*(1:100000))
y <- lsoda(c(1,0,0), times, dydt, params)

png(filename="rode-ctrnn.png")
scatterplot3d(y[,2:4], type='l')
dev.off()

じゃあ、基本に戻ってVan der Polだ！

初期値を色々変えてプロットしてるけど、全部ひとつのアトラクタに落ち込むで しょ、という例。力学系やるとまず初めに感動するやつ（多分）。もう少し綺麗 にかけないのかなこれw

library(odesolve)

dydt <- function(t, y, p){
  c <- p[['c']]
  F <- p[['F']]

  list( c(y[2],
          c * (1 - y[1]^2) *y[2] - y[1] + F * cos(t) ))
}

params <- list(c=1.0, F=0)
times <- c(0, 0.05*(1:1000))

epg = expand.grid(seq(-3,3,length=6),seq(-3,3,length=6))
y0s = matrix(c(epg[,1],epg[,2]), 36,2)

png(filename="rode-vdp.png")
plot(0, type='l', xlim=c(-3,3), ylim=c(-3,3))
apply(y0s, 1,
      function(y0){
        print(y0)
        y <- lsoda(y0, times, dydt, params)
        lines(y[,2],y[,3], type='l')
      })
dev.off()

同じく非線形振動子のDuffing方程式

こいつは、リミットサイクル持たずに平衡点を持つ。ハミルトン系+減衰項って かたちなので当然。

library(odesolve)

dydt <- function(t, y, p){
  c <- p[['c']]
  F <- p[['F']]

  list( c( y[2],
          - c * y[2] + y[1] - y[1]^3 + F * cos(t)
          ))
}

params <- list(c=0.25, F=0.0)
times <- c(0, 0.1*(1:1000))

len = 4
epg = expand.grid(seq(0.5,2.5,length=len),seq(0.5,2.5,length=len))
y0s = matrix(c(epg[,1],epg[,2]), 16,2)

png(filename="rode-duffing.png")
plot(0, type='l', xlim=c(-3,3), ylim=c(-3,3))
apply(y0s, 1,
      function(y0){
        #print(y0)
        y <- lsoda(y0, times, dydt, params)
        lines(y[,2],y[,3], type='l')
      })
dev.off()

エネルギーもプロットしてみた例

頑張ったけど、線ではプロットできなかった。。。 二つの凹みがある容器みたいな ポテンシャル面を沿うような形で降りてきている（減衰）していることが分かる。

library(odesolve)
library(scatterplot3d)

dydt <- function(t, y, p){
  c <- p[['c']]
  F <- p[['F']]

  list( c( y[2],
          - c * y[2] + y[1] - y[1]^3 + F * cos(t)
          ))
}

eng <- function(y1,y2){
  return( y2^2 / 2.0 - y1^2 / 2.0 + y1^4 / 4.0 )
}

params <- list(c=0.25, F=0.0)
times <- c(0, 0.1*(1:1000))

len = 4
epg = expand.grid(seq(0.5,2.5,length=len),seq(0.5,2.5,length=len))
y0s = matrix(c(epg[,1],epg[,2]), 16,2)

y1 <- c()
y2 <- c()
ye <- c()
clr <- c()
cls <- colors()
dcls <- round(length(cls)/(len*len))
for (i in 1:(len*len)){
  y <- lsoda(y0s[i,], times, dydt, params)
  y1 <-  c(y1, y[,2])
  y2 <-  c(y2, y[,3])
  ye <-  c(ye, eng(y[,2],y[,3]))
  clr <- c(clr, rep(cls[(i-1)*dcls+1], length(times)))
}

png(filename="rode-duffing-we.png")
scatterplot3d(y1,y2,ye, type='p', color=clr, angle=20)
dev.off()

周期外力を加えるとカオスになるね！

library(odesolve)

dydt <- function(t, y, p){
  c <- p[['c']]
  F <- p[['F']]

  list( c( y[2],
          - c * y[2] + y[1] - y[1]^3 + F * cos(t)
          ))
}

params <- list(c=0.25, F=0.4)
times <- c(0, 0.1*(1:5000))

y <- lsoda(c(0.5,0.5), times, dydt, params)

png(filename="rode-duffing-chaos.png")
plot(y[,2:3], type='l')
dev.off()

疲れたので、 Lorentz Attractor でも書いて今日はおしまいにするか

ちゃんちゃん。

library(odesolve)
library(scatterplot3d)

dydt <- 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] ))
}

params <- list(dl = 10,
               bt = 8/3,
               ro = 28)
times <- c(0, 0.01*(1:10000))
y <- lsoda(c(1,0,0), times, dydt, params)

png(filename="rode-lorenz.png")
scatterplot3d(y[,2:4], type='l')
dev.off()

参考

R:
力学系:

Possibly related posts: (automatically generated)

Written by tkf

February 27, 2009 at 10:27 pm

Posted in PC, 数学

Tagged with ,

«
»

One Response

Subscribe to comments with RSS.

  1. 色々コメントしてみるテスト。

    3次元plotにscatterplot3dを使ってるけど、persp関数を使う手もあります。こんな感じでplotできます。

    x <- seq(-3,3,length.out=30)
    y <- seq(-3,3,length.out=30)

    persp(x,y,outer(x,y,function(i,j){
    exp(-(i^2+j^2))}),
    theta=320,phi=20,col=rainbow(50),ticktype=”detailed”,
    xlab=”x”,ylab=”y”,zlab=”density”)

    scatterplot3dが悪いというわけではないけど、なんとなく俺はpersp関数を使ってますwww。あとpersp関数なら以下のリンクのような感じで線を引くことができます。
    http://d.hatena.ne.jp/syou6162/20080914/1221307975

    あとなんか細かいところを。
    times <- c(0, 0.1*(1:100000))
    って書いてあるんだけど、望みのものはたぶんこれでどんぴしゃだと思います。
    seq(from=0,to=10,length.out=1000)
    …あれ、でも
    seq(-3,3,length=6))
    とか書いてあるのか。。。

    てか、Rの練習でこのレベルとかすげえwww俺はapplyまともに理解したのはR使い始めてから一年半後でした＞＜Tsukuba.Rでは期待してます＞＜

    syou6162

    February 28, 2009 at 4:17 am


Leave a Reply