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

This page is the "walk through every door" companion to the parent note. There we built the tools. Here we stress-test them: every sign, every degenerate input, every limiting case, one word problem, one exam trap.

Intuition Assumed imports for every code cell

To keep each snippet short, assume these two lines already ran at the top of the notebook:

import numpy as np
from scipy.integrate import quad, solve_ivp, odeint

So every np.sin, np.pi, np.inf, quad, solve_ivp, odeint below refers to these.

Definition The four classic pitfalls (self-contained, referenced below)

Because this page must stand alone, here are the four solve_ivp/odeint traps we will point back to:

P1 — signature swap. solve_ivp wants f(t, y) (time first); odeint wants f(y, t) (state first).
P2 — trusting the default tolerance. The default rtol=1e-3 can be visibly wrong for oscillators or long integrations; set rtol/atol yourself when accuracy matters.
P3 — confusing t_eval with step size. t_eval only chooses where the answer is reported; the solver picks its own internal adaptive steps.
P4 — explicit method on a stiff problem. RK45/DOP853 are forced into microscopic steps for stability on stiff systems; switch to implicit BDF/Radau/LSODA.

One more fact we will lean on, restated so this page stands alone: solve_ivp, odeint and every Runge–Kutta method only accept first-order systems — a 2nd-order equation like $\overset{x}{¨} = \dots$ must first be rewritten as a pair of 1st-order equations (see Runge-Kutta methods).

The scenario matrix

Two engines live in this module, so we enumerate cases for each. Examples are numbered flat 1–12 so nothing repeats.

Cell	What makes it special	Example
Q-finite	`quad` on a plain finite interval $[a, b]$	Ex 1
Q-infinite	`quad` with one $\infty$ limit	Ex 2
Q-endpoint-singular	integrand blows up at an endpoint, area finite	Ex 3
Q-signed	integrand goes negative → signed area, cancellation	Ex 4
Q-degenerate	$a = b$ (zero width) and $a > b$ (reversed)	Ex 5
Q-interior-singular	blow-up at an interior point $c$ → split / `points=`	Ex 6
Q-doubly-infinite	interval $(- \infty, \infty)$ , symmetric	Ex 7
ODE-scalar-decay	one equation, solution shrinks	Ex 8
ODE-system-oscillate	2nd order → system, energy conserved	Ex 9
ODE-stiff	two wildly different time-scales	Ex 10
ODE-word-problem	real-world model (cooling)	Ex 11
ODE-exam-trap	`odeint` vs `solve_ivp` signature swap	Ex 12

Every cell below is worked to the end and every number is machine-checked in the verify block.

Ex 1 — Cell Q-finite: a polynomial, no surprises

Forecast: guess the number before reading. (Hint: the antiderivative is $x^{3} + x$ .)

Recognise the exact answer first. $\int_{0}^{2} (3 x^{2} + 1) d x = [x^{3} + x]_{0}^{2} = 8 + 2 = 10$ . Why this step? Whenever an exact answer exists, compute it — it becomes the yardstick for the numeric one.
Call quad.
```
val, err = quad(lambda x: 3*x**2 + 1, 0, 2)
# val = 10.0, err ≈ 1.1e-13
```
Why this step? quad returns a pair: the estimate and its own guess of the error. Never throw away err.

Ex 2 — Cell Q-infinite: one limit goes to infinity

Forecast: an infinite interval — will quad even survive? What's the exact area?

Exact value. $\int_{0}^{\infty} e^{- x} d x = [- e^{- x}]_{0}^{\infty} = 0 - (- 1) = 1$ . Why this step? The tail $e^{- x}$ decays to zero fast, so the area is finite even though the interval is not. That's the condition for quad to work with np.inf.
Call it with np.inf.
```
val, err = quad(lambda x: np.exp(-x), 0, np.inf)
# val ≈ 1.0, err ≈ 5e-11
```
Why this step? quad internally maps the infinite ray onto a finite interval (a change of variable) before applying Gauss–Kronrod. You never see it; you just pass np.inf.

Ex 3 — Cell Q-endpoint-singular: integrand blows up at an endpoint

Forecast: at $x = 0$ the integrand is $+ \infty$ . Does the area explode too, or converge?

