Rで力学系の軌道探索

先日のTokyoRで「(続)Rで力学系」というタイトルで発表させていただきました。

発表では残念ながら時間切れでStandard mapのホモクリニック軌道を重み付け(lyapunov weighted sampling)によって探索する手法の紹介までには至りませんでした。

やや不完全な部分がありますが使用したRのコードは以下になります。

library(deSolve)

standardmap_step<- function(y,K){
  p <- y[1]
  q <- y[2]
  p <- (p+K*sin(q)) %% (2*pi)
  q <- (q+p ) %% (2*pi)
  c(p,q)
}

standardmap_step<- function(y,K,del){
  p <- y[1]
  q <- y[2]
  p <- (p-K*del*sin(2*pi*q)/(2*pi)) 
  q <- (q+del*p ) %% 1
  c(p,q)
}

source("standardmap.f.R")

#debug <- T
times <- 100
fname <- "stdm_walker"


propotional_p<- function(inits,costs){
   inits.n <- c()

    l <- length(costs)
    for (j in 1:l){
      nn <- costs[j]

      if(nn>=1)
        for(kk in 1:nn){
          inits.n <- rbind(inits.n,inits[j,])
        }

      if((nn%%1)>runif(1,0,1)){
        inits.n <- rbind(inits.n,inits[j,])
      }
    }
   inits.n <- rbind(inits.n,inits.n)
   inits.n <-inits.n[1:l,]
}

calclyap <- function(x0,dx,times,f,a=1){
  x <- x0
  xd <-  x0+dx
  for (t in times){
    r <- runif(1,-dx,dx)
    x <- f(x)+c(r,0)
    xd <- f(xd)+c(r,0)
  }
  m <- abs(matrix(c(x-xd,1-(x-xd))))
  dxt <- max(apply(m,2,min))
#  dxt <- max(abs((x-xd)))
  p <- dxt/abs(dx)**a
#  pp <- log(p)/times
}

main <- function(){
  d <- 2
  N <- 200
  del <- 1#0.51
  K <- 1
  var <- 1e-4
  ites <- 200
  a <- -1
  
  f <- function(y){
    standardmap_step(y,K,del)
  }

  std.lyap.run <- function(i){
    calclyap(i,var,times,f,a)
  }

  lps <- c()
  inits <- c()

  inits <- matrix(runif(d*N,-1,1),nc=d)
  
#  for (i in 1:N){
#    inits <-rbind(inits,c(0,i/N))
#    inits <-rbind(inits,c(0.1,i/N))
#  }
#  inits <- as.matrix(inits)
  
  for (t in 0:ites){
    weight <- apply (inits,1,std.lyap.run)
    costs <- weight/sum(weight)*N

    print(max(weight))
    lps <-  rbind(lps,mean(weight))

    suf <- paste("N",N,"a",a,"K",K,"del",del,"v",var,sep="_")
    
    if(showplot && (t%%round(ites/10)==0)){
      print(weight)    

      ps <- c()
      qs <- c()
      ll <- min(N,40)
      for(i in 1:ll){
         y <- inits[i,]
         for (tt in 1:1000){
           ps <- rbind(ps,y[1])
           qs <- rbind(qs,y[2])
           y <- f(y)
         }
       }
      col <- ceiling((1:times)/(times/ll))%%255

      png(paste(fname,suf,"t",t,".png",sep="_"))
      plot(ps,qs,col=col,cex=0.2,
           xlim=c(-1,1), ylim=c(0,1),
           type='p',
           xlab="p",ylab="q(θ)")
      dev.off()

      rr <- 0.001
      png(paste(fname,"small",suf,"t",t,".png",sep="_"))

      plot(ps,qs,col=col,cex=0.2,
           xlim=c(-rr,rr), ylim=c(.5-rr,.5+rr),
           type='l',
           xlab="p",ylab="q(θ)")
      dev.off()

    }
    inits <- propotional_p(inits,costs)
    inits <- inits+matrix(rnorm(length(inits),0,var),nc=ncol(inits))
  } 

  png(filename=paste(fname,suf,"lyp",".png",sep="_"))
  plot(lps)
  dev.off()

元ネタ論文で紹介されている手法はparticle filterに似ているのですが、RのパッケージRcppSMCでは限られた時間発展モデルしか対応していないようなので、再サンプリングの部分は自作してしまいました。
また時間発展した座標を取り直すのではなく初期値を再サンプリングする手法をとっています。