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

quad. You have the function and want a number (an area) — that is the definition of a definite integral. Call:

from scipy.integrate import quad
value, error = quad(f, 0, 1)

It returns a tuple: value = the estimate, error = an estimate of the absolute error. If you only see one number printed, you forgot to unpack the tuple.

• solve_ivp wants time first: f(t, y) → so rhs_A(t, y).
• odeint is the old interface (wraps LSODA) and is reversed: f(y, t) → so rhs_B(y, t).

Mnemonic from the parent: "solve_ivp = time is super first." Old = reversed.

• (a) quad — a fixed area.
• (b) solve_ivp — you know the slope law $\ddot{\theta }=-\frac{g}{L}\sin \theta$ and an initial state, want the path.
• (c) quad — you HAVE $i(t)$ and want the accumulated number.
• (d) solve_ivp — acceleration is the slope of velocity; integrate the path from initial conditions.

Rule of thumb: "Do I already know the whole curve and want its area?" → quad. "Do I only know the rate of change and a start?" → solve_ivp.

Antiderivative is $x^3+x$, so exact value $=(8+2)-(0)=10$.

from scipy.integrate import quad
val, err = quad(lambda x: 3*x**2 + 1, 0, 2)
# val ≈ 10.0, err ≈ 1.1e-13

quad is exact for polynomials up to high degree (Gauss quadrature with $n$ nodes nails degree $2n-1$), so a cubic integrand comes back essentially at machine precision — see Gaussian quadrature.

$y(5)=3e^{-10}\approx 3\times 4.5400\times 10^{-5}\approx 0.00013620$.

from scipy.integrate import solve_ivp
import numpy as np
sol = solve_ivp(lambda t, y: -2*y, [0, 5], [3.0],
                rtol=1e-9, atol=1e-12)
sol.y[0, -1]   # ≈ 0.00013620

• y0=[3.0] is a list because the state is always a vector, even a 1-element one.
• sol.y[0] picks the first (only) state component; [-1] is the last time — the value at $t=5$.

from scipy.integrate import odeint
import numpy as np
t = np.linspace(0, 5, 200)
y = odeint(lambda y, t: -2*y, 3.0, t)   # note (y, t)!
y[-1, 0]   # ≈ 0.00013620

Two things flip vs solve_ivp: the RHS is f(y, t), and the output shape is (len(t), n_states) — so the final value is y[-1, 0] (last row, first column), not y[0, -1].

Simpson on $[a,b]$ with midpoint $m$ and half-width $h=\frac{b-a}{2}=1$: ${\int }_a^{b}f\approx \frac{h}{3}[f(a)+4f(m)+f(b)]=\frac{1}{3}[0+4(1)+8]=\frac{12}{3}=4.$ Exact: ${\int }_0^2{x}^{3}dx=\frac{x^4}{4}{|}_0^2=\frac{16}{4}=4$. Zero error. Why exact? Simpson fits a parabola (degree 2) through 3 points, yet the error term is $-\frac{h^5}{90}{f}^{(4)}(\xi )$. For a cubic, $f^{(4)}=0$ everywhere — so Simpson is secretly exact for cubics too. This is the same "free extra order" that makes Gauss rules so cheap (see Riemann sums and numpy.trapz and simpson).

Neither a bug nor real physics — it's accumulated numerical error. Each step has a tiny error; over thousands of oscillations these add up, and with a loose rtol=1e-3 the solver is allowed large local slips that can slowly pump energy in or out. The true solution conserves amplitude (see Harmonic oscillator). Fix: tighten tolerances (rtol=1e-10, atol=1e-12) or use a higher-order method (DOP853). The figure shows loose vs tight tolerance over many periods.

Because t_eval only chooses where the result is reported, not how the solver steps. Internally solve_ivp takes as many adaptive steps as accuracy demands (governed by rtol/atol) and then interpolates onto your t_eval points. Accuracy is set by tolerances and method — never by t_eval. Reporting at 5 points vs 500 points gives the same underlying integration.

Why reduce? Every solver accepts only first-order systems. Introduce state vector $y=(y_0,y_1)=(x,\dot{x})$: ${\dot{y}}_0=y_1,{\dot{y}}_1=-{\omega }^2{y}_0.$ Exact solution $x(t)=\cos (\omega t)=\cos (2t)$, so $x(\pi )=\cos (2\pi )=1$.

