常微分方程式

Python で常微分方程式の数値解 (近似解) を求める (フェールベルク法)

Python3 で常微分方程式の数値解 (近似解) を求めてみる。古典的ルンゲクッタ法による解法はどこにでもあるのでここではフェールベルク法 (6 段 4 次 & 5 次) を用いた。

任意の常微分方程式に適用しやすくする手段として解法全体をクラス化した。計算手順については『パソコンで見る天体の動き』(長沢 工, 檜山 澄子, 地人書館, 1992/10/1) に依った。

計算手順を大まかに述べる。
1. 最初に設定した刻み幅でとりあえず 1 刻みだけ計算してみる。
2. 4 次の解と 5 次の解とが同時に求まる。
3. その両者の差が真値との誤差であると見なす。
4. 真値との誤差が指定許容誤差に収まっていた場合は、その時点で 1 刻み分の計算を終え、そのときの 5 次の解を次の 1 刻みの初期値として使う。
5. 真値との誤差が指定許容誤差からはみ出ていた場合は、たった今おこなった計算をなかったものとし、指定許容誤差を満足するように刻み幅を狭くして計算し直す (新しい刻み幅は下のコード内の __hNew() 函数で計算している。理窟はわからない)。そうして得られた 5 次の解を次の 1 刻みの初期値として使う (下に示したコードは、刻み幅を変える必要のある場合には、最初に設定した 1 刻みを内部的に分割している。だから見かけ上の刻み幅は、最初に設定した刻み幅のままである)。

例として次の微分方程式を解いてみる。
x'(t) = x(t)^2 - t^2 - 2t + 2

確認の手段として t = 0 から 1 までを刻み幅 1 で解く。許容誤差は 0.0001 とする。実行結果を見ると、最初に設定した 1 刻みが 4 分割されたことと、指定許容誤差を満たしていることとがわかる。

fehlberg.py
# Python による常微分方程式の数値解法 / Fehlberg 法
from math import log, floor

class Fehlberg:
    def __init__(self, funcs, t0, inits, h, tol=1e-3):
        self.funcs    = funcs
        self.t0       = t0
        self.inits    = inits
        self.h        = h
        self.numOfDiv = None
        self.dim      = len(funcs)
        self.tol      = tol
        self.f        = [[None for i in range(self.dim)] for i in range(6)]
        self.temp     = [[None for i in range(self.dim)] for i in range(5)]
        self.temp2    = [None for i in range(self.dim)]
        self.inits4   = [None for i in range(self.dim)]
    def __floorB(self, num):
        if 0 < num < 1:
            return 2**floor(log(num)/log(2))
        else:
            return floor(num)
    def __hNew(self, h, err, tol):
        if err > tol:
            return 0.9 * h * (tol * h/(h * err))**(1/4)
        else:
            return h
    def __maxOfErr(self, listA, listB):
        for i in range(self.dim):
            self.temp2[i] = abs(listA[i] - listB[i])
        return max(self.temp2)
    def __preCalc(self, numOfDiv):
        preInits = [*self.inits]
        preT0    = self.t0
        preH     = self.h/numOfDiv
        for i in range(numOfDiv):
            for j in range(self.dim):
                self.f[0][j] = self.funcs[j](preT0 , *preInits)
                self.temp[0][j] = preInits[j] + preH * (1/4) * self.f[0][j]
            for j in range(self.dim):
                self.f[1][j] = self.funcs[j](preT0 + preH * (1/4) , *self.temp[0])
                self.temp[1][j] = preInits[j] + preH * (1/32) * (3 * self.f[0][j] + 9 * self.f[1][j])
            for j in range(self.dim):
                self.f[2][j] = self.funcs[j](preT0 + preH * (3/8) , *self.temp[1])
                self.temp[2][j] = preInits[j] + preH * (1/2197) * (1932 * self.f[0][j] - 7200 * self.f[1][j] + 7296 * self.f[2][j])
            for j in range(self.dim):
                self.f[3][j] = self.funcs[j](preT0 + preH * (12/13), *(self.temp[2]))
                self.temp[3][j] = preInits[j] + preH * ((439/216) * self.f[0][j] - 8 * self.f[1][j] + (3680/513) * self.f[2][j] - (845/4104) * self.f[3][j])
            for j in range(self.dim):
                self.f[4][j] = self.funcs[j](preT0 + preH , *self.temp[3])
                self.temp[4][j] = preInits[j] + preH * ((-8/27) * self.f[0][j] + 2 * self.f[1][j] - (3544/2565) * self.f[2][j] + (1859/4104) * self.f[3][j] - (11/40) * self.f[4][j])
            for j in range(self.dim):
                self.f[5][j] = self.funcs[j](preT0 + preH * (1/2) , *self.temp[4])
            if numOfDiv == 1:
                for j in range(self.dim):
                    self.inits4[j] = preInits[j] + preH * ((25/216) * self.f[0][j] + (1408/2565) * self.f[2][j] + (2197/4104) * self.f[3][j] - (1/5) * self.f[4][j])
            for j in range(self.dim):
                preInits[j] = preInits[j] + preH * ((16/135) * self.f[0][j] + (6656/12825) * self.f[2][j] + (28561/56430) * self.f[3][j] - (9/50) * self.f[4][j] + (2/55) * self.f[5][j])
            preT0 = preT0 + preH
        return preT0, self.inits4, preInits
    def update(self):
        self.numOfDiv = 1
        t, a4, a5     = self.__preCalc(self.numOfDiv)
        err           = self.__maxOfErr(a4, a5)
        if err < self.tol:
            self.t0, self.inits = t, a5
        else:
            self.numOfDiv = int(self.h/self.__floorB(self.__hNew(self.h, err, self.tol)))
            t, a4, a5     = self.__preCalc(self.numOfDiv)
            self.t0, self.inits = t, a5
        return self
    def print(self):
        print(self.t0, self.inits, self.numOfDiv)
        return self


#############################################################
# 解くべき微分方程式を定義する。リストで括っておく。
# 聯立方程式にも適用可能である。
def xDot(t, x):  # x'(t) = x(t)^2 - t^2 - 2t + 2
    return x**2 - t**2 - 2 * t + 2
funcs = [xDot]

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

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

# 刻み幅を指定する。
h = 1

# 許容誤差を指定する。
tol = 0.0001

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

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

#############################################################
trueVal = tMax + 1 - 1/(tMax + 1)
print("真値: " + str(trueVal))
print("真値との差: " + str(sol.inits[0] - trueVal))
input("Push any key to exit ......")

実行結果 (t [x(t)] 自動内部分割数):
image.png

許容誤差 tol を 0.1 にした場合:
image.png

関連:
https://qiita.com/souuuuta/items/722c9238a60c101183a2
http://ti-nspire.hatenablog.com/