🛰️Ordinary Differential Equations: Initial Value Problems

Euler's Method and Error

Initial value problems, Euler's method, local and global error, order, consistency, stability and convergence

Cape Canaveral, February 1962

John Glenn is days from becoming the first American to orbit the Earth. The trajectory of Friendship 7 — launch, three orbits, re-entry, splashdown — has been computed on NASA's new IBM 7090 electronic computers. They are fast and, for the time, astonishingly capable. Glenn does not entirely trust them. He makes a request that has become famous: have Katherine Johnson, a mathematician at NASA's Langley Research Center, check the numbers by hand. If she says they are good, he says, he is ready to go. Johnson works through the same equations of motion on a desk calculator, step by step, for a day and a half. Her results agree with the machine's. On 20 February 1962 Glenn flies three orbits of about 88.5 minutes each and comes home. What both Johnson and the IBM were doing is the subject of this unit: taking a law that says how a state changes — Newton's second law for a spacecraft — and turning it, step by step, into where the spacecraft is.

An initial value problem asks for

y(t)

given

y'(t) = f(t, y), \qquad y(t_0) = y_0,

where

y

may be a vector — position and velocity of a spacecraft, temperatures and concentrations in a reactor, currents and voltages in a circuit. If

f

is continuous and Lipschitz in

y

, a unique solution exists (the Picard–Lindelöf theorem). Only a handful of such problems have formula solutions; for the rest we compute

y

at discrete times

t_n = t_0 + nh

y_{n+1} = y_n + h\,f(t_n, y_n)

Euler's method (1768): follow the tangent line from the current point for one step h. It is the first-order Taylor expansion with the h² term dropped.

💡Local error, global error and order

Each Euler step makes a local truncation error — the error of one step started from the exact solution — of size

\tfrac{h^2}{2}y''

. Reaching a fixed time

T

takes

T/h

steps, so the local errors accumulate into a global error of order

h

|y_N - y(T)| \le C\,h.

A method whose global error is

O(h^p)

has order

p

. Euler is first order: on Newton's law of cooling, halving

h

from 1 to 0.0625 min halves the error each time (ratios 1.95, 1.98, 1.99, 2.00). Ten times the accuracy costs ten times the steps.

Derivation

From local errors to the global error

Let

e_n = y_n - y(t_n)

. By Taylor's theorem the exact solution satisfies Euler's formula up to the local error

\tau_n = \tfrac{h^2}{2}y''(\xi_n)

y(t_{n+1}) = y(t_n) + h\,f\big(t_n, y(t_n)\big) + \tau_n.

Subtract this from Euler's step

y_{n+1} = y_n + h f(t_n, y_n)

e_{n+1} = e_n + h\Big[f(t_n, y_n) - f\big(t_n, y(t_n)\big)\Big] - \tau_n.

The bracket is the slope evaluated at the wrong point. If

f

is Lipschitz with constant

L

(its slope in

y

is at most

L

), the bracket is at most

L|e_n|

in size. With

|\tau_n| \le \tau = Mh^2/2

, where

M

bounds

|y''|

|e_{n+1}| \le (1 + hL)\,|e_n| + \tau.

Start from

e_0 = 0

and unroll:

|e_N| \le \tau\big(1 + (1 + hL) + \dots + (1 + hL)^{N-1}\big) = \tau\,\dfrac{(1 + hL)^N - 1}{hL}

. Since

1 + hL \le e^{hL}

and

Nh = T

|e_N| \le \frac{Mh^2}{2}\cdot\frac{e^{LT} - 1}{hL} = \frac{M\,(e^{LT} - 1)}{2L}\;h.

First order, as claimed. The derivation shows both halves of the next box at work: the

h^2

\tau

is consistency, and the bounded growth factor

(1 + hL)^N \le e^{LT}

— earlier errors are amplified, but by a bounded amount — is stability.

💡Consistency + stability = convergence

