古典的ルンゲ＝クッタ法のPythonプログラムをWikiPediaの記述と同じ変数名で書いてみた

教育用の覚書
Python だと汎用の Runge-Kutta 法のルーチンが６行で書けてしまう*1。以下のコードではRunge-Kutta法ルーチン rk4() を書き換えることなく、1階, 2階, 3階の常微分方程式(ordinary differential equation(s), ODE)を解いている*2*3。関数 rk4() で使っている変数名はルンゲ＝クッタ法 - Wikipediaに揃えた。C言語*4でコーディング/教育する気力が一気に失われてしまったｗｗｗｗｗｗ
次のコンプリートコードは次の微分方程式を積分して、あらかじめ求めておいた解の関数の値と比較している*5：

1. class getExpt():
   [1階]dydt=y${\displaystyle \frac{dy}{dt}=y}$　(比較する解：y=et$y=e^t$);
2. class getSint():
   [2階]d2ydt2=−y${\displaystyle \frac{d^{2}y}{d{t}^2}=-y}$　(比較する解：(y,dydt)=(sint,cost)${\displaystyle (y,\frac{dy}{dt})=(\sin t,\cos t)}$);
   [1階連立ODE] ddt(yy0,dydty1)=(dydt,−y)f(t,(y0,y1))${\displaystyle \frac{d}{dt}(\underset{y_0}{\underbrace{y}},\underset{y_1}{\underbrace{\frac{dy}{dt}}})=\underset{f(t,(y_0,y_1))}{\underbrace{(\frac{dy}{dt},-y)}}}$
3. class getExptSint():
   [3階]d3ydt3+d2ydt2+dydt+y=0${\displaystyle \frac{d^{3}y}{d{t}^3}+\frac{d^{2}y}{d{t}^2}+\frac{dy}{dt}+y=0}$
   (比較する解：(y,dydt,d2ydt2)=(e−t+sint,−e−t+cost,e−t−sint)${\displaystyle (y,\frac{dy}{dt},\frac{d^{2}y}{d{t}^2})=(e^{-t}+\sin t,-e^{-t}+\cos t,e^{-t}-\sin t)}$);
   [1階連立ODE] ddt(yy0,dydty1,d2ydt2y2)=(dydt,d2ydt2,−y−dydt−d2ydt2)f(t,(y0,y1,y2))${\displaystyle \frac{d}{dt}(\underset{y_0}{\underbrace{y}},\underset{y_1}{\underbrace{\frac{dy}{dt}}},\underset{y_2}{\underbrace{\frac{d^{2}y}{d{t}^2}}})=\underset{f(t,(y_0,y_1,y_2))}{\underbrace{(\frac{dy}{dt},\frac{d^{2}y}{d{t}^2},-y-\frac{dy}{dt}-\frac{d^{2}y}{d{t}^2})}}}$
4. class spring():
   [2階]単振動(調和振動子)の問題 ma=−kx${\displaystyle ma=-kx}$, すなわち
   位置：x(t)${\displaystyle x(t)}$, 速度：v(t)=dxdt${\displaystyle v(t)=\frac{dx}{dt}}$, 加速度：a(t)=d2xdt2${\displaystyle a(t)=\frac{d^{2}x}{d{t}^2}}$ よりd2xdt2=−kmx${\displaystyle \frac{d^{2}x}{d{t}^2}=-\frac{k}{m}x}$
   (比較する解：(x(t),v(t))=(sinωt,ωcosωt)$(x(t),v(t))=(\sin \omega t,\omega \cos \omega t)$, ω=km−−−√${\displaystyle \omega =\sqrt{\frac{k}{m}}}$);
   [1階連立ODE] ddt(xy0,vy1)=(v,−kmy)f(t,(y0,y1))${\displaystyle \frac{d}{dt}(\underset{y_0}{\underbrace{x}},\underset{y_1}{\underbrace{v}})=\underset{f(t,(y_0,y_1))}{\underbrace{(v,-\frac{k}{m}y)}}}$