Test convergence with the exact integral. $\int_{0}^{1} x^{- 1/2} d x = [2 x]_{0}^{1} = 2 - 0 = 2$ . Why this step? An integrable singularity: the spike is infinitely tall but infinitely thin, and the areas add to a finite $2$ . Contrast $1/ x$ , whose integral would diverge.
Call quad — it handles the endpoint singularity.
```
val, err = quad(lambda x: x**-0.5, 0, 1)
# val ≈ 2.0, err ≈ 3e-11
```
Why this step? The adaptive part of "adaptive Gauss–Kronrod" clusters sample points near $x = 0$ where the function wiggles violently, so it still converges.

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

Figure (s01): the red curve $1/ x$ rockets to infinity at the left wall, yet the blue shaded area is bounded — it fills up to the yellow dashed line at height $2$ , the exact value of the integral.

Look at the figure: the red curve rockets up near the left wall, yet the shaded blue area is bounded by the dashed line at height $2$ .

Ex 4 — Cell Q-signed: the integrand dips negative

Forecast: $sin$ is positive on $[0, π]$ and negative on $[π, 2 π]$ . Guess both answers.

Full period first. $\int_{0}^{2 π} sin x d x = [- cos x]_{0}^{2 π} = - 1 - (- 1) = 0$ . Why this step? quad computes signed area. The hill on $[0, π]$ and the equal valley on $[π, 2 π]$ cancel exactly.
Half period. $\int_{0}^{π} sin x d x = [- cos x]_{0}^{π} = 1 - (- 1) = 2$ . Why this step? This isolates the positive hump so nothing cancels — shows the difference between signed area and "how much curve there is."
```
v_full, _ = quad(np.sin, 0, 2*np.pi)   # ≈ 0 (like 1e-16)
v_half, _ = quad(np.sin, 0, np.pi)      # ≈ 2.0
```

Figure (s02): sin(x) over one full period. The green lobe above the axis contributes +2, the red lobe below contributes −2; the signed total is 0.

Green area (above axis) $= + 2$ , red area (below axis) $= - 2$ ; together $0$ .

Ex 5 — Cell Q-degenerate: zero width and reversed limits

Forecast: one has zero width, one has the limits backwards. Predict both signs.

Zero-width interval. $\int_{3}^{3} x^{2} d x = 0$ . Why this step? Width $b - a = 0$ ⇒ no area regardless of the function. quad(lambda x: x**2, 3, 3) returns $0.0$ . This is the degenerate edge case that catches people who feed in equal limits by accident.
Reversed limits. By the definition $\int_{a}^{b} = - \int_{b}^{a}$ : $\int_{5}^{1} x^{2} d x = - \int_{1}^{5} x^{2} d x = - [\frac{x ^{3}}{3}]_{1}^{5} = - (\frac{125}{3} - \frac{1}{3}) = - \frac{124}{3}$ . Why this step? quad honours the orientation of the interval, so a backwards interval flips the sign. $- 124/3 \approx - 41.333$ .
```
quad(lambda x: x**2, 3, 3)   # (0.0, 0.0)
quad(lambda x: x**2, 5, 1)   # (-41.333..., err)
```

Ex 6 — Cell Q-interior-singular: blow-up inside the interval

Forecast: the spike is now in the middle, not at an endpoint. Will a plain quad(f, -1, 1) notice it, and what is the exact area?

Exact value by symmetry. The integrand is even, so $\int_{- 1}^{1} ∣ x ∣^{- 1/2} d x = 2 \int_{0}^{1} x^{- 1/2} d x = 2 \cdot 2 = 4$ . Why this step? We reuse Ex 3's endpoint result on each half, doubling by symmetry. The area is finite: $4$ .
Warn quad where the trouble is with points=.
```
val, err = quad(lambda x: abs(x)**-0.5, -1, 1, points=[0])
# val ≈ 4.0, err ≈ 1e-9
```
Why this step? quad places its samples adaptively, but it can miss a razor-thin interior spike because no sample landed near $x = 0$ . Passing points=[0] forces it to split the interval at the kink, treating $0$ as an endpoint of each half where singularities are handled cleanly. Equivalently you could do quad(f,-1,0)[0] + quad(f,0,1)[0] by hand. Why not just quad(f,-1,1)? Without points, quad may return a value with a large err or a warning, because it never adapted near the interior blow-up.

Ex 7 — Cell Q-doubly-infinite: the whole real line

Forecast: both limits are infinite. Does quad accept -np.inf and np.inf together, and what famous number comes out?

