Integrating Ordinary Differential Equations (ODEs) with SciPy's solve_ivp

This guide provides a comprehensive overview of how to solve Ordinary Differential Equations (ODEs) in Python using the solve_ivp function from the scipy.integrate module. solve_ivp is a versatile and powerful tool for tackling Initial Value Problems (IVPs) of ODE systems.

What is an Ordinary Differential Equation (ODE)?

An Ordinary Differential Equation (ODE) describes the relationship between a dependent variable and its rate of change with respect to a single independent variable (commonly time).

The general form of an ODE is:

$$ \frac{dy}{dt} = f(t, y) $$

Where:

• $y(t)$: Represents the state variable, which is a function of the independent variable $t$.
• $t$: Is the independent variable, often representing time.
• $f(t, y)$: Is a function that defines the rate of change of $y$ with respect to $t$.

Key Concepts in ODEs

SciPy solve_ivp Function: Syntax and Parameters

The solve_ivp function offers a robust interface for solving ODEs.

Syntax

from scipy.integrate import solve_ivp

sol = solve_ivp(fun, t_span, y0, method='RK45', t_eval=None, dense_output=False, events=None, vectorized=False, args=None, **options)

Parameters

• fun:
  • Type: Callable.
  • Description: A function that computes the vector of the time derivative of the state vector $y$ from the current time $t$ and state vector $y$. It must have the signature fun(t, y, *args).
• t_span:
  • Type: Tuple of floats (t0, tf).
  • Description: The interval of integration. t0 is the initial time and tf is the final time.
• y0:
  • Type: ndarray.
  • Description: Initial state vector $y(t_0)$. The shape should be (n,) for a single ODE or (n, k) for $k$ systems of ODEs.
• method:
  • Type: String, optional.
  • Description: The integration method to use. Available methods include:
    • 'RK45' (default): Explicit Runge-Kutta method of order 5(4). Good general-purpose solver for non-stiff problems.
    • 'RK23': Explicit Runge-Kutta method of order 3(2). Also for non-stiff problems.
    • 'DOP853': Explicit Runge-Kutta method of order 8(5,3). High accuracy for non-stiff problems.
    • 'LSODA': Adams method for non-stiff and Gear method for stiff problems. Automatically switches.
    • 'BDF': Implicit backward differentiation formulas. Suitable for stiff problems.
    • 'Radau': Implicit Runge-Kutta method. Suitable for stiff problems with a high order of accuracy.
    • 'CDIPC': Implicit trapezoidal rule. Suitable for stiff problems.
• t_eval:
  • Type: ndarray, optional.
  • Description: Times at which to store the computed solution.
• dense_output:
  • Type: bool, optional.
  • Description: If True, a callable sol.sol is returned, which interpolates the solution at any time.
• events:
  • Type: Callable or list of callables, optional.
  • Description: Functions to evaluate at each time step. They are called with fun(t, y, *args) and must return a float. The solver stops when an event function crosses zero.
• vectorized:
  • Type: bool, optional.
  • Description: If True, fun is called with a 2D array of states y and the derivative is returned as a 2D array. Useful for systems where fun can be vectorized.
• args:
  • Type: Tuple, optional.
  • Description: Extra arguments to pass to the fun and events functions.
• **options:
  • Type: Dict, optional.
  • Description: Additional arguments passed to the chosen method, such as:
    • rtol: Relative tolerance.
    • atol: Absolute tolerance.
    • max_step: Maximum step size allowed for integration.
    • first_step: Initial step size.

Examples

Example 1: First-Order Exponential Decay

This example solves the ODE $\frac{dy}{dt} = -0.5y$ with initial condition $y(0) = 2$.

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt

def decay(t, y):
    """Defines the ODE for exponential decay."""
    return -0.5 * y

# Time span for integration
t_span = (0, 10)
# Initial condition y(0) = 2
y0 = [2]
# Time points where the solution is evaluated
t_eval = np.linspace(0, 10, 100)

# Solve the ODE
sol = solve_ivp(decay, t_span, y0, t_eval=t_eval)

# Plot the solution
plt.plot(sol.t, sol.y[0])
plt.title("Exponential Decay: $dy/dt = -0.5y$")
plt.xlabel("Time")
plt.ylabel("y(t)")
plt.grid(True)
plt.show()

Example 2: First-Order Linear ODE

This example solves the ODE $\frac{dy}{dt} = e^{-t} - 2y$ with initial condition $y(0) = 1$.

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt

def linear_ode(t, y):
    """Defines a first-order linear ODE."""
    return np.exp(-t) - 2 * y

# Time span and initial condition
t_span = (0, 5)
y0 = [1]
t_eval = np.linspace(0, 5, 100)

