Ordinary Differential Equations in 1D

Numerical Integration with SciPy

Numerical integration is fundamental for solving ordinary differential equations (ODEs) that don’t have analytical solutions. In this first session, you will learn how to:

Contents

The Initial Value Problem

An initial value problem (IVP) consists of:

\[\frac{dy}{dt} = f(t, y), \quad y(t_0) = y_0\]

Where: - \(f(t, y)\) is the rate of change function - \(y_0\) is the initial condition at time \(t_0\)

SciPy’s solve_ivp

The scipy.integrate.solve_ivp function is the standard tool for solving ODEs in Python. Try the following code:

import numpy as np
from scipy.integrate import solve_ivp

def exponential_decay(t, y):
     return -0.5 * y

t_span = [0, 10]
y0 = [2, 4, 8]
sol = solve_ivp(
    fun=exponential_decay,
    t_span=t_span, 
    y0=y0)

print(sol.t)

print(sol.y)

How does solve_ivp work? Let’s understand its parameters (copied from the documentation):

  • fun: Right-hand side of the system: the time derivative of the state \(y\) at time \(t\). The calling signature is fun(t, y), where t is a scalar and y is an ndarray with len(y) = len(y0). Additional arguments need to be passed if args is used (see documentation of args argument). fun must return an array of the same shape as y. In our example, exponential_decay defines the ODE \(\frac{dy}{dt} = -0.5y\), which models exponential decay \(y(t) = y_0 e^{-0.5t}\).
  • t_span: Interval of integration \((t_0, t_f)\). The solver starts with \(t=t_0\) and integrates until it reaches \(t=t_f\). Both \(t_0\) and \(t_f\) must be floats or values interpretable by the float conversion function.
  • y0: Initial state. For problems in the complex domain, pass y0 with a complex data type (even if the initial value is purely real).

You can also specify additional parameters: - method: Integration method to use. Common choices include 'RK45' (default), 'RK23', 'DOP853', 'Radau', 'BDF', and 'LSODA'. - t_eval: Times at which to store the computed solution, must be sorted and lie within t_span. If None (default), use points selected by the solver. - args: Additional arguments to pass to the user-defined functions. If, for example, fun has the signature fun(t, y, a, b, c), then args=(a, b, c).

Systems of ODEs

For multiple coupled equations, return a list or array of derivatives:

def sir_model(t, y, beta, gamma):
    S, I, R = y
    N = S + I + R
    
    dSdt = -beta * S * I / N
    dIdt = beta * S * I / N - gamma * I
    dRdt = gamma * I
    
    return [dSdt, dIdt, dRdt]

Best Practices

  1. Always check convergence: Plot solutions at different tolerances
  2. Use appropriate methods: 'RK45' (default) works well for most problems
  3. Vectorize when possible: Makes code faster and cleaner
  4. Document parameters: Keep track of units and meanings
  5. Validate against known solutions: Test your implementation

Next Steps

Apply these techniques to the classical models in Session 1: - SIR epidemiological model - Spruce budworm population dynamics - Michaelis–Menten enzyme kinetics

Resources