Exact value. This is the classic Gaussian integral: $\int_{- \infty}^{\infty} e^{- x^{2}} d x = π \approx 1.7725$ . Why this step? The tails decay like $e^{- x^{2}}$ (faster than any exponential), so the area over the whole line is finite — the precondition for a doubly-infinite quad.
Pass both infinities.
```
val, err = quad(lambda x: np.exp(-x**2), -np.inf, np.inf)
# val ≈ 1.7724539, err ≈ 1e-08
```
Why this step? quad maps both infinite tails onto finite intervals internally. You literally write -np.inf, np.inf; no manual splitting needed for this smooth, symmetric integrand.

Ex 8 — Cell ODE-scalar-decay: solution shrinks toward zero

Forecast: decay with rate 2. Roughly how big is $y (1)$ ?

Exact solution. $\overset{y}{˙} = - 2 y \Rightarrow y (t) = 3 e^{- 2 t}$ , so $y (1) = 3 e^{- 2} \approx 0.406$ . Why this step? Linear decay has a closed form; use it as the ground truth for the numeric run.

Solve numerically with tight tolerances.

sol = solve_ivp(lambda t, y: -2*y, [0,5], [3.0],
                method='RK45', t_eval=[0,1,2,3,4,5],
                rtol=1e-10, atol=1e-12)
sol.y[0][1]   # y at t=1 ≈ 0.4060

Why this step? We set rtol/atol tightly (pitfall P2 above) so the numeric

y (1)

matches the analytic

0.406

closely. t_eval only picks where we report, not the internal step size (pitfall P3).

Ex 9 — Cell ODE-system-oscillate: 2nd order → system, energy conserved

Forecast: with those initial conditions the motion is a pure cosine. What is $x (π)$ ?

Reduce to first order. Let $y_{0} = x$ , $y_{1} = \overset{x}{˙}$ . Then $\overset{y}{˙}_{0} = y_{1}, \overset{y}{˙}_{1} = - ω^{2} y_{0} .$ Why this step? As restated in the pitfalls box, solvers accept only first-order systems. We trade one 2nd-order equation for two 1st-order ones stacked in a vector.
Analytic solution. $x (t) = cos (2 t)$ , so $x (π) = cos (2 π) = 1$ . Why this step? $ω = 2$ , zero initial velocity ⇒ pure cosine amplitude 1. Gives us the checkpoint.

Energy invariant.

E = \frac{1}{2} \overset{x}{˙}^{2} + \frac{1}{2} ω^{2} x^{2}

. At

t = 0

E = \frac{1}{2} (0) + \frac{1}{2} (4) (1) = 2

. It must stay

2

forever. Why this step? Conservation is a physics-level sanity check independent of the solver's internals — if

E

drifts, the integration is inaccurate.

w = 2.0
rhs = lambda t, y: [y[1], -w**2*y[0]]
sol = solve_ivp(rhs, [0,10], [1,0], method='DOP853',
                t_eval=[np.pi], rtol=1e-11, atol=1e-13)
x = sol.y[0][0]                       # ≈ 1.0
E = 0.5*sol.y[1][0]**2 + 0.5*w**2*x**2 # ≈ 2.0

Figure (s03): blue is position $x (t) = cos 2 t$ , yellow is velocity $\overset{x}{˙} (t) = - 2 sin 2 t$ , green is total energy — a flat line at 2, showing the high-order solver conserves it.

Blue is $x (t) = cos 2 t$ , yellow is $\overset{x}{˙} (t) = - 2 sin 2 t$ ; the green energy line stays flat.

Ex 10 — Cell ODE-stiff: two clashing time-scales

Forecast: the $- 1000$ factor pulls $y$ toward $cos t$ almost instantly. Which solver copes, and which one takes thousands of steps?

Spot the exact particular solution. Try $y = cos t$ : then $\overset{y}{˙} = - sin t$ and the RHS becomes $- 1000 (cos t - cos t) - sin t = - sin t$ . It matches, and $y (0) = cos 0 = 1$ … but our start is $y (0) = 0$ , so there's a fast transient that snaps up to $cos t$ within $\sim 1/1000$ of a second, after which $y (t) \approx cos t$ . Why this step? This is the definition of stiff (see Stiff differential equations): a superfast decay ( $1/1000$ s) riding on a slow oscillation ( $2 π$ s).
Run RK45 and inspect its work. An explicit method needs $h ≲ 2/1000$ for stability, so it takes thousands of tiny steps.
```
f = lambda t, y: -1000*(y - np.cos(t)) - np.sin(t)
s1 = solve_ivp(f, [0,2], [0.0], method='RK45', rtol=1e-8, atol=1e-10)
print(s1.nfev, s1.t.size)   # e.g. ~7000+ RHS evals, ~1000+ accepted steps
print(s1.status, s1.message) # 0, 'The solver successfully reached ...'
```
Why this step? sol.nfev (number of RHS evaluations) and sol.t.size (accepted steps) are the diagnostics that substantiate the "grinds to tiny steps" claim — this is pitfall P4 in numbers, not just words. It succeeds, but at huge cost.

