Python による常微分方程式の数値解法 / scipy.integrate.odeint

spODE.py

import scipy as sp
from scipy.integrate import odeint

class spODE:
    def __init__(self, funcs, t0, inits, h=1/2**5, numOfDiv=1):
        self.funcs    = lambda v, t: [func(t, *v) for func in funcs]
        self.t0       = t0
        self.inits    = inits
        self.numOfDiv = numOfDiv
        self.h        = h
        self.tList    = sp.linspace(self.t0, self.t0 + self.h, self.numOfDiv + 1)

    def update(self):
        self.inits  = odeint(self.funcs, self.inits, self.tList)[-1]
        self.t0     = self.tList[-1]
        self.tList += self.h
        return self

    def print(self):
        print(self.t0, self.inits)
        return self

########
# test #
########
if __name__ == "__main__":
    def xDot(t, x, y):
        return y
    def yDot(t, x, y):
        return t - x
    funcs = [xDot, yDot]

    x0       = 0.0
    y0       = 0.0
    inits    = [x0, y0]
    
    t0       = 0.0
    tMax     = 7.5

    h        = 1/2**1
    numOfDiv = 2**3
    
    sol = spODE(funcs, t0, inits, h, numOfDiv)
    while sol.t0 < tMax:
        sol.update().print()
        
    print("真値:", str(tMax - sp.sin(tMax)), str(1 - sp.cos(tMax)))
    input("Push any key to exit ......")