あらきけいすけのメモ帳

あらきけいすけの雑記帳2

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

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

  1. class getExpt():
    [1階]\displaystyle\frac{dy}{dt}=y (比較する解:y=e^t);
  2. class getSint():
    [2階]\displaystyle\frac{d^2y}{dt^2}=-y (比較する解:\displaystyle\left(y,\frac{dy}{dt}\right)=(\sin t,\cos t));
    [1階連立ODE] \displaystyle\frac{d}{dt}\Big(\underbrace{y}_{y_0},\underbrace{\frac{dy}{dt}}_{y_1}\Big)=\underbrace{\left(\frac{dy}{dt},-y\right)}_{f(t,(y_0,y_1))}
  3. class getExptSint():
    [3階]\displaystyle\frac{d^3y}{dt^3}+\frac{d^2y}{dt^2}+\frac{dy}{dt}+y=0
    (比較する解:\displaystyle\left(y,\frac{dy}{dt},\frac{d^2y}{dt^2}\right)=(e^{-t}+\sin t,-e^{-t}+\cos t,e^{-t}-\sin t));
    [1階連立ODE] \displaystyle\frac{d}{dt}\bigg(\underbrace{y}_{y_0},\underbrace{\frac{dy}{dt}}_{y_1},\underbrace{\frac{d^2y}{dt^2}}_{y_2}\bigg)=\underbrace{\left(\frac{dy}{dt},\frac{d^2y}{dt^2},-y-\frac{dy}{dt}-\frac{d^2y}{dt^2}\right)}_{f(t,(y_0,y_1,y_2))}
  4. class spring():
    [2階]単振動(調和振動子)の問題 \displaystyle m a=- k x, すなわち
    位置:\displaystyle x(t), 速度:\displaystyle v(t)=\frac{dx}{dt}, 加速度:\displaystyle a(t)=\frac{d^2x}{dt^2} より\displaystyle \frac{d^2x}{dt^2}=-\frac{k}{m}x
    (比較する解:(x(t),v(t))=(\sin\omega t,\omega\cos\omega t), \displaystyle\omega=\sqrt{\frac{k}{m}});
    [1階連立ODE] \displaystyle\frac{d}{dt}\Big(\underbrace{x}_{y_0},\underbrace{v}_{y_1}\Big)=\underbrace{\left(v,-\frac{k}{m}x\right)}_{f(t,(y_0,y_1))}

初期条件は比較する解に 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

# each of the following four classes has
# __init__(): constructor
# ref(t): reference solution
# rhsODE(t,y): f(t,y) of the 1st order ODE dy/fx=f(t,y)
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=1,springConstant=1):
        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=2,springConstant=3);

    # 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個の関数 \boldsymbol{y}=(y,y',...,y^{(n-1)}) に対する連立方程式に書き換える。

*3:Pythonなので関数 rk4() には識別子さえ渡せばいいので、任意サイズの配列 y で計算ができる。例えば3次元流体をスペクトル法で解くとか(配列のメモリアロケーションどうなっとるか知らんけどwwwwww)。

*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= の行を適宜、書き換えればよい。