BDF (implicit) sails through in a handful of steps.

s2 = solve_ivp(f, [0,2], [0.0], method='BDF', rtol=1e-8, atol=1e-10)
print(s2.nfev, s2.t.size)   # dramatically fewer than RK45
s2b = solve_ivp(f, [0,2], [0.0], method='BDF',
                t_eval=[1.0], rtol=1e-8, atol=1e-10)
s2b.y[0][0]                 # ≈ cos(1) = 0.5403

Why this step? After the transient,

y \approx cos t

; at

t = 1

cos 1 = 0.5403

. Implicit BDF is stable at large

h

regardless of the

1000

, so its nfev/step count is a small fraction of RK45's.

Ex 11 — Cell ODE-word-problem: Newton's law of cooling

Forecast: it heads toward room temperature $20$ . After 10 min will it be near 20, or still hot?

Model in solver form. State is scalar $T$ ; RHS is $f (t, T) = - 0.1 (T - 20)$ . Why this step? Word problems must be translated into $\dot{T} = f (t, T)$ with a numeric initial condition $T (0) = 90$ .
Exact solution for the checkpoint. $T (t) = 20 + (90 - 20) e^{- 0.1 t} = 20 + 70 e^{- 0.1 t}$ . At $t = 10$ : $T = 20 + 70 e^{- 1} = 20 + 70 (0.36788) = 45.750 4^{\circ}$ C. Why this step? Newton cooling is linear ⇒ closed form; use it to validate.
```
f = lambda t, T: -0.1*(T - 20)
sol = solve_ivp(f, [0,10], [90.0], t_eval=[10],
                rtol=1e-9, atol=1e-11)
sol.y[0][0]   # ≈ 45.7504
```

Ex 12 — Cell ODE-exam-trap: the `odeint` vs `solve_ivp` signature swap

Forecast: which library wants f(t, y) and which wants f(y, t)? Which output is transposed?

solve_ivp — time first.
```
sol = solve_ivp(lambda t, y: -2*y, [0,5], [3.0], t_eval=[1.0], rtol=1e-10)
sol.y[0][0]   # value at t=1, rows = states
```
Why this step? Pitfall P1 (see the definition box up top): "solve_ivp = time super first" → f(t, y). Output sol.y has shape (n_states, n_times).
odeint — state first, old = reversed.
```
y = odeint(lambda y, t: -2*y, 3.0, [0.0, 1.0])
y[1][0]        # value at t=1, rows = TIMES here
```
Why this step? odeint is the legacy LSODA wrapper: f(y, t) and output shape (n_times, n_states) — both orderings are transposed relative to solve_ivp.
Both must give the same physics. Exact $y (1) = 3 e^{- 2} = 0.4060$ . Why this step? The trap is purely about interface, not answer — get the argument order right and both return $0.406$ .

Flashcards

Does quad compute signed or geometric area over

[0, 2 π]

sin

Signed — the answer is

0

because the lobes cancel; geometric area would be

4

What does quad(f, 3, 3) return?

(0.0, 0.0)

— zero-width interval, zero area.

What does quad(f, 5, 1) give relative to quad(f, 1, 5)?

The negative of it — reversed limits flip the sign.

How do you make quad handle a singularity inside the interval at

x = c

Pass points=[c] (or split the integral at

c

) so it treats

c

as an endpoint of each piece.

What is

\int_{- \infty}^{\infty} e^{- x^{2}} d x

and how do you call it?

π \approx 1.77245

; quad(f, -np.inf, np.inf).

Which diagnostics show RK45 grinding on a stiff problem?

sol.nfev (RHS evaluations) and sol.t.size (accepted steps) balloon vs BDF.

Newton cooling: temperature after one time-constant?

The gap to room temperature has shrunk by factor

e

(here from

70

\approx 25.75