A method converges (global error → 0 as

h \to 0

) when two things hold:

consistency — the local error per step vanishes faster than $h$ , so the method approximates the right differential equation;
stability — errors made in early steps are not amplified without bound by later ones.

This 'consistency plus stability equals convergence' principle (Dahlquist, 1956, for ODEs; Lax, for PDEs) runs through the rest of the course. Consistency is usually easy to check; stability is where methods fail.

I like to work problems. If you do your best, nobody can ask you to do it over again.

– Katherine Johnson, interview with AARP (2018)

📐Euler around the Earth

A circular orbit at 200 km altitude — period 88.49 minutes, like Friendship 7's — integrated for one revolution with Euler's method. The spacecraft should return to its starting point. It misses by 4,213 km with 10 s steps, 463 km with 1 s steps, 47 km with 0.1 s steps — first order, and hopeless: even 53,000 steps leave a 47 km error. Worse, the error has a direction: the orbit gains 9.8 km of altitude per revolution. Each Euler step moves along the tangent, which always points slightly outside the circle, so the computed spacecraft spirals outward and gains energy it never had.

Energy drift. For the oscillator

x'' = -x

, explicit Euler multiplies the energy by exactly

(1 + h^2)

every step: after 2000 steps of 0.05 it has grown by a factor of 147. Change one thing — update the velocity first, then use the new velocity to update the position (semi-implicit or symplectic Euler) — and the energy stays within 2% of its true value for ever, at no extra cost. Same order, completely different long-term behaviour. Lesson 9-5 explains why.

Derivation

Why explicit Euler's energy grows by exactly

1 + h^2

Write the oscillator as

x' = v

v' = -x

, with energy

E = \tfrac12(x^2 + v^2)

. One explicit Euler step uses the old values on both right-hand sides:

x_{n+1} = x_n + hv_n, \qquad v_{n+1} = v_n - hx_n.

Square and add. The cross terms

2hx_nv_n

and

-2hx_nv_n

cancel:

x_{n+1}^2 + v_{n+1}^2 = x_n^2 + h^2v_n^2 + v_n^2 + h^2x_n^2 = (1 + h^2)\,(x_n^2 + v_n^2).

Every step multiplies the energy by

1 + h^2

, whatever the state: after

n

steps by

(1 + h^2)^n

. With

h = 0.05

and 2000 steps,

1.0025^{2000} = 147.5

— the lab's growth factor. The orbit of the cold-open example spirals outwards for the same reason.

💡Where you will use this

Orbital Mechanics and Rocket trajectories: every propagator is an IVP solver. Attitude Control: the attitude and angular-velocity equations of a spacecraft, stepped at the flight computer's control rate. Circuits and Power Electronics: transient simulation. Structural dynamics: time-history analysis under earthquake records. In real-time flight software Euler-type steps still appear — because they are cheap — and their errors must be budgeted.

📚Further reading

Arieh Iserles, A First Course in the Numerical Analysis of Differential Equations (2nd ed., Cambridge University Press, 2009), chapter 1, 'Euler's method and beyond' — the convergence proof above, carefully.
Ernst Hairer, Syvert Nørsett and Gerhard Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems (2nd ed., Springer, 1993) — the standard reference for everything in this unit except stiffness.
Michael Heath, Scientific Computing: An Introductory Survey (2nd ed.; SIAM Classics 2018), chapter 9, 'Initial Value Problems for Ordinary Differential Equations'.

↓But there is a problem

Euler is too inaccurate for anything that matters: first order means a hundred-fold effort for two more digits. The obvious cure is to use more information per step — not just the slope at the start, but slopes sampled across the step. Carl Runge (1895), Karl Heun (1900) and Wilhelm Kutta (1901) worked out how. Eighty years later two numerical analysts in England, John Dormand and Peter Prince, found a way to make such a method check its own accuracy at every step — and their 1980 method is what runs when you call scipy.integrate.solve_ivp or MATLAB's ode45.