# Solve the ODE
sol = solve_ivp(linear_ode, t_span, y0, t_eval=t_eval)

# Plot the solution
plt.plot(sol.t, sol.y[0])
plt.title("First-Order Linear ODE: $dy/dt = e^{-t} - 2y$")
plt.xlabel("Time")
plt.ylabel("y(t)")
plt.grid(True)
plt.show()

Example 3: Second-Order ODE (Converted to First-Order System)

A second-order ODE can be converted into a system of two first-order ODEs. Consider $ \frac{d^2y}{dt^2} + 3\frac{dy}{dt} + 2y = 0 $.

Let $y_1 = y$ and $y_2 = \frac{dy}{dt}$. Then the system becomes:

1. $\frac{dy_1}{dt} = y_2$
2. $\frac{dy_2}{dt} = -3y_2 - 2y_1$

With initial conditions $y(0) = 1$ and $\frac{dy}{dt}(0) = 0$, which translates to $y_1(0) = 1$ and $y_2(0) = 0$.

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt

def system(t, Y):
    """Defines a system of first-order ODEs for a second-order ODE."""
    y, v = Y  # Y[0] is y, Y[1] is dy/dt
    dydt = v
    dvdt = -3 * v - 2 * y
    return [dydt, dvdt]

# Time span and initial conditions [y(0), dy/dt(0)]
t_span = (0, 10)
y0 = [1, 0]
t_eval = np.linspace(0, 10, 100)

# Solve the system of ODEs
sol = solve_ivp(system, t_span, y0, t_eval=t_eval)

# Plot the solutions
plt.plot(sol.t, sol.y[0], label="y(t)")
plt.plot(sol.t, sol.y[1], label="v(t) = dy/dt")
plt.title("Second-Order ODE: $d^2y/dt^2 + 3dy/dt + 2y = 0$")
plt.xlabel("Time")
plt.ylabel("State Variable")
plt.legend()
plt.grid(True)
plt.show()

Example 4: Predator-Prey Model (Lotka-Volterra Equations)

This example demonstrates solving a system of ODEs for a predator-prey model.

The Lotka-Volterra equations are:

• $\frac{dx}{dt} = \alpha x - \beta xy$ (Prey population change)
• $\frac{dy}{dt} = \delta xy - \gamma y$ (Predator population change)

Where:

• $x$: Prey population
• $y$: Predator population
• $\alpha$: Natural growth rate of prey
• $\beta$: Predation rate of prey
• $\delta$: Efficiency of converting prey into predator biomass
• $\gamma$: Natural death rate of predators

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt

def predator_prey(t, z, alpha, beta, delta, gamma):
    """Defines the Lotka-Volterra predator-prey equations."""
    x, y = z  # z[0] is prey, z[1] is predator
    dxdt = alpha * x - beta * x * y
    dydt = delta * x * y - gamma * y
    return [dxdt, dydt]

# Model parameters
alpha = 0.1
beta = 0.02
delta = 0.01
gamma = 0.1

# Time span and initial populations [prey(0), predator(0)]
t_span = (0, 200)
y0 = [40, 9]
t_eval = np.linspace(0, 200, 1000)

# Solve the system, passing parameters via 'args'
sol = solve_ivp(predator_prey, t_span, y0, args=(alpha, beta, delta, gamma), t_eval=t_eval)

# Plot the population dynamics
plt.plot(sol.t, sol.y[0], label="Prey")
plt.plot(sol.t, sol.y[1], label="Predator")
plt.title("Lotka-Volterra Predator-Prey Model")
plt.xlabel("Time")
plt.ylabel("Population")
plt.legend()
plt.grid(True)
plt.show()

Applications of ODE Integration

Numerical integration of ODEs is fundamental in many scientific and engineering disciplines:

• Physics: Simulating projectile motion, analyzing electrical circuits, modeling planetary orbits, quantum mechanics.
• Biology: Population dynamics (e.g., Lotka-Volterra), spread of diseases (epidemiological models), biochemical reaction kinetics.
• Finance: Option pricing, interest rate modeling, portfolio optimization.
• Engineering: Control systems design, signal processing, fluid dynamics, structural analysis, chemical reaction engineering.

Summary

SciPy's solve_ivp function provides a powerful and flexible framework for solving Ordinary Differential Equations in Python. It supports a variety of numerical integration methods, allowing users to select the most appropriate solver based on the characteristics of their ODEs (e.g., RK45 for non-stiff problems, BDF for stiff problems). By converting higher-order ODEs into systems of first-order equations, solve_ivp can handle a wide range of dynamic system models across various scientific and engineering domains.