2

I wrote a code for calculating the velocity and force of a pendulum using 4th order Runge-Kutta integration, unfortunately I cannot run it since I get this error:

 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)

I am not sure what I'm missing or what is my mistake. Any idea would be helpful. Thanks!

This is my code:

 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()
Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
Estefy
  • 424
  • 2
  • 6
  • 20
  • 1
    Well first copy paste seems to run fine for me... small note: things are approximately 2 zillion times easier to debug when variable names are less one-lettery. M makes a lot of sense when you're writing it 100 times in a physics equation. Less sense when you're trying to remember what it stood for 1000 lines later – en_Knight Mar 04 '16 at 20:57
  • I think you can find an answer to your problem here: http://stackoverflow.com/questions/7559595/python-runtimewarning-overfow-encountered-in-long-scalars – Dimitri Merejkowsky Mar 04 '16 at 21:05
  • I recently updated my python to 3.5.1, I wonder if that has to do with the source of my RuntimeWarnings – Estefy Mar 04 '16 at 21:35

1 Answers1

0

There is a misunderstanding in the implementation of RK4. The state values for the intermediate points of the stages are not incremental to each other, but all based on the base state of the step. Thus it should be, for the ODE system th' = om, om'=f(t,th,om)

def f(t,th,om): return -q*om - sin(th) + b*cos(od*t)
for j in range(te):
     k_th1 = om[j]
     k_om1 = f(t[j],th[j],om[j])   

     #1 at mid-interval
     k_th2 = om[j] + (ts/2)*k_om1
     k_om2 = f(t[j]+(dt/2),th[j]+(dt/2)*k_th1,om[j]+(dt/2)*k_om1)

     #2 at mid-interval
     k_th3 = om[j] + (ts/2)*k_om2
     k_om3 = f(t[j]+(dt/2),th[j]+(dt/2)*k_th2,om[j]+(dt/2)*k_om2)

     #next time step
     k_th4 = om[j] + ts*k_om3
     k_om4 = f(t[j]+dt,th[j]+dt*k_th3,om[j]+dt*k_om3)

     dth=(k_om1 + 2*k_om2 + 2*k_om3 + k_om4)/6
     dom=(k_th1 + 2*k_th2 + 2*k_th3 + k_th4)/6
     th.append(th[j] + ts*dth) #integral of angle 
     om.append(om[j] + ts*dom) #integral of velocity
Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51