from scipy.integrate import solve_ivp
import numpy as np
w = 2.0
sol = solve_ivp(lambda t, y: [y[1], -w**2*y[0]],
                [0, np.pi], [1.0, 0.0],
                method='DOP853', rtol=1e-11, atol=1e-13)
sol.y[0, -1]   # ≈ 1.0

At $t=0$: $\dot{x}=0,x=1\Rightarrow E=\frac{1}{2}(0)+\frac{1}{2}(4)(1)=2$. At $t=\pi$ the true state is $x=1,\dot{x}=0$ again → $E=2$. Numerically with DOP853 at tight tolerance, $E$ stays $\approx 2$ to $\sim 10^{-9}$. Energy drift is a direct read-out of solver error, which is why oscillators are the classic accuracy stress-test.

Charge as one number — this is an area, use quad: $Q={\int }_0^{\infty }{e}^{-t}\sin tdt=\frac{1}{1^2+1^2}=\frac{1}{2}=\mathrm{0.5.}$ (Standard result ${\int }_0^{\infty }{e}^{-at}\sin (bt)dt=\frac{b}{a^2+b^2}$ with $a=b=1$.)

from scipy.integrate import quad
import numpy as np
Q, err = quad(lambda t: np.exp(-t)*np.sin(t), 0, np.inf)  # Q ≈ 0.5

Charge over time — this is a path, use solve_ivp:

from scipy.integrate import solve_ivp
sol = solve_ivp(lambda t, q: np.exp(-t)*np.sin(t), [0, 20], [0.0],
                rtol=1e-10, atol=1e-12)
sol.y[0, -1]   # ≈ 0.5  (q(∞) matches Q)

The two answers must agree: the ODE's final value is the definite integral. This is the Fundamental Theorem of Calculus made numerical — quad accumulates area directly, solve_ivp accumulates it as a running total.

Why it crawls: the $-1000y$ term is a very fast decay (time-scale $1/1000$ s) sitting alongside the slow $e^{-t}$ drive (time-scale $1$ s). That coexistence of fast + slow scales is the definition of stiffness (see Stiff differential equations). An explicit method like RK45 is only stable if its step $h$ resolves the fastest scale: roughly $h\lesssim 2/1000$. Even once the fast transient has died and the solution is smooth, RK45 is forced to keep tiny steps just to avoid blowing up — this is a Numerical stability limit, not an accuracy limit. Fix: use an implicit solver that is stable at large $h$:

from scipy.integrate import solve_ivp
sol = solve_ivp(lambda t, y: -1000*y + 3000 - 2000*__import__('numpy').exp(-t),
                [0, 5], [0.0], method='BDF', rtol=1e-8, atol=1e-10)

BDF, Radau, or LSODA take big smooth steps here. See Runge-Kutta methods for why explicit RK cannot.

Exact: ${\int }_0^1{x}^{-1/2}dx=[2\sqrt{x}{]}_0^1=2$ — an integrable singularity (the area is finite even though the height isn't).

from scipy.integrate import quad
import numpy as np
val, err = quad(lambda x: 1/np.sqrt(x), 0, 1)  # val ≈ 2.0

quad's adaptive Gauss–Kronrod clusters points near the trouble spot and returns $\approx 2$ with a small err. Warning signals: a large returned err, or an IntegrationWarning. For a known kink or singularity, pass points=[...] (or split the integral) so the sampler is told where to be careful.

Use an event function: a function g(t, y) whose sign change the solver watches for, locating the crossing by root-finding between steps.

def hit_zero(t, y): return y[0]
hit_zero.terminal = True     # stop at the crossing (optional)
sol = solve_ivp(rhs, [0, 10], y0, events=hit_zero, rtol=1e-9)
sol.t_events[0]              # exact crossing time(s)

Why not t_eval: t_eval only samples the solution at pre-chosen times — the true crossing almost never lands on a grid point. Events do a root solve to pinpoint the moment to full precision. This is one of the concrete reasons the parent recommends solve_ivp over odeint.

quad returns a tuple ::: (value, estimated_absolute_error) solve_ivp RHS order ::: f(t, y) — time first odeint RHS order and output shape ::: f(y, t); output (len(t), n_states) Cure for a stiff ODE ::: switch to an implicit method (BDF/Radau/LSODA) What t_eval controls ::: where the answer is reported, never the internal step size