Rで微分方程式2

• ２要素が作る微分方程式
• 

p <- 1
q <- 1
r <- 1
s <- 1

n.time <- 1000
xs <- ys <- rep(0,n.time)
xs[1] <- 2
ys[1] <- 1

delta.t <- 0.01

for(i in 2:n.time){
	dx <- p * xs[i-1] -r *xs[i-1]*ys[i-1]
	dy <- s * xs[i-1] * ys[i-1] - s * ys[i-1]
	
	xs[i] <- xs[i-1] + dx*delta.t
	ys[i] <- ys[i-1] + dy*delta.t
}

matplot(cbind(xs,ys),type="l")

plot(xs,ys,type="l")

• ロトカ=ヴォルテラ