Python / 常微分方程式の数値解法 / 補外法 / 應用例 / 単振り子 / 時間領域のグラフを描く

関連: ログイン - はてな

from myODE.Extrapo import *               # myODE.Extrapo は補外法で 1 階常微分方程式を解くためのプログラム
from math          import sin    as Sin
from math          import pi     as Pi
from matplotlib    import pyplot as plt

class Pend:
    def __init__(self, length, radian, velocity=0, h=1/2**5, k=0):
        g = 9.80665

        def dotRad(t, rad, vel):  # θ'(t) = velocity(t)
            return vel
        def dotVel(t, rad, vel):  # velocity'(t) = -(9.8/length) * sin(θ)
            return (-g/length) * Sin(rad) - k * vel                    # k は速度に比例する減衰係数
        
        self.sol = Extrapo([dotRad, dotVel], 0, [radian, velocity], h) # これが初期値を 1 刻みだけ更新する函数。
        
        self.__convData()
        
    def update(self):
        self.sol.update()
        self.__convData()
        return self
    
    def print(self):
        print(self.t0, self.rad)
        return self

    def __convData(self):
        self.t0  = self.sol.t0
        self.rad = self.sol.inits[0]


########
# test #    
########
if __name__ == "__main__":
    pend = Pend(length = 1, radian = Pi*0.999) # 振り子を実体化して、
    tList   = [pend.t0 ]
    radList = [pend.rad]
    while pend.t0 < 15:
        pend.update().print()                  # 振り子の状態をどんどん更新して、
        tList  .append(pend.t0 )               # データを溜める。
        radList.append(pend.rad)
        
    plt.plot(tList, radList)
    plt.xlabel("Second")
    plt.ylabel("Radian")
    plt.show()

実行結果:

(支点のほぼ真上から振っているので等時性の破れがはっきりと現れている)