5.4.9 · D2Scientific Computing (Python)

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

2,123 words10 min readBack to topic

Before we start, three words we will use constantly. Nothing else is assumed.

We write the whole problem as

  • — read "the rate at which changes as changes." It is the slope of the curve.
  • — the rule that hands you that slope, given where and when you are.
  • — the starting dot: at the start time we are told exactly where we are.

Our whole job: knowing only the arrows and the starting dot, draw the curve.


Step 1 — The exact truth we are allowed to trust

WHAT. We do not guess. We start from a statement that is perfectly true, the Fundamental Theorem of Calculus:

  • — where we are now (the -th stop on our march).
  • — the step size: a small chunk of time we jump forward. Picture the width of one stride.
  • — our height now.
  • — our height after the stride. This is the unknown we want.
  • — the total change piled up during the stride. It is the area under the slope curve over that little time window.

WHY this equation. It says something obvious once you see it: your new height = old height + (all the change that happened in between). And "all the change" is literally the area under the slope. Nothing is approximate yet.

PICTURE. The slope curve over one stride; the shaded area is exactly the amount climbs.

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

Step 2 — The laziest estimate: Euler (one flat guess)

WHAT. We approximate the area by the simplest shape: a rectangle. Its height is the slope at the start of the step, , and its width is . Area . So:

  • — height now (shorthand for ).
  • — the slope arrow read off right where we stand.
  • — "speed × time = distance travelled," the rectangle's area.
  • — our estimate of the next height.

WHY this tool and not another. A rectangle is the crudest Riemann rectangle: assume the slope stays frozen for the whole stride. It costs one slope evaluation — the cheapest possible.

PICTURE. We walk in a straight line along the starting arrow. If the true curve bends away, we drift off it.

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

Step 3 — Peek at the middle (why one extra sample helps so much)

WHAT. The rectangle used the slope only at the left edge. What if we instead use the slope at the middle of the step? First take a half-Euler step to reach the midpoint, read the slope there, and use that slope for the full stride.

  • — slope at the start (same arrow as Euler).
  • — the midpoint in time.
  • — a cheap guess of the height at the midpoint, made by half-stepping along .
  • — the slope arrow read at that guessed midpoint.

Then step the whole width using : (this is the "midpoint method").

WHY the midpoint. For a curve that bends steadily, the average slope over the step is best represented by the slope halfway through — the left edge is too shallow, the right edge too steep, the middle splits the difference. This is the same instinct that made Simpson's rule fit a parabola instead of a flat top: sample where the shape is best summarised.

PICTURE. Grey dashed = the throwaway half-step to find where the middle is; red = the actual full stride taken along the midpoint slope. It hugs the true curve far tighter than Euler.

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

Step 4 — Sample four slopes and blend them: RK4 is born

WHAT. Two samples are good; four are excellent. Runge–Kutta of order 4 evaluates the slope at four trial points inside the step and blends them:

  • The four 's are four different arrows, probed at start / middle / middle / end.
  • — a weighted average slope. The middle samples count twice as much (weights vs ).

WHY these exact weights . Look closely: they are Simpson's rule (, since ). We are literally applying the best 3-point area rule to the little integral of Step 1. That is why RK4 nails the area so accurately: its error shrinks like per step — halve and the error drops by a factor of 32.

PICTURE. All four sampled arrows drawn where they are read; the final blended arrow lands almost exactly on the true curve.

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

Step 5 — How the solver picks by itself (the "45" in RK45)

WHAT. A fixed is wasteful on smooth stretches and dangerous on sharp ones. RK45 (Dormand–Prince, solve_ivp's default) uses its slope samples to compute two answers at once: a 5th-order estimate and a 4th-order estimate , from the same 's. Their gap estimates how wrong the step is.

  • If is too big → the step was risky → shrink and redo.
  • If is comfortably small → the road is smooth → grow to save work.

WHY compare two orders. You cannot know your error without knowing the true answer — unless you have a second, independent estimate to compare against. The difference between a good guess and a better guess is a reliable stand-in for "how far off am I?" This is the embedded-pair trick.

PICTURE. Big lazy strides on the flat part, tiny cautious strides through the sharp bend — same target accuracy everywhere.

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

Step 6 — The degenerate cases (never leave the reader stranded)

WHAT & WHY, each case shown:

Case A — flat slope, . Every arrow points horizontally. All 's are , so . The curve is a flat line; the solver correctly does nothing but coast, and can take a giant single step.

Case B — constant slope, . All four 's equal , blend . Even lowly Euler is exact here — a straight line has no curvature to miss. No method beats another on a straight road.

Case C — stiff decay, with huge . The true curve plummets then flattens. An explicit stride overshoots below zero and bounces bigger each step → the numbers explode even though the true answer is calmly near zero. This is the instability that forces implicit methods.

PICTURE. The three regimes side by side: flat (trivial), straight (exact), stiff (explicit blows up while implicit stays glued).

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

The one-picture summary

Everything above is one idea repeated: estimate the shaded area of Step 1, take a stride, repeat, and shrink the stride only where the road curves.

Figure — scipy.integrate — odeint, solve_ivp (RK45, DOP853), quad
Recall Feynman retelling — the whole walkthrough in plain words

You know how fast something is moving at every point on the map, but not the actual path it traces. To find the path you cheat cleverly. Step 1: the exact rule is "new position = old position + all the distance covered in this little moment," and that distance is the area under the speed curve — but you can't measure that area directly, so you estimate it. Step 2 (Euler): the dumbest estimate — assume your speed never changes for the whole stride and walk in a straight line. It drifts off on any bend. Step 3: peek at your speed halfway through the stride instead — much better, because the middle speed represents the whole stride fairly. Step 4 (RK4): peek four times — start, middle, middle, end — and take a weighted average, counting the middles double. Those weights are exactly Simpson's area rule, which is why it's so accurate. Step 5 (RK45): do two quality guesses at once; if they disagree a lot the road is tricky so take baby steps, if they agree the road is smooth so take giant steps. Step 6: on flat or perfectly straight roads any method is trivially right, but on "stiff" roads the explicit stride bounces and explodes unless you switch to a stabler implicit solver. That's the entire family of ODE solvers — Euler, RK4, RK45, DOP853, BDF — all born from estimating one little shaded area.

Recall

Why is RK4's blend so accurate? ::: The weights are Simpson's rule applied to the little integral in Step 1 — it fits a parabola to the slope, exact through cubics. What does the RK45 error estimate actually compare? ::: A 5th-order and an embedded 4th-order answer built from the same slope samples; their difference is the local error used to adapt . On (constant slope), which method is exact? ::: All of them, even Euler — a straight line has no curvature to miss. Why does explicit RK blow up on a stiff decay ? ::: Stability, not accuracy: the stride must be microscopic to avoid the step-to-step bounce amplifying; switch to implicit BDF/Radau.