ナンクル力学系

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

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

with 2 comments

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() [/sourcecode]

じゃあ、基本に戻って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() [/sourcecode]

同じく非線形振動子の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() [/sourcecode]

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

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

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() [/sourcecode]

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

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() [/sourcecode]

疲れたので、 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() [/sourcecode]

Advertisements

Written by tkf

February 27, 2009 at 10:27 pm

Posted in PC, 数学

Tagged with ,

2 Responses

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

  2. 液晶が割れていますが、正常に起動します。付属品や箱等すべてあります。査定できますか?

    wodyxwnqk@gmail.com

    October 19, 2013 at 3:48 pm


Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: