SciPy - 3 (solving ordinary differential equations)

SciPy provides two different ways to solve ODEs: An API based on the function odeint, and object-oriented API based on the class ode. Usually odeint is easier to get started with, but the ode class offers some fi ner level of control. To solve an ordinary differential equation of the type-



with a given y(t0) = y0 we can use scipy's odeint function. The following program (useodeint.py) to find y(t) for t 2 [0; 2] given this differential equation:



See the code below:

from scipy . integrate import odeint
import numpy as N
def f(y, t):
    """ this is the rhs of the ODE to integrate , i.e. dy/dt=f(y,t)"""
    return -2 * y * t
y0 = 1 # initial value
a = 0 # integration limits for t
b = 2
t = N. arange (a, b, 0.01) # values of t for
# which we require
# the solution y(t)
y = odeint (f, y0 , t) # actual computation of y(t)
import pylab # plotting of results
pylab . plot (t, y)
pylab . xlabel ('t'); pylab . ylabel ('y(t)')
pylab . show ()

The output of the program is as shown below:



The odeint command takes a number of optional parameters to change the default error tolerance
of the integration (and to trigger the production of extra debugging output).  Here we have used the odeint functions. For more information about the class ode, try help(ode). It does pretty much the same thing as odeint, but in an object-oriented fashion.

ODE problems are important in computational physics, so we will look at one more example: the damped harmonic oscillation.

Damping is an influence within or upon an oscillatory system that has the effect of reducing, restricting or preventing its oscillations. In physical systems, damping is produced by processes that dissipate the energy stored in the oscillation. Examples include viscous drag in mechanical systems, resistance in electronic oscillators, and absorption and scattering of light in optical oscillators. Damping not based on energy loss can be important in other oscillating systems such as those that occur in biological systems and bikes.


The damping ratio is a dimensionless measure describing how oscillations in a system decay after a disturbance. Many systems exhibit oscillatory behavior when they are disturbed from their position of static equilibrium. A mass suspended from a spring, for example, might, if pulled and released, bounce up and down. On each bounce, the system tends to return to its equilibrium position, but overshoots it. Sometimes losses (e.g. frictional) damp the system and can cause the oscillations to gradually decay in amplitude towards zero or attenuate. The damping ratio is a measure describing how rapidly the oscillations decay from one bounce to the next.

The damping ratio is a system parameter, denoted by ζ (zeta), that can vary from undamped (ζ = 0), underdamped (ζ < 1) through critically damped (ζ = 1) to overdamped (ζ > 1).

The behaviour of oscillating systems is often of interest in a diverse range of disciplines that include control engineering, chemical engineering, mechanical engineering, structural engineering, and electrical engineering. The physical quantity that is oscillating varies greatly, and could be the swaying of a tall building in the wind, or the speed of an electric motor, but a normalised, or non-dimensionalised approach can be convenient in describing common aspects of behavior.

The equation of motion for the damped oscillator is:



where x is the position of the oscillator, ω0 is the frequency, and ζ is the damping ratio. To write this
second-order ODE on standard form we introduce p = dx/dt :



In the implementation of this example we will add extra arguments to the RHS function for the ODE,
rather than using global variables as we did in the previous example. As a consequence of the extra arguments to the RHS, we need to pass an keyword argument args to the odeint function. See the following code:

from scipy . integrate import odeint
import numpy as N
from matplotlib import pyplot as plt

def dy(y, t, zeta, w0):
    """
    The right-hand side of the damped oscillator ODE
    """
    x, p = y[0], y[1]
    dx = p
    dp = -2 * zeta * w0 * p - w0**2 * x
    return [dx, dp]

# initial state:
y0 = [1.0, 0.0]

# time coodinate to solve the ODE for
t = N.linspace(0, 10, 1000)
w0 = 2*N.pi*1.0

# solve the ODE problem for three different values of the damping ratio
y1 = odeint(dy, y0, t, args=(0.0, w0)) # undamped
y2 = odeint(dy, y0, t, args=(0.2, w0)) # under damped
y3 = odeint(dy, y0, t, args=(1.0, w0)) # critial damping
y4 = odeint(dy, y0, t, args=(5.0, w0)) # over damped

fig, ax = plt.subplots()
ax.plot(t, y1[:,0], 'k', label="undamped", linewidth=0.25)
ax.plot(t, y2[:,0], 'r', label="under damped")
ax.plot(t, y3[:,0], 'b', label=r"critical damping")
ax.plot(t, y4[:,0], 'g', label="over damped")
plt.show()
ax.legend()

In this program we have use matplotlib a mathematical plotting library which is used to make
simple plots, such as line graphs and scatter plots., which we shall discuss separately. The output of the above program is as shown:



Here I am ending this post regarding SciPy solving ordinary differential equations. Next post will discuss about Fourier transforms. Till we meet next keep learning and practicing Python as Python is easy to learn!