Python numeric integrators can be used to solve equation system of the form

\[ \frac{d(y)}{dt} = f(t, y). \]

scipy.integrate.ode provides a generic interface to number of integrators, including

• vode;
• zvode;
• isoda;
• dopri5 (explicit runge-kutta of order 4); and
• dop853 (explicit runge-kutta of order 8).

For an up-to-date list check out scipy.integrate.ode documentation.

Example Usage

Consider a particle in free-fall. We can describe the motion of this particle using Newton’s equations of motion as follows:

\[ \frac{dx}{dt} = v \]

and

\[ \frac{dv}{dt} = \frac{f}{m}. \]

Here, \(x\) is the position (height) of the particle, \(v\) its velocity, \(f\) is the force acting on the particle, and \(m\) is its mass. Since the particle is in free fall, the force acting on it is \(m \times g\). Here, \(g\) is the acceleration due to gravity, and its value is \(= 9.8\ (meters)({seconds})^{-2}\). Say that the particle of mass \(100\ kilogram\) was released at a height of \(100\ meters\) at time \(0\ seconds\) (i.e., the initial velocity is \(0\ (meters)(second)^{-1}\)). We can set up the particle as follows.

import numpy as np

state = np.array([100.0, 0.0]) # Initial state
                               # state[0] = pos, state[1] = velocity
t = 0.0                        # time
mass = 100
g = 9.8

Next, we need to define the gradient function which captures the equations above. This function needs to return the derivatives of velocity and position w.r.t. time.

def f(t, state, arg1, arg2):
    # t = time
    # state[0] = position
    # state[1] = velocity
    # arg = extra arguments, where needed
    #    here, we are passing the value of g and mass m
    #    arg1 = mass
    #    arg2 = g
    
    force = - arg1 * arg2                        # -ve, since gravity points downwards

    dstate = np.array([ state[1], force/arg1)    # dstate[0] = velocity
                                                 # dstate[1] = force / mass
    
    return dstate

Now lets set up the ODE solver.

from scipy.integrate import ode

solver = ode(f)
solver.set_integrator('dop853')

We are almost there. Lets set the initial conditions. Don’t forget to also set up the parameters that needs to be passed to the grad function above.

solver.set_initial_value(state, t) # Recall that state contains
                                   # the initial state, t contains
                                   # the initial time

solver.set_f_params(m, g)          # These parameters will become available to
                                   # function f as arg1 and arg2

Now we integrate. Say we want to compute the location of the particle at time \(t + \delta t\). We do the following.

dt = 0.1   # seconds

if solver.successful():
    solver.integrate(solver.t + dt)
    
state = solver.y   # The new state of the particle
t = solver.t       # The new time (t + dt)

Done.

Code Listing

import numpy as np

state = np.array([100.0, 0.0]) # Initial state
                               # state[0] = pos, state[1] = velocity
t = 0.0                        # time
m = 100                        # mass
g = 9.8                        # acceleration due to gravity

def f(t, state, arg1, arg2):
    # t = time
    # state[0] = position
    # state[1] = velocity
    # arg = extra arguments, where needed
    #    here, we are passing the value of g and mass m
    #    arg1 = mass
    #    arg2 = g
    
    force = - arg1 * arg2                        # -ve, since gravity points downwards

    dstate = np.array([ state[1], force/arg1 ])  # dstate[0] = velocity
                                                 # dstate[1] = force / mass
    
    return dstate

from scipy.integrate import ode

solver = ode(f)
solver.set_integrator('dop853')

solver.set_initial_value(state, t) # Recall that state contains
                                   # the initial state, t contains
                                   # the initial time

solver.set_f_params(m, g)          # These parameters will become available to
                                   # function f as arg1 and arg2

dt = 0.1   # seconds

for i in range(10):
    if solver.successful():
        solver.integrate(solver.t + dt)
        
    state = solver.y   # The new state of the particle
    t = solver.t       # The new time (t + dt)

    print('Time', t)
    print('State', state)