2014-09-14
ロトカ=ヴォルテラの向きの分離

my.Moving.Frame.2 <- function(X,k=length(X[1,])){
	n <- length(X[,1]) # No. points
	d <- length(X[1,]) # Dimension
	# inter-point distance
	L <- rep(0,n-1)
	for(i in 1:(n-1)){
		L[i] <- sqrt(sum((X[i,]-X[i+1,])^2))
	}
	diff.X <- list()
	#diff.X[[1]] <- X
	diff.X[[1]] <- apply(X,2,diff)
	diff.L <- list()
	diff.L[[1]] <- L
	for(i in 2:k){
		diff.X[[i]] <- apply(diff.X[[i-1]],2,diff)
		tmp.n.pt <- length(diff.X[[i]][,1])
		tmp.L <- rep(0,tmp.n.pt-1)
		for(j in 1:(tmp.n.pt-1)){
			tmp.L[j] <- sqrt(sum((diff.X[[i]][j,]-diff.X[[i]][j+1,])^2))
		}
		diff.L[[i]] <- tmp.L
	}
	mf <- list() # Moving Frame
	mf[[2]] <- diff.X[[2]]/sqrt(apply(diff.X[[2]]^2,1,sum))
	if(k>2){
		for(i in 3:k){
			n.pt <- length(diff.X[[i]][,1])
			ip <- apply(diff.X[[i]] * mf[[i-1]][1:n.pt,],1,sum)
			mf[[i]] <- diff.X[[i]] - ip * mf[[i-1]][1:n.pt,]
			mf[[i]] <- mf[[i]]/sqrt(apply(mf[[i]]^2,1,sum))
		}
	}
	n.pt <- length(diff.X[[k]][,1])
	mf[[1]] <- matrix(0,n.pt,d)
	for(i in 1:n.pt){
		tmp.mat <- matrix(0,d-1,d-1)
		tmp.vec <- rep(0,d-1)
		for(j in 2:d){
			tmp <- mf[[j]][i,]
			tmp.mat[j-1,] <- tmp[1:(d-1)]
			tmp.vec[j-1] <- tmp[d]
		}
		mf[[1]][i,] <- c(-solve(tmp.mat) %*% tmp.vec,1)
	}
	mf[[1]] <- mf[[1]]/sqrt(apply(mf[[1]]^2,1,sum))
	
	Cs <- matrix(0,d-1,length(diff.L[[1]]))
	for(i in 1:(d-1)){
		n.pt <- length(diff.L[[i+1]])
		Cs[i,1:n.pt] <- diff.L[[i+1]][1:n.pt]/L[1:n.pt]
	}
	mf.array <- diff.array <- array(0,c(n,k,k))
	for(i in 1:k){
		tmp.n <- length(mf[[i]][,1])
		delta.n <- n-tmp.n
		init.n <- ceiling(delta.n/2)
		if(init.n==0)init.n <- 1
		r.id <- init.n:(init.n+tmp.n-1)
		mf.array[r.id,i,] <- mf[[i]]
		
		tmp.n <- length(diff.X[[i]][,1])
		delta.n <- n-tmp.n
		init.n <- ceiling(delta.n/2)
		if(init.n==0)init.n <- 1
		r.id <- init.n:(init.n+tmp.n-1)
		diff.array[r.id,i,] <- diff.X[[i]]
	}
	
	return(list(L=L,mf=mf,mf.array=mf.array,diff.X=diff.X,diff.array=diff.array,diff.L=diff.L,Cs=Cs))
}
library(deSolve)
# 微分方程式
# この関数my.ODEに 引数 t は使われていないが
# runge-kuttaの本体関数 rk() との関係で、この引数は必要(らしい)
# 時刻を表す引数
# xは対象となっている変数の値のベクトル
# parmsは微分方程式の係数
# 差分を計算して返す
# この連立微分方程式はロトカ=ヴォルテラ
# parms[1],parms[4]は負、parms[2],parms[3]は正
my.ODE <- function(t, x, parms)  {

	dX <- (parms[1] + parms[2]*x[2])*x[1]
	dY <- (parms[3] + parms[4]*x[1])*x[2]

	res <- c(dX, dY)
	list(res)

}


## The parameters 
# 連立微分方程式の４係数を決める
# ２タイプあって、回転方向が逆
library(MCMCpack)
n <- 2
A1 <- rdirichlet(1,rep(1,4))*5
A1[c(1,4)] <- -A1[c(1,4)]
#A2 <- rdirichlet(1,rep(1,4))*5
#A2[c(1,4)] <- -A2[c(1,4)]
A2 <- -A1 + rnorm(4)*0.01
## vector of timesteps
# 計算する時刻
times <- seq(0, 10, length = 1001)

## 初期値
# 固定点の周囲から始まるとうまく周期性を出すので
# 固定点に若干のずれを入れる

N <- 100
X <- list()
dets <- matrix(0,length(times),N)
types <- c()
for(i in 1:N){


## Euler method
# hini 計算するための時間差分
if(runif(1)<0.5){
	xstart <- c(-A1[3]/A1[4], -A1[1]/A1[2]) + rnorm(2)*0.1

out1  <- rk(xstart, times, my.ODE, A1, hini = 0.1, method = "rk4")
types <- c(types,1)
}else{
	xstart <- c(-A2[3]/A2[4], -A2[1]/A2[2]) + rnorm(2)*0.1

out1  <- rk(xstart, times, my.ODE, A2, hini = 0.1, method = "rk4")
types <- c(types,2)
}
X[[i]] <- out1[,2:(n+1)]
mf1 <- my.Moving.Frame.2(X[[i]])
dets1 <- c()
for(ii in 1:length(mf1$diff.array[,1,1])){
	dets1 <- c(dets1,det(mf1$diff.array[ii,,]))
}

dets[,i] <- dets1
}
a <- which(types==1)[1]
b <- which(types==2)[1]
plot(X[[1]],type="l",xlim=range(c(X[[a]][,1],X[[b]][,1])),ylim=range(c(X[[a]][,2],X[[b]][,2])),col=types[1])
for(i in 2:N){
	lines(X[[i]],col=types[i])
}
matplot(dets,type="l",col=types)