Runge-Kutta-Gill法のPythonコード - あらきけいすけのメモ帳

4段4次のルンゲ＝クッタ法で記憶領域の量が少ないRunge-Kutta-Gill法*1のコードをPythonで書いてみたのでメモ*2。関数rk4()のインターフェースは古典的ルンゲ＝クッタ法のPythonプログラムをWikiPediaの記述と同じ変数名で書いてみた - あらきけいすけのメモ帳に合わせてある。解いている微分方程式は $\ddot{y}=-y$ , すなわち $\frac{d}{dt}\left(y(t),\dot{y}(t)\right)=\left(\dot{y}(t),-y(t)\right)$ , 初期条件は $\left(y(0),\dot{y}(0)\right)=\left(\cos0,-\sin0\right)=(1,0)$ , 解は $\left(y(t),\dot{y}(t)\right)=\left(\cos t,-\sin t\right)$ .

import numpy

def rk4 ( y, t, h, f ):
  # Gill's low strage method is implemented.
  p= numpy.zeros_like(y); q= numpy.zeros_like(y)
  p=f(t,y); p=p*h*(1./2.); rt2=2**-.5; y+=p; q=p*2
  p=f(t,y); p=p*h*(1-rt2) - q*(1-rt2); y+=p; q=p*2 + q*rt2
  p=f(t,y); p=p*h*(1+rt2) - q*(1+rt2); y+=p; q=p*2 - q*rt2
  p=f(t,y); p=p*h*(1./6.) - q*(2./6.); y+=p;
  return y, t+h

def RHS(t,y):
  return numpy.array([y[1],-y[0]])

y=numpy.array([1.0,0.0]); t=0; h=0.001
for i in range(100):
  y, t =rk4( y, t, h, RHS )
  print(y[0]-numpy.cos(t))

*1:ストレージは3個（このコードでは y, p, q）で、典型的なRunge-Kutta法より係数がちょっと複雑。[参考]ルンゲ・クッタ・ギル法
4次の公式の最小ストレージは2個らしい: David I.Ketcheson, "Runge–Kutta methods with minimum storage implementations", J. Comput. Phys. Vol. 229, Issue 5, 1 March 2010, Pages 1763-1773 https://doi.org/10.1016/j.jcp.2009.11.006 http://www.davidketcheson.info/assets/papers/2010_LSRK_postprint.pdf.

*2:意外にPythonでのチートシート的なページが無い。fortranなどの他言語に移植することを少しだけ意識した書き方にしたつもり。