Reputation: 11
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
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