Differential Equations: Initial & Boundary Value Problems

Dynamic systems in physics and engineering—heat flow, chemical kinetics, vibrations, and orbital mechanics—are modeled as ordinary differential equations (ODEs).

1. Initial Value Problems: solve_ivp

scipy.integrate.solve_ivp is the standard solver for systems of the form: dy/dt = f(t, y), y(t0) = y0

Selecting the Right Integration Method

Solving a Damped Harmonic Oscillator

Convert the second-order equation $m x'' + c x' + k x = 0$ into a first-order system with state $\mathbf{z} = [x, v]$: $z_0' = z_1$ $z_1' = -(k/m) z_0 - (c/m) z_1$

2. Event Detection: Stopping at Exact Crossings

Events allow you to locate when a condition $g(t, y) = 0$ occurs—such as a bouncing ball hitting the ground or a cooling tank reaching target temperature:

3. Stiff Differential Equations (Radau, BDF)

A system is stiff when fast and slow dynamics occur simultaneously. An explicit solver like 'RK45' is forced to take tiny steps to stay numerically stable, requiring hundreds of thousands of steps. Implicit solvers ('Radau', 'BDF') step over fast decaying transients safely:

4. Boundary Value Problems: solve_bvp

In BVPs, boundary conditions are specified at multiple spatial locations. scipy.integrate.solve_bvp uses fourth-order collocation on an adaptive mesh:

Solve $y'' + y = 0$ with $y(0) = 0$ and $y(\pi / 2) = 1$ (analytical solution $y(x) = \sin(x)$):

5. Try It Yourself

Problem: Radioactive Decay with solve_ivp

A radioisotope decays according to $d N / dt = -\lambda N$ with decay constant $\lambda = 0.2\text{ s}^{-1}$. Starting from $N_0 = 500$ atoms at $t = 0$, compute the population at $t = 5\text{ s}$.

import numpy as np
from scipy.integrate import solve_ivp
 
lam = 0.2
sol = solve_ivp(lambda t, N: [-lam * N[0]], (0.0, 5.0), [500.0], method='RK45')
print(f"Population at t=5s: {sol.y[0, -1]:.2f}")