初期条件は比較する解に t=0 を代入して作っている。

import numpy
from pylab import *

# classical 4-th order Runge-Kutta method subroutine
# variables letters are the same as those in
# https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods
def rk4 ( y, t, h, f ):
    k1 = h * f( t, y )
    k2 = h * f( t + 0.5*h, y + 0.5*k1 )
    k3 = h * f( t + 0.5*h, y + 0.5*k2 )
    k4 = h * f( t + h, y + k3 )
    return y + ( k1 + 2*k2 + 2*k3 + k4 )/6, t + h

class getExpt():
    def __init__ (self):
        print(self.__class__,' y(t)=exp(t)')
    def ref ( self, t ):
        return exp(t)
    def rhsODE ( self, t, y ):
        return y

class getSint():
    def __init__ (self):
        print(self.__class__,' y(t)=sin(t)')
    def ref ( self, t ):
        return numpy.array([sin(t),cos(t)])
    def rhsODE ( self, t, y ):
        return numpy.array([y[1],-y[0]])

class getExptSint():
    def __init__ (self):
        print(self.__class__,' y(t)=exp(-t)+sin(t)')
    def ref ( self, t ):
        et=exp(-t); st=sin(t);
        return numpy.array([et+st,-et+cos(t),et-st])
    def rhsODE ( self, t, y ):
        return numpy.array([y[1],y[2],-y[0]-y[1]-y[2]])

class spring():
    def __init__ (self,mass,springConstant):
        self.mass= mass
        self.springConstant= springConstant
        m= mass
        k= springConstant
        self.angularFrequency= (k/m)**0.5
        print(self.__class__,' linear spring motion m=',m,' k=',k)
    def ref ( self, t ):
        w= self.angularFrequency
        x= sin(w*t)
        v= w*cos(w*t)
        return numpy.array([x,v])
    def rhsODE ( self, t, y ):
        k= self.springConstant
        m= self.mass
        x= y[0]
        v= y[1]
        a= -(k/m)*x
        return numpy.array([v,a])

if __name__=="__main__":

    # choose class to e[x]ecute
    #x= getExpt();
    #x= getSint();
    #x= getExptSint();
    x= spring(mass=1,springConstant=1);

    # set [ref]erence [sol]ution function
    #   and [f]unction appears [r]ight [h]and [s]ide of
    #   [O]rdinary [D]ifferential [E]quation
    sol= x.ref; f= x.rhsODE;

    # set step-size and initial condition(s)
    h = 0.01; t = 0; y = sol(t); print(t,y,sol(t),y-sol(t))

    # set terminal condition(s)
    tEnd = 0.1;

    # integration and comparison to reference solution
    while t < tEnd:
        y, t = rk4( y, t, h, f )
        print(t,y,sol(t),y-sol(t))

*1:ルンゲクッタ python - Google 検索で検索すると上位に上がるページに意外と汎用性のあるコードが無い。

*2:ルンゲクッタ法は1階の連立常微分方程式を解く計算スキームなので、n階の微分方程式はn個の関数 y=(y,y′,...,y(n−1))$\mathit{y}=(y,y^{\prime },...,y^{(n-1)})$ に対する連立方程式に書き換える。

*3:Pythonコードなので関数 rk4() は任意の配列で計算ができる。例えば3次元流体をスペクトル法で解くとか（配列のメモリアロケーションどうなっとるか知らんけどｗｗｗｗｗｗ）。

*4:C言語だと構造体へのポインタ
　str *rk4( str *y, double t, double h, (str*)f(double,*str) ) { str *k1, *k2, *k3, *k4; ...; return k?; }
で宣言すればいいのかな？

*5:上掲のプログラムは[4]の調和振動子を解いている。[1][2][3]の方程式を解くためには x= の行を適宜、書き換えればよい。