Python による常微分方程式の数値解法 / 補外法

extrapo.py

# Python による常微分方程式の数値解法 / 補外法
class Extrapo:
    def __init__(self, funcs, t0, inits, h, numOfDiv=1):
        self.funcs    = funcs
        self.t0       = t0
        self.inits    = inits
        self.h        = h / numOfDiv
        self.numOfDiv = numOfDiv

    def __midPoints(self, funcs, t0, inits, H, S=7):
        dim       = len(funcs)
        RANGE     = 2**S+2
        listRange = [None for i in range(RANGE)]
        temp      = [None for i in range(dim)]
        listS     = [None for i in range(S)]
        t         = listRange.copy()
        t[0]      = t0
        x         = [[] for i in range(dim)]
        X         = x.copy()
        
        for i in range(dim):
            x[i]    = listRange.copy()
            x[i][0] = inits[i]
            X[i]    = listS.copy()

        for s in range(S+1):
            N = 2**s
            hs = H/N
            t[1] = t[0] + hs

            for i in range(dim):
                x[i][1] = x[i][0] + hs * funcs[i](t[0], *inits)
            
            for n in range(0, N):
                t[n+2] = t[n] + 2 * hs
                
                for i in range(dim):
                    temp[i] = x[i][n+1]
                        
                for i in range(dim):
                    x[i][n+2] = x[i][n] + 2 * hs * funcs[i](t[n+1], *temp)
        
            for i in range(dim):
                X[i][s-1] = (1/4) * (x[i][N-1] + 2 * x[i][N] + x[i][N+1])
            
        return X

    def __neville(self, listOfLists):
        dim    = len(listOfLists)
        S      = len(listOfLists[0])
        output = [None for i in range(dim)]
        x      = [v for v in listOfLists] # x = [[xList],[yList],[List]]
        
        for i in range(dim):
            x[i].insert(0, None) # x = [[0,xList],[0,yList],[0,List]] # 添字をテキストと揃える手段として、使わない [0] 番目の要素を設けておく。
    
        table = [[[None for i in range(S+1)] for i in range(S+1)] for i in range(dim)] # 添字をテキストと揃える手段として、使わない [0] 行目、[0] 列目を設けておく。
    
        for i in range(1,S+1): # テーブルの [1] 列目に点列を並べる。
            for j in range(1,S+1):
                for k in range(dim):
                    table[k][i][j] = x[k][i]

        for i in range(2,S+1):
            for j in range(2,i+1):
                for k in range(dim):
                    table[k][i][j] = table[k][i][j-1]+(table[k][i][j-1]-table[k][i-1][j-1])/(4**(j-1)-1)
                    output[k] = table[k][S][S]

        return output

    def update(self):
        for i in range(self.numOfDiv):
            self.inits = self.__neville(self.__midPoints(self.funcs, self.t0, self.inits, self.h))
            self.t0    = 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 = Extrapo(funcs, t0, inits, h)

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