ロトカ・ボルテラの方程式

ロトカ・ボルテラの方程式
倉橋貴彦
ロトカ・ボルテラの方程式（捕食者・非捕食者の関係を表した連立微分方程式）
dx
 x(  y )  f x, y 
dt
dy
  y (  x)  g x, y 
dt
X:捕食者
Y：非捕食者

k x1  tf x n , y n


 txn   y n
修正オイラー法（２次のルンゲ・クッタ法）
による時間進展式
1
k x1  k x 2 
2
1
y n 1  y n  k y1  k y 2
2
x n 1  x n 




k x 2  tf x n 1, y n 1


 txn 1   y n 1






 t x n  tf x n , y n    y n  tg x n , y n

k y1  tg x n , y n


  tyn (  x n )

k y 2  tg x n 1, y n 1


  tyn 1   x n 1



 

  t y n  tg x n , y n    x n  tf x n , y n

Start
フローチャート例
n=0，Δt ，x0=x|t=0 ，y0=y|t=0 ，T
(※T/Δtが整数になるようにΔtを選択する。）
imax = T/Δt
T:全時間
n：ステップ数
imax：最大ステップ数
n*Δt, x0, y0

k x1  tf x n , y n


 txn   y n
n=n+1


n 1
n 1
k x 2  tf x n 1, y n 1
  y 
 t x  tf x , y    y
 tx
kx1, kx2, ky1, ky2の計算
n
1
k x1  k x 2 
2
1
y n 1  y n  k y1  k y 2
2
n

x n 1  x n 


k y1  tg x n , y n


 tg x n , y n

  tyn (  x n )


  x 
  t y  tg x , y    x
yn+1
  ty
n 1
n 1
n
n
Yes
No
n

k y 2  tg x n 1, y n 1
n*Δt, xn+1,
n
n=imax
end
n
n

 tf x n , y n

Fortranによるプログラム例
implicit double precision (a-h,o-z)
c
open ( 10,file='output.dat')
c
x0 = 0.3d0
y0 = 0.7d0
delt = 0.2d0
imax = 500
c
alpha = 0.2d0
beta = 0.4d0
gamma = 0.3d0
delta = 0.3d0
c
c
write(10,*) 0.d0, x0, y0
c
do istep = 1,imax
c
c----- value of right hand side term at n step
c
call fx
& ( f0, alpha, beta, x0, y0 )
call fy
& ( g0, gamma, delta, x0, y0 )
c
c----- value update by Euler method
c
ax1 = x0 + delt * f0
ay1 = y0 + delt * g0
c
c----- value of right hand side term at n+1 step
c
call fx
& ( f1, alpha, beta, ax1, ay1 )
c
call fy
& ( g1, gamma, delta, ax1, ay1 )
c
c----c
akx1 = delt * f0
akx2 = delt * f1
c
aky1 = delt * g0
aky2 = delt * g1
c
c----- value update by modified Euler method
c
x1 = x0 + 0.5d0 * ( akx1 + akx2 )
y1 = y0 + 0.5d0 * ( aky1 + aky2 )
c
write(10,*) delt*istep, x1, y1
c
c----- step change
c
x0 = x1
y0 = y1
c
end do
c
stop
end
c
c==============================================
c
subroutine fx
& ( f0, alpha, beta, x0, y0 )
c
c==============================================
c
implicit double precision (a-h,o-z)
c
f0 = x0 * ( alpha - beta * y0 )
c
return
end
c
c==============================================
c
subroutine fy
& ( g0, gamma, delta, x0, y0 )
c
c==============================================
c
implicit double precision (a-h,o-z)
c
g0 = - y0 * ( gamma - delta * x0 )
c
return
end
計算条件・計算結果
計算条件
x0 = 0.3
y0 = 0.7
Δt= 0.2
time step = 500
α = 0.2
β = 0.4
γ = 0.3
δ = 0.3
x,yの値の時系列変化

Download Report