scipy.integrate — odeint, solve_ivp (RK45, DOP853), quad

1. quad — numerical definite integral

Deriving the simplest rule from scratch (trapezoid → Simpson)

We want $\int_a^b f\,dx$ using only $f(a), f(m), f(b)$ with $m=\frac{a+b}{2}$. Fit a parabola $p(x)=\alpha + \beta x + \gamma x^2$ through these 3 points and integrate it exactly. Shift so $a=-h,\ m=0,\ b=h$. Then:

$\int_{-h}^{h} (\alpha+\beta x+\gamma x^2)\,dx = 2\alpha h + \tfrac{2}{3}\gamma h^3$

The odd term vanishes (Why? $\int_{-h}^h x\,dx=0$ by symmetry). Now express $\alpha,\gamma$ via the samples: $f(0)=\alpha$, and $f(h)+f(-h)=2\alpha+2\gamma h^2 \Rightarrow \gamma h^2 = \tfrac{f(h)+f(-h)}{2}-f(0)$. Substitute:

$\boxed{\int_{-h}^{h} f \approx \frac{h}{3}\big[f(-h)+4f(0)+f(h)\big]}$

That is Simpson's rule. quad is this idea, but adaptive and higher order.

from scipy.integrate import quad
import numpy as np
 
val, err = quad(lambda x: np.exp(-x**2), 0, np.inf)
# val ≈ 0.8862269  (= sqrt(pi)/2), err ≈ 7.1e-09

• val = estimate, err = estimated absolute error.
• Pass extra params with args=(...); use points=[...] to flag kinks.

2. solve_ivp — the modern ODE solver

Deriving a single step (forward Euler, then RK)

Start from the exact statement (fundamental theorem of calculus): $y(t_{n}+h) = y(t_n) + \int_{t_n}^{t_n+h} f(t,y)\,dt$

WHY this is the seed: if we knew that integral exactly we'd be done; every method is just a quadrature rule for that little integral.

• Euler (left-rectangle): approximate the integral by $h\,f(t_n,y_n)$: $y_{n+1}=y_n + h f(t_n,y_n)\quad\text{(error }O(h^2)\text{ per step)}$

• RK4 / RK45 (Runge–Kutta): sample the slope at several trial points inside the step, like Simpson did inside an interval, and blend: $y_{n+1}=y_n+\tfrac{h}{6}(k_1+2k_2+2k_3+k_4)$ with $k_1=f(t_n,y_n)$, $k_2=f(t_n+\tfrac h2, y_n+\tfrac h2 k_1)$, etc.

Worked example — exponential decay $\dot y=-2y,\ y(0)=3$

Exact solution: $y(t)=3e^{-2t}$.

from scipy.integrate import solve_ivp
import numpy as np
 
sol = solve_ivp(lambda t, y: -2*y, [0, 5], [3.0],
                method='RK45', t_eval=np.linspace(0, 5, 100),
                rtol=1e-8, atol=1e-10)
sol.t      # times
sol.y[0]   # y values, shape (n_states, n_times)

• Why y0=[3.0] (a list)? solve_ivp always treats the state as a vector; even a scalar must be 1-element.
• Why sol.y[0]? rows are state components, columns are time points.

Worked example — a system (SHM / 2nd-order → 1st-order)

$\ddot x = -\omega^2 x$. Why convert? Solvers only take first-order systems. Let $y_0=x,\ y_1=\dot x$: $\dot y_0 = y_1,\qquad \dot y_1=-\omega^2 y_0$

w = 2.0
def rhs(t, y):
    return [y[1], -w**2 * y[0]]
sol = solve_ivp(rhs, [0, 10], [1, 0], method='DOP853', rtol=1e-10)
# x(t) = sol.y[0] should match cos(2t); energy conserved

3. odeint — the legacy classic

from scipy.integrate import odeint
y = odeint(lambda y, t: -2*y, 3.0, np.linspace(0,5,100))  # note (y, t)!

• Returns an array of shape (len(t), n_states) — also transposed vs solve_ivp.
• Modern advice: prefer solve_ivp (events, dense output, method choice). Keep odeint knowledge for reading old code.

Flashcards

Connections

• Riemann sums and Fundamental Theorem of Calculus — the seed of all integration here.
• Runge-Kutta methods · Stiff differential equations · Numerical stability
• Gaussian quadrature — the math behind quad.
• numpy.trapz and simpson — fixed-grid integration when you only have samples.
• Harmonic oscillator — canonical solve_ivp test system.

Concept Map

Hinglish (regional understanding)