Python による常微分方程式の数値解法 / Shanks による 12 段 8 次の Runge-Kutta 法

Shanks8.py

# Python による常微分方程式の数値解法 / Shanks による 12 段 8 次の Runge-Kutta 法
class Shanks8:
    def __init__(self, funcs, t0, inits, h, numOfDiv=1):
        self.funcs    = funcs
        self.t0       = t0
        self.inits    = inits
        self.dim      = len(funcs)
        self.numOfDiv = numOfDiv
        self.h        = h / self.numOfDiv
        self.f        = [[None for i in range(self.dim)] for i in range(12)]
        self.temp     = [[None for i in range(self.dim)] for i in range(11)]
    def update(self):
        for i in range(self.numOfDiv):
            for j in range(self.dim):
                self.f[0][j] = self.funcs[j](self.t0 , *self.inits)
                self.temp[0][j] = self.inits[j] + (self.h/9) * (self.f[0][j])
            for j in range(self.dim):
                self.f[1][j] = self.funcs[j](self.t0 + self.h/9 , *self.temp[0])
                self.temp[1][j] = self.inits[j] + (self.h/24) * (self.f[0][j] + 3 * self.f[1][j])
            for j in range(self.dim):
                self.f[2][j] = self.funcs[j](self.t0 + self.h/6 , *self.temp[1])
                self.temp[2][j] = self.inits[j] + (self.h/16) * (self.f[0][j] + 3 * self.f[2][j])
            for j in range(self.dim):
                self.f[3][j] = self.funcs[j](self.t0 + self.h/4 , *self.temp[2])
                self.temp[3][j] = self.inits[j] + (self.h/500) * (29 * self.f[0][j] + 33 * self.f[2][j] - 12 * self.f[3][j])
            for j in range(self.dim):
                self.f[4][j] = self.funcs[j](self.t0 + self.h/10 , *self.temp[3])
                self.temp[4][j] = self.inits[j] + (self.h/972) * (33 * self.f[0][j] + 4 * self.f[3][j] + 125 * self.f[4][j])
            for j in range(self.dim):
                self.f[5][j] = self.funcs[j](self.t0 + self.h/6 , *self.temp[4])
                self.temp[5][j] = self.inits[j] + (self.h/36) * (-21 * self.f[0][j] + 76 * self.f[3][j] + 125 * self.f[4][j] - 162 * self.f[5][j])
            for j in range(self.dim):
                self.f[6][j] = self.funcs[j](self.t0 + self.h/2 , *self.temp[5])
                self.temp[6][j] = self.inits[j] + (self.h/243) * (-30 * self.f[0][j] - 32 * self.f[3][j] + 125 * self.f[4][j] + 99 * self.f[6][j])
            for j in range(self.dim):
                self.f[7][j] = self.funcs[j](self.t0 + self.h*2/3 , *self.temp[6])
                self.temp[7][j] = self.inits[j] + (self.h/324) * (1175 * self.f[0][j] - 3456 * self.f[3][j] - 6250 * self.f[4][j] + 8424 * self.f[5][j] + 242 * self.f[6][j] - 27 * self.f[7][j])
            for j in range(self.dim):
                self.f[8][j] = self.funcs[j](self.t0 + self.h/3 , *self.temp[7])
                self.temp[8][j] = self.inits[j] + (self.h/324) * (293 * self.f[0][j] - 852 * self.f[3][j] - 1375 * self.f[4][j] + 1836 * self.f[5][j] - 118 * self.f[6][j] + 162 * self.f[7][j] + 324 * self.f[8][j])
            for j in range(self.dim):
                self.f[9][j] = self.funcs[j](self.t0 + self.h*5/6 , *self.temp[8])
                self.temp[9][j] = self.inits[j] + (self.h/1620) * (1303 * self.f[0][j] - 4260 * self.f[3][j] - 6875 * self.f[4][j] + 9990 * self.f[5][j] + 1030 * self.f[6][j] + 162 * self.f[9][j])
            for j in range(self.dim):
                self.f[10][j] = self.funcs[j](self.t0 + self.h*5/6 , *self.temp[9])
                self.temp[10][j] = self.inits[j] + (self.h/4428) * (-8595 * self.f[0][j] + 30720 * self.f[3][j] + 48750 * self.f[4][j] - 66096 * self.f[5][j] + 378 * self.f[6][j] - 729 * self.f[7][j] - 1944 * self.f[8][j] - 1296 * self.f[9][j] + 3240 * self.f[10][j])
            for j in range(self.dim):
                self.f[11][j] = self.funcs[j](self.t0 + self.h , *self.temp[10])
            for j in range(self.dim):
                self.inits[j] += (self.h/840) * (41 * self.f[0][j] + 216 * self.f[5][j] + 272 * self.f[6][j] + 27 * self.f[7][j] + 27 * self.f[8][j] + 36 * self.f[9][j] + 180 * self.f[10][j] + 41 * self.f[11][j])
            self.t0 += self.h
        return self
    def print(self):
        print(self.t0, self.inits)
        return self

########
# test #
########
if __name__ == "__main__":
    # 解くべき聯立微分方程式を定義する。リストで括っておく。
    def xDot(t, x, y):  # x'(t) = y(t)
        return y
    def yDot(t, x, y):  # y'(t) = t - x(t)
        return t - x
    funcs = [xDot, yDot]

    # 独立変数の開始値と終了値とを指定する。
    t0 = 0
    tMax = 7.5

    # 従属変数の初期値を指定する。リストで括っておく。
    x0 = 0
    y0 = 0
    inits = [x0, y0]

    # 刻み幅を指定する。2 の整数乗分の 1 にすることが望ましい。
    h = 1/2**2

    # 1 刻みの内部分割数を指定する場合は 2 の整数乗にすることが望ましい。
    # numOfDiv = 2**4

    # 1 刻みだけ計算する函数を実体化して、
    sol = Shanks8(funcs, t0, inits, h)

    # 初期値を更新しながら必要な回数だけ実行を繰り返す。
    while sol.t0 < tMax:
        sol.update().print()