Rで微分方程式

  • まずは離散的な現象を…
    • 時刻t=0に1個の細菌が、単位時間1ごとに2倍になるという
    • x(t+1) = 2 x(t);x(0) = 1
max.t <- 20
xs <- rep(0,max.t+1)
xs[1] <- 1
for(i in 1:max.t){
	xs[i+1] <- 2 * xs[i]
}
plot(0:max.t,xs)
    • 差分を取って同じことをする
    • \Delta x(t) = x(t+1)-x(t) = 2x(t)-x(t)=x(t)
xs.2 <- rep(0,max.t+1)
xs.2[1] <- 1
for(i in 1:max.t){
	delta <- xs.2[i]
	xs.2[i+1] <- xs.2[i] + delta
}

plot(0:max.t,xs.2)
  • 練習
    • 単位時間ごとに k(k=3,0.9) 倍になる例をやってみよう
  • 単位時間 \Delta t = 1ごとにすべての細菌が2倍になるのではなく、単位時間\Delta t = \frac{1}{100}ごとに、全体の3%だけが2倍になるとのこと
    • x(t+0.01) = x(t) + 0.03 x(t)
    • 差分で考えると
    • x(t+0.01) - x(t) = 0.03 x(t)
    • Rでやってみる
unit.time <- 0.01
fraction.growth <- 0.03

times <- seq(from=0,to=10,by=0.01)

xs.3 <- rep(0,length(times))
xs.3[1] <- 1
for(i in 1:(length(xs.3)-1)){
	delta <- xs.3[i] * fraction.growth
	xs.3[i+1] <- xs.3[i] + delta
}
plot(times,xs.3)
    • プロットすると、「ぐんぐん増えて」「指数関数的に」増えている
    • 本当に指数関数的?
plot(times,log(xs.3),type="l")
abline(0,fraction.growth/delta.t,col=2)
  • どうして何となくうまく行っているか…
    • x(t) = e^{at}のとき
    • \frac{dx(t)}{dt} = a e^{at} = a x(t)
  • 単位時間を1、0.01にしたけれど、0.1や0.05だと…
unit.time <- 0.1
fraction.growth <- 3 * unit.time

times <- seq(from=0,to=10,by=0.01)

xs.3 <- rep(0,length(times))
xs.3[1] <- 1
for(i in 1:(length(xs.3)-1)){
	delta <- xs.3[i] * fraction.growth
	xs.3[i+1] <- xs.3[i] + delta
}
plot(times,xs.3)

plot(times,log(xs.3),type="l")
abline(0,fraction.growth/delta.t,col=2)
    • すごく細かく!
unit.time <- 0.0001
fraction.growth <- 3 * unit.time

times <- seq(from=0,to=10,by=0.01)

xs.3 <- rep(0,length(times))
xs.3[1] <- 1
for(i in 1:(length(xs.3)-1)){
	delta <- xs.3[i] * fraction.growth
	xs.3[i+1] <- xs.3[i] + delta
}
plot(times,xs.3)

plot(times,log(xs.3),type="l")
abline(0,fraction.growth/delta.t,col=2)
  • 速度が座標で決まるようなとき
    • \frac{d x}{dt} = function(x)
    • たとえばfunction(x) = 0.001 * \cos{x}
max.t <- 10000
times <- 0:max.t
xs <- rep(0,length(times))
xs[1] <- 4

speed.fx <- function(x){
	return(sin(x)*0.001)
}


for(i in 1:(length(times)-1)){
	delta <- speed.fx(xs[i])
	xs[i+1] <- xs[i] + delta
}

plot(times,xs,type="l")
    • 初期位置が変わると…
max.t <- 10000
times <- 0:max.t
xs <- rep(0,length(times))
xs[1] <- 1

speed.fx <- function(x){
	return(sin(x)*0.001)
}


for(i in 1:(length(times)-1)){
	delta <- speed.fx(xs[i])
	xs[i+1] <- xs[i] + delta
}

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