python - 振り子の統合。オーバーフロー。パイソン

Question

4次のルンゲ・クッタ積分を使用して振り子の速度と力を計算するコードを書きましたが、残念ながら次のエラーが発生するため実行できません。

 30: RuntimeWarning: overflow encountered in double_scalars
   th2 = th[j+1] + (ts/2)*om[j+1]
 33: RuntimeWarning: overflow encountered in double_scalars
   om2 = om[j+1] + (ts/2)*f2
 39: RuntimeWarning: invalid value encountered in double_scalars
   th3 = th2 + (ts/2)*om2
 40: RuntimeWarning: invalid value encountered in sin
   f3 = (-q*om2 - sin(th2) + b*cos(od)*ts)

何が欠けているのか、何が間違いなのかわかりません。どんなアイデアでも役に立ちます。ありがとう！

これは私のコードです:

 from scipy import *
 from matplotlib.pyplot import *
 ##A pendulum simulation using fourth order 
 ##Runge-Kutta differentiation

 ts=.05 #time step size, dt_a
 td=20 #trial duration, dt_b
 te=int(td/ts) #no of timesteps

 q=0.5 #scaled damping coefficient or friction factor
 m=1 #mass
 g=9.81 #grav. acceleration
 l=9.81 #length of pendulum

 th=[((rand()*2)-1)*pi] #random initial angle, [-pi,pi]
 om=[0] #initial angular velocity
 b=0.9 #scaled initial (force)/ml or driving coef.
 od=0.66 #driving force freq


 for j in range(te):
     #Euler approximation
     th.append(th[j] + ts*om[j]) #storing theta values
     f1 = (-q*om[j] - sin(th[j]) + b*cos(od)*ts)
     #ts += ts 
     om.append(om[j] + ts*f1)    

     #ts=.05
     #1 at mid-interval
     th2 = th[j+1] + (ts/2)*om[j+1]

     f2 = (-q*om[j+1] - sin(th[j+1]) + b*cos(od)*ts)
     om2 = om[j+1] + (ts/2)*f2
     #ts += ts

     #ts=.05
     #2 at mid-interval

     th3 = th2 + (ts/2)*om2
     f3 = (-q*om2 - sin(th2) + b*cos(od)*ts)
     om3 = om2 + (ts/2)*f3
     #ts += ts


     #next time step

     th4 = th3 + (ts)*om3
     f4 = (-q*om3 - sin(th3) + b*cos(od)*ts)
     om4 = om3 + (ts)*f4

     ts += ts


     dth=(om[j] + 2*om[j+1] + 2*om2 + om3)/6
     dom=(f1 + 2*f2 + 2*f3 + f4)/6
     th[j+1] = th[j] + ts*dth #integral of angle 
     om[j+1] = om[j] + ts*dom #integral of velocity

 plot(th,om),xlabel('Angle'),ylabel('')
 show()