How to map a roller coaster differential equation to scipy?

Question

I want to replicate this fantastic blog post in Scipy, which shows how to simulate roller coaster physics using this ODE system

with

p = position on the track (measured by path length along the track)
v = velocity
k(p) = slope at point p (in code denoted as slope(p))
g = gravity
(-b/m) v = dampening term (not in code)

I don't care about the dampening. Just the general concept.

I tried to replicate it in scipy with a sample track that just goes at a constant slope down (slope = -1)

import numpy as np
from scipy import integrate

def slope(p):
    return -1

def f(p, U):
    return [
            U[1],
            (-9.81 * slope(p))/ np.sqrt(1 + slope(p)**2)
            ]

sol = integrate.solve_ivp(f, (0, 5), [0, 0])
print(sol)

As I understand it, the returned function f maps time to postion p and velocity: f(time) -> [position, velocity].

Result:

t: [0.0000e+00, 1.0000e-04, 1.1000e-03, 1.1100e-02, 1.1110e-01, 1.1111e+00, 5.0000e+00]
y[0]: [0.00000000e+00, 3.46835876e-08, 4.19671410e-06, 4.27336483e-04, 4.28106806e-02, 4.28183876e+00, 8.67089690e+01]
y[1]: [0.00000000e+00, 6.93671752e-04, 7.63038928e-03, 7.69975645e-02, 7.70669317e-01, 7.70738684e+00, 3.46835876e+01]]

However, this is not what I expected.

Using the potential/kinectic energy equations and rearranging to v: v = sqrt(2*g*h)

I would expect that after 5m at slope of -1 (meaning going down by 5 meters), the velocity to be:

v = sqrt(2 * 9.81 * 5) = 9.904 (m/s)

but the ode says it's 34 (m/s)

What am I doing wrong?

P.S. I am aware that g might need to changed to -g but the numbers don't add up either way.

Answer 1

Since k(p) = -1, your acceleration is g/sqrt(2) not g

Also, you are calling a function passing a time interval of 5 seconds. The 34 m/s are attained after 5 seconds not 5 meters = 9.81 / sqrt(2) x 5, which is correct.