lunazagor
lunazagor

Reputation: 11

Pendulum-like odeint integration

I'm trying to solve a pendulum-like differential equation that goes like d2(phi)/dt = -(g/R) *sin(phi) (it's the skateboard problem in Taylor's Classical Mechanics). I'm new to scipy and odeint and the like, so I'm doing this to prepare for more sophisticated numerical solutions in the future.

I've used code from here to try and navigate the coding, but all I'm coming up with is a flat line for phi(t). I think it stems from the fact that I'm trying to split a second order differential equation into two first orders, where one doesn't depend on the other (because d(phi)/dt to make an appearance); but I'm not sure how to go about fixing it.

Anyone know what's wrong with it?

    # integrate skateboard problem, plot result
    from numpy import *
    from scipy.integrate import odeint 
    import matplotlib.pyplot as plt

    def skate(y, t, params):
        phi, omega = y
        g, R = params
        derivs = [omega, -(g/R)*np.sin(phi)]
        return derivs

   # Parameters
   g = 9.81
   R = 5
   params = [g, R]

   #Initial values
   phi0 = 20
   omega0 = 0
   y0 = [phi0, omega0]

   t = linspace(0, 20, 5000)

   solution = odeint(skate, y0, t, args=(params,))

   plt.plot(t, solution[:,0])
   plt.xlabel('time [s]')
   plt.ylabel('angle [rad]')
   plt.show()

Upvotes: 1

Views: 1043

Answers (1)

pptaszni
pptaszni

Reputation: 8345

I suspect a bug here: -(g/R)*np.sin(phi). Perhaps you forget to define the alias for numpy lib import (for example: import numpy as np). Instead you just did (from numpy import *). Try this:

def skate(y, t, params):
    phi, omega = y
    g, R = params
    derivs = [omega, -(g/R)*sin(phi)]
    return derivs

Upvotes: 1

Related Questions