The Tangent Line Approximation
The Geometric Idea
An ordinary differential equation dy/dx = f(x,y) assigns a slope to every point in the (x,y) plane — a direction field. The true solution y(x) is a curve that everywhere follows those assigned slopes.
Euler's method converts the continuous direction field into a discrete walk:
> (xₙ, yₙ) → (xₙ + h, yₙ + h·f(xₙ, yₙ))
From point (xₙ, yₙ), move distance h along the tangent line. Arrive at an approximate next point. Repeat.
Geometric error: the tangent line at (xₙ, yₙ) has slope f(xₙ, yₙ), but the true curve has a different slope at every point along the interval [xₙ, xₙ + h]. The Euler step uses the slope at the left endpoint throughout — 'the slope that was.' The error per step grows as h².
Accumulated Error
Over N steps to reach a fixed endpoint x = a, with h = a/N:
- Local truncation error per step: O(h²)
- Number of steps: N = a/h
- Global error: O(h²) × (a/h) = O(h) — first-order accuracy
Euler's method is first-order: halving h halves the global error.
Running Euler's Method
Consider dy/dx = y, with initial condition y(0) = 1. True solution: y(x) = eˣ, so y(1) = e ≈ 2.71828.
Apply Euler's method with h = 0.5, from x = 0 to x = 1 (2 steps):
Step 1: y₁ = y₀ + h·f(x₀, y₀) = 1 + 0.5·(1) = 1.5. New point: (0.5, 1.5).
Step 2: y₂ = y₁ + h·f(x₁, y₁) = 1.5 + 0.5·(1.5) = 2.25. New point: (1.0, 2.25).
Euler gives 2.25 vs true value 2.71828. Error: 0.468. Relative error: ~17%.
Deriving Euler's Stability Region
For the test equation dy/dx = λy (where λ is a complex number), Euler's method gives:
> yₙ₊₁ = yₙ + h·λ·yₙ = yₙ·(1 + hλ)
The amplification factor per step: z = 1 + hλ.
Stability condition: the computed solution stays bounded if and only if |z| ≤ 1, i.e., |1 + hλ| ≤ 1.
This is a geometric condition in the complex hλ-plane: the point hλ must lie inside the circle of radius 1 centered at (-1, 0).
Euler's stability region: { hλ ∈ ℂ : |1 + hλ| ≤ 1 }
For real, negative λ (a decaying ODE like dy/dx = -2y): hλ must lie in the interval (-2, 0] on the real axis. With λ = -2 and h = 0.5: hλ = -1. This is exactly on the stability boundary — the method is marginally stable, which explains the qualitative failure in the earlier example.
With h = 1 and λ = -2: hλ = -2, putting us outside the stability region. The solution oscillates with growing amplitude.
Finding the Stability Boundary
Runge-Kutta 4 (RK4) has a larger stability region than Euler, which is one reason it is preferred for most problems.
For real negative λ, RK4 allows hλ down to approximately -2.785 on the real axis (vs Euler's -2 limit).
For stiff equations with eigenvalues λ at very different magnitudes — say λ₁ = -1 and λ₂ = -1000 — stability requires hλ₂ to stay inside the region. For RK4 on the real axis: h·(-1000) ≥ -2.785, so h ≤ 0.002785.
This tiny step size, dictated by the stiff eigenvalue λ₂, makes the simulation expensive even though the slow component λ₁ could use h = 2.
Fixed Points & Attraction Basins
Euler's method applied to dy/dx = f(y) defines a discrete map: yₙ₊₁ = g(yₙ) = yₙ + h·f(yₙ).
A fixed point of this map: y such that g(y) = y. For Euler on dy/dx = f(y), fixed points satisfy f(y) = 0 — equilibria of the ODE.
Stability of a fixed point: if |g'(y)| < 1, nearby iterates converge to y. If |g'(y*)| > 1, they diverge.
g'(y) = 1 + h·f'(y). At a fixed point y: |1 + h·f'(y)| < 1 for stability.
This is exactly the Euler stability condition with λ = f'(y*) — the linearization of the ODE at the equilibrium.
Attraction basin: the set of initial conditions that converge to y* under the Euler map. For nonlinear systems, the basin boundary defines where the simulation will reliably track the ODE equilibrium vs diverge to another attractor.
The simulation loop is a discrete dynamical system. Its qualitative behavior — convergence, oscillation, divergence — depends on the step size h relative to the geometry of the ODE's direction field.
Connecting Geometry to Simulation Design
The geometry of numerical simulation comes down to three questions:
1. Where is the stability region? For Euler: the disk |1 + hλ| ≤ 1. Larger for RK4. Unbounded (entire left half-plane) for implicit methods.
2. Where are the ODE's eigenvalues? The eigenvalues λ of the Jacobian of f at each point determine which stability region must contain hλ.
3. What h keeps all hλ inside the region? The maximum allowable h = (stability region radius) / max|λ|.
For stiff systems: max|λ| is huge, forcing tiny h for explicit methods. Implicit methods are expensive per step but allow large h — they trade per-step cost for stability.
Hamming's insight translates: the choice of numerical method encodes a bet about the geometry of the ODE's eigenvalue spectrum. Making that bet explicit is the first design decision in any simulation.