Search code examples
pythonscipyodedifferential-equations

How to solve a differential equation with a non-constant coefficient using Python ？

Suppose there is an equation about the length of the bungee cord (denoted by x), which is dependent on the mass of the object e.g. a player (denoted by m).

Assume that the natural length of the bungee cord is 30 meters, in other words, the starting position is x(0)=-30.

The equation of the length of the bungee cord is given by:

    x''(m) = g + b/m*x(m) -a1/m*x'(m) - a2*|x'(m)|*x'(m)

where g, a1, a2 are constants; b is a step function: b = -k (another constant) when x<0 and b = 0 when x>=0.

import numpy as np
from scipy.integrate import odeint

g = 9.8
a1, a2 = 0.6, 0.8
k = 20
b = [-k, 0]

def dx_dt(x, m):
    if x[0]>=0:
        return [x[1], g+b[0]/m*x[0]-a1/m*x[1]-a2/m*np.abs(x[1])*x[1]]
    else:
        return [x[1], g+b[1]/m*x[0]-a1/m*x[1]-a2/m*np.abs(x[1])*x[1]]

init = [[-30, 0], [-40, 0.0001]]

m = np.linspace(1, 100, 10000)

fig, ax = plt.subplots(1, 2, sharey=True, figsize=(6, 4))

for i in range(len(init)):
    xs = odeint(dx_dt, init[i], m)
    ax[i].plot(m, xs[:, 0], 'r-')
    ax[i].set_xlabel('mass (m)')
    ax[i].set_ylabel('length (x)')
    ax[i].set_xlim(xmin, xmax)
    ax[i].set_ylim(ymin, ymax)
    ax[i].set_title('init={}'.format(init[i]))

the right answer should be a sine-like curve

but the result from above codes turns out to be

enter image description here

Is there something wrong with the codes?

Solution

  • Change the length coordinate x to point upwards, the cord without jumper being at rest at position 0 so that for x<0 the cord behaves like a spring. Then also gravity points downwards. The modified ODE function for this is 

    def dx_dt(x, t, m):
    acc = -g-a1/m*x[1]-a2/m*np.abs(x[1])*x[1]
    if x[0] < 0: acc -= k/m*x[0]
    return [x[1], acc]

    Plotting for 3 different masses 

    for i,ini in enumerate(init):
    for m in masses: 
        xs = odeint(dx_dt, ini, t, args=(m,))
        ax[i].plot(t, xs[:, 0], label="m=%.2f"%m)

    gives then the picture

    enter image description here

    If, as example, the ground level in the first situation is at -80m, giving a total height of 110m then the mass of the jumper has to be less than 90kg. For more precise statements use inter- and extrapolation or a numerical solver to find the first time where x'(t)=0 and the critical mass where x(t)=ground at that time.

    It appears clear that in no case the jumper does re-enter the "free-fall" phase x>0.