Numerical Solution of Differential Equations

Orientation

When there is no formula to find

Four units have solved differential equations exactly: separation and integrating factors, the characteristic equation, and eigenvalues for systems. Each produced a formula.

Those equations are a thin slice. Change $y^{'} = 2 x y$ to

y^{'} = t^{2} + y^{2},

and no elementary formula exists. Nothing is wrong with the equation. The right side is a polynomial, smooth everywhere, and the solution through $(0, 0)$ is a perfectly ordinary curve. It simply cannot be written with the functions we have names for. That is the normal case; the solvable equations are the exception.

What remains available is the equation itself. It reports the slope at any point, and that is enough to walk the curve out numerically: from a known point, follow the slope a short distance, ask the equation again at the new point, and repeat. That is Euler's method, and it turns a differential equation into arithmetic.

The answer changes character. Instead of $y = 3 e^{x^{2}}$ there is a list of numbers, $y$ at $t = 0.1$ , at $0.2$ , and so on. Applied to the equation above, this yields $y (1) \approx 0.3502318443$ , a value no formula in the corpus can produce.

Two questions follow.

How wrong is it? Every step follows a slope that is exactly right only at its start, while the true solution curves away during the step. Errors accumulate. Halving Euler's step only halves the error, so a tenfold increase in work buys about one extra digit, poor value. Sampling the slope more than once per step does far better: two evaluations quarter the error, four divide it by sixteen.

When does it break? Not always where one would expect. Applied to a rapidly decaying equation with too large a step, Euler returns values that oscillate and grow. One case in this unit produces $- 243$ where the true answer is about $2 \times 10^{- 9}$ . That is the method failing, not the solution misbehaving, and there is an exact condition that says when it will happen.

Why this matters

Why one slope per step is not enough

Euler's method advances each step using the slope at the step's start. The following figures give its error behaviour before any refinement is introduced.

Take $y^{'} = y$ with $y (0) = 1$ , whose exact solution is $y = e^{t}$ , so $y (1) = 2.7182818285$ . With $h = 0.25$ the method takes four steps, each multiplying by $1 + h = 1.25$ :

1 \to 1.25 \to 1.5625 \to 1.953125 \to 2.441406 .

The answer is $2.441406$ , low by $0.276876$ , about 10% off. Every step used the slope at its own start, and since this solution curves upward, each tangent line fell below the curve and the shortfalls accumulated.

Refining the step helps slowly. Halving $h$ repeatedly and recording the error at $t = 1$ :

$n$ steps	$h$	Euler's $y (1)$	Error	Ratio to previous
1	$1$	$2.00000000$	$0.71828183$	—
2	$0.5$	$2.25000000$	$0.46828183$	$1.53$
4	$0.25$	$2.44140625$	$0.27687558$	$1.69$
8	$0.125$	$2.56578451$	$0.15249731$	$1.82$
16	$0.0625$	$2.63792850$	$0.08035333$	$1.90$
32	$0.03125$	$2.67699013$	$0.04129170$	$1.95$
64	$0.015625$	$2.69734495$	$0.02093688$	$1.97$
128	$0.0078125$	$2.70773902$	$0.01054281$	$1.99$

The ratios climb toward 2: halving the step halves the error. This is what first order means. Going from 1 step to 128, which is 128 times the work, reduced the error from $0.72$ to $0.011$ , about two digits.

Notice also that the ratio approaches 2 rather than equalling it. Order is a statement about the trend as $h \to 0$ , and at coarse steps the observed ratio can be well short of it.

The remedy is to look at the slope more than once. The slope at the start of a step is not the slope across it. Evaluating $f$ again at the far end, using an Euler step to guess where that end is, and averaging the two slopes is the improved Euler method. The averaging cancels the leading error term, and the error now falls with $h^{2}$ :

$n$	$h$	Improved Euler $y (1)$	Error	Ratio
1	$1$	$2.5000000000$	$2.183 \times 10^{- 1}$	—
2	$0.5$	$2.6406250000$	$7.766 \times 10^{- 2}$	$2.81$
4	$0.25$	$2.6948556900$	$2.343 \times 10^{- 2}$	$3.32$
8	$0.125$	$2.7118412386$	$6.441 \times 10^{- 3}$	$3.64$
16	$0.0625$	$2.7165935225$	$1.688 \times 10^{- 3}$	$3.81$
32	$0.03125$	$2.7178496740$	$4.322 \times 10^{- 4}$	$3.91$

The ratios head for 4. At $n = 32$ improved Euler has error $4.3 \times 10^{- 4}$ against Euler's $4.1 \times 10^{- 2}$ . A hundredfold improvement for twice the work per step.

Four evaluations do better still. RK4 reaches error $3.3 \times 10^{- 7}$ at $n = 16$ , with ratios climbing toward 16:

$n$	RK4 $y (1)$	Error	Ratio
1	$2.708333333333$	$9.948 \times 10^{- 3}$	—
2	$2.717346191406$	$9.356 \times 10^{- 4}$	$10.6$
4	$2.718209939201$	$7.189 \times 10^{- 5}$	$13.0$
8	$2.718276844417$	$4.984 \times 10^{- 6}$	$14.4$
16	$2.718281500341$	$3.281 \times 10^{- 7}$	$15.2$

With a single step, RK4 is more accurate than Euler with 128. On this problem, four evaluations per step buy more accuracy than the same total work spent on additional first-order steps.

This comparison uses one smooth, non-stiff equation with a fixed step. It establishes the order of each method on that problem. It does not establish which method to use in general: stiff equations, adaptive step-size control, and error tolerances change the comparison, and solver libraries select among several methods on those grounds.

Definition

Reading the formulas: where each slope is evaluated

Almost every error in applying these formulas is an error about where a slope is evaluated.

Euler: one slope, at the current point. In $y_{n + 1} = y_{n} + h f (t_{n}, y_{n})$ both arguments are current values. The formula is the tangent line at $(t_{n}, y_{n})$ , extended a distance $h$ . It uses the equation once per step and nothing else.

**Improved Euler: the second slope is at a predicted point.** The formula

k_{2} = f (t_{n} + h, y_{n} + h k_{1})

evaluates $f$ at the far end of the step, but at what $y$ ? Not $y_{n + 1}$ , which is not known yet; at $y_{n} + h k_{1}$ , the Euler prediction of where the far end lies. This is the step most often got wrong, by evaluating $k_{2}$ at $y_{n}$ instead. The method is a predictor followed by a corrector: guess the endpoint with Euler, then re-estimate the slope there and average.

The averaging is what earns the second order. Overestimating the slope at one end tends to be compensated by the other, and the leading error term cancels.

RK4: two of the four slopes are at the midpoint. Reading the weights $\frac{h}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4})$ : the coefficients sum to $1 + 2 + 2 + 1 = 6$ , matching the divisor, so the bracket is a weighted average of slopes and the formula still has the shape "value plus step times a slope". The midpoint slopes $k_{2}$ and $k_{3}$ carry double weight, because the average slope across an interval is better represented by the middle than by the ends.

Note that $k_{3}$ uses $k_{2}$ , not $k_{1}$ : each slope is evaluated at a point built from the previous one. The four evaluations are sequential, not independent.

Order: what $O (h^{p})$ claims. Saying the error at a fixed final time is $O (h^{p})$ means it behaves like $C h^{p}$ for small $h$ , so replacing $h$ by $h / 2$ multiplies it by $2^{- p}$ . Three consequences:

It is a statement about small $h$ . The observed ratios in this unit climb toward their limits, Euler's reach $1.99$ at $n = 128$ , RK4's $15.2$ at $n = 16$ , rather than attaining them.
It says nothing about the constant $C$ . A higher-order method can be less accurate at coarse steps, though not usually.
It counts steps, not work. Halving $h$ doubles the number of steps, so comparisons between methods must be made per slope evaluation, which is what makes RK4's factor of 16 worth its four evaluations.

Stability: a separate constraint entirely. Applying Euler to $y^{'} = λ y$ gives $y_{n + 1} = (1 + h λ) y_{n}$ , so after $n$ steps $y_{n} = (1 + h λ)^{n} y_{0}$ . A geometric sequence. It decays only when $| 1 + h λ | < 1$ .

For real $λ < 0$ this is $h < 2 / | λ |$ . Above that bound the factor has magnitude greater than 1 and the computed values grow, alternating in sign once $1 + h λ < - 1$ . The true solution decays throughout. Accuracy and stability are different requirements, and a step can be accurate enough in principle yet unstable in practice.

Procedure

Approximating a solution numerically

Input. An initial value problem $y^{'} = f (t, y)$ , $y (t_{0}) = y_{0}$ , a target time $T$ , and an accuracy requirement.

Step 1 — Check whether an exact method applies first. If the equation is separable, linear, or constant-coefficient, solve it exactly; a formula is better than a table of numbers. Numerical methods are for when that fails, which is most of the time, but not always.

Step 2 — Choose a method. Use RK4 unless there is a reason not to: four evaluations per step buy fourth-order accuracy, which almost always costs less total work than Euler at the same accuracy. Euler and improved Euler are worth doing by hand for understanding, and Euler remains useful where $f$ is very expensive and accuracy demands are low.

Step 3 — Check the stability limit if the equation decays rapidly. Estimate $λ = \partial f / \partial y$ along the solution. If $λ$ is negative and large in magnitude, Euler requires

h < \frac{2}{| λ |},

and this bound can be far more restrictive than accuracy alone would suggest. Ignoring it produces output that grows and alternates in sign.

Step 4 — Choose a step size. Set $h = (T - t_{0}) / n$ for a whole number of steps $n$ . If an accuracy is required and a method of order $p$ is in use, an error $E$ at step $h$ predicts an error of roughly $E / 2^{p}$ at $h / 2$ , so the number of halvings needed can be estimated in advance rather than discovered by trial.

Step 5 — Iterate the update. Keep a table of $(t_{n}, y_{n})$ . For each step:

Euler: $y_{n + 1} = y_{n} + h f (t_{n}, y_{n})$ .
Improved Euler: $k_{1} = f (t_{n}, y_{n})$ ; $k_{2} = f (t_{n} + h, y_{n} + h k_{1})$ ; $y_{n + 1} = y_{n} + \frac{h}{2} (k_{1} + k_{2})$ .
RK4: the four slopes in order, each built from the previous, then $y_{n + 1} = y_{n} + \frac{h}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4})$ .

Advance $t_{n + 1} = t_{n} + h$ after each step, and carry more digits than the answer needs.

Step 6 — Estimate the error by halving. With no exact solution to compare against, run the method at $h$ and again at $h / 2$ . The difference between the two answers estimates the error of the coarser one; if the two agree to the digits required, those digits are likely sound. Agreement across three step sizes with the right ratio is stronger evidence still.

Step 7 — Sanity-check the output. Does it respect what the equation says? A solution of $y^{'} = - 20 y$ starting positive must stay positive and decrease; output that alternates sign is instability, not behaviour. Where an exact solution is known, compare directly.

Step 8 — Report it as what it is. Values at chosen points, with an error estimate and the step size used, not a formula, and not exact.

Where it goes wrong.

Evaluating $k_{2}$ at $y_{n}$ instead of at the predicted $y_{n} + h k_{1}$ . This silently reduces improved Euler to something no better than Euler.
Forgetting to advance $t$ , so every step uses $f (t_{0}, y_{n})$ . Harmless when $f$ does not depend on $t$ , wrong otherwise.
Assuming smaller is always better. Below the stability limit refinement helps, but round-off eventually sets a floor and the total error stops falling.
Reading order as a promise at the step in use. Ratios approach $2$ , $4$ and $16$ as $h \to 0$ ; at coarse steps they are lower.
Continuing past a blow-up. For $y^{'} = y^{2}$ with $y (0) = 1$ the solution ceases to exist at $t = 1$ , but the arithmetic will happily return numbers for $t > 1$ . They mean nothing.

Worked example

Four Euler steps, then the same problem improved

Problem. Approximate $y (2)$ for

y^{'} = y - t^{2} + 1, y (0) = 0.5,

using Euler's method with $h = 0.5$ . This equation is linear, so an exact solution exists for comparison: $y = (t + 1)^{2} - \frac{1}{2} e^{t}$ , giving $y (2) = 5.3054719505$ .

Step 1 — Set up. $f (t, y) = y - t^{2} + 1$ , $t_{0} = 0$ , $y_{0} = 0.5$ , and $n = (2 - 0) / 0.5 = 4$ steps.

Step 2 — Iterate $y_{n + 1} = y_{n} + h f (t_{n}, y_{n})$ .

First step, at $(0, 0.5)$ :

f (0, 0.5) = 0.5 - 0^{2} + 1 = 1.5, y_{1} = 0.5 + 0.5 (1.5) = 1.25 .

The exact value is $y (0.5) = 1.425639$ , so this is already low by $0.176$ .

Second step, at $(0.5, 1.25)$ , where both arguments have advanced:

f (0.5, 1.25) = 1.25 - 0.25 + 1 = 2.0, y_{2} = 1.25 + 0.5 (2.0) = 2.25 .

Exact $y (1) = 2.640859$ .

Third step, at $(1, 2.25)$ :

f (1, 2.25) = 2.25 - 1 + 1 = 2.25, y_{3} = 2.25 + 0.5 (2.25) = 3.375 .

Exact $y (1.5) = 4.009155$ .

Fourth step, at $(1.5, 3.375)$ :

f (1.5, 3.375) = 3.375 - 2.25 + 1 = 2.125, y_{4} = 3.375 + 0.5 (2.125) = 4.4375 .

Result. $y (2) \approx 4.4375$ against the exact $5.3054719505$ , an error of $0.86797195$ , about 16%.

$n$	$t_{n}$	$y_{n}$ (Euler)	Exact	Error
0	$0$	$0.500000$	$0.500000$	$0$
1	$0.5$	$1.250000$	$1.425639$	$0.175639$
2	$1.0$	$2.250000$	$2.640859$	$0.390859$
3	$1.5$	$3.375000$	$4.009155$	$0.634155$
4	$2.0$	$4.437500$	$5.305472$	$0.867972$

The error grows at every step: each new error is added to the ones already carried, and the solution's upward curvature means every tangent line falls short.

Step 3 — Refine and watch the ratio. Repeating with smaller steps:

$n$	$h$	Euler $y (2)$	Error
4	$0.5$	$4.43750000$	$0.86797195$
8	$0.25$	$4.77965164$	$0.52582031$
16	$0.125$	$5.01046864$	$0.29500331$
32	$0.0625$	$5.14824995$	$0.15722200$

Successive error ratios are $1.65$ , $1.78$ and $1.88$ , climbing toward 2 as first order predicts. Even at 32 steps the error is still $0.157$ .

Step 4 — Improved Euler on the first step, for contrast. Back at $(0, 0.5)$ with $h = 0.5$ :

k_{1} = f (0, 0.5) = 1.5,

predicted endpoint = y_{0} + h k_{1} = 0.5 + 0.5 (1.5) = 1.25,

k_{2} = f (0.5, 1.25) = 1.25 - 0.25 + 1 = 2.0,

y_{1} = 0.5 + \frac{0.5}{2} (1.5 + 2.0) = 0.5 + 0.25 (3.5) = 1.375 .

Against the exact $1.425639$ the error is $0.050639$ , versus Euler's $0.175639$ : better by a factor of about $3.5$ , for one extra evaluation of $f$ .

Note where $k_{2}$ was taken: at $t = 0.5$ and $y = 1.25$ , the Euler prediction. Evaluating it at $y = 0.5$ instead would give $f (0.5, 0.5) = 1.25$ and the wrong answer $1.1875$ .

Step 5 — Verify what can be verified. The exact solution checks out: $y (0) = (0 + 1)^{2} - 0.5 = 0.5$ , and differentiating $y = (t + 1)^{2} - \frac{1}{2} e^{t}$ gives $y^{'} = 2 (t + 1) - \frac{1}{2} e^{t}$ , while $y - t^{2} + 1 = (t + 1)^{2} - \frac{1}{2} e^{t} - t^{2} + 1 = 2 t + 2 - \frac{1}{2} e^{t}$ , the same expression.

Example

An equation with no closed-form solution

The worked example used an equation with a known exact solution, so the error could be measured. That is the exception. Here is the ordinary case.

Problem. Approximate $y (1)$ for

y^{'} = t^{2} + y^{2}, y (0) = 0 .

Why no formula is coming. The right side is a polynomial in two variables, about as simple as a nonlinear equation gets. But it does not separate: $t^{2} + y^{2}$ is a sum, not a product $f (t) g (y)$ , so the first-order unit's method does not apply. It is not linear in $y$ either, because of the $y^{2}$ , so no integrating factor exists. It has no constant coefficients to feed a characteristic equation. Every exact method in the corpus fails, and in fact no solution exists in elementary functions.

None of that is a defect of the equation. The right side is continuous and differentiable everywhere, so a unique solution through $(0, 0)$ exists on some interval. It simply has no name.

Solve it numerically with RK4. Since there is no exact answer to compare against, run the method at several step counts and look for agreement:

Steps $n$	$h$	RK4 estimate of $y (1)$
10	$0.1$	$0.3502337418$
100	$0.01$	$0.3502318445$
1000	$0.001$	$0.3502318443$
10000	$0.0001$	$0.3502318443$

The estimates stop changing. From $n = 10$ to $n = 100$ the value shifts by about $1.9 \times 10^{- 6}$ ; from $100$ to $1000$ by $2 \times 10^{- 10}$ ; from $1000$ to $10000$ not at all in the digits shown. So

y (1) = 0.3502318443

to ten decimal places, and the last two rows agreeing is the evidence for it.

This is the error-estimation technique of the procedure block put to work: with no exact solution available, agreement between step sizes is the check. The rapid convergence, three extra digits from $n = 10$ to $n = 100$ , is the fourth-order behaviour showing itself.

What has and has not been produced. There is now a number, good to ten digits, for a quantity no formula in this corpus can express. What there is not is a function: to know $y (0.5)$ requires running the method again, and nothing here reveals how $y$ depends on $t$ in general, or how the answer would change with a different starting value. A closed form would give all of that at once, which is why step 1 of the procedure says to look for one first.

A caution about going further. Solutions of nonlinear equations can cease to exist in finite time, $y^{'} = y^{2}$ with $y (0) = 1$ blows up at $t = 1$ , as the first-order unit showed. This equation's solution is well behaved on $[0, 1]$ , but the arithmetic would return numbers past a blow-up just as readily as before one, and those numbers would be meaningless. A numerical method never reports that the solution has ended.

Theorem

Order of accuracy, and where the exponents come from

Theorem (order of accuracy). Let $y$ solve $y^{'} = f (t, y)$ , $y (t_{0}) = y_{0}$ , with $f$ sufficiently smooth on the interval of interest. At a fixed final time $T$ , the error of the approximation computed with step $h$ satisfies

| y (T) - y_{N} | \leq C h^{p}, N h = T - t_{0},

with $p = 1$ for Euler, $p = 2$ for improved Euler and $p = 4$ for RK4. Consequently, replacing $h$ by $h / 2$ multiplies the error bound by $2^{- p}$ .

Where the exponent comes from. Expand the true solution about $t_{n}$ by Taylor's theorem:

y (t_{n} + h) = y (t_{n}) + h y^{'} (t_{n}) + \frac{h^{2}}{2} y^{″} (t_{n}) + \frac{h^{3}}{6} y^{‴} (t_{n}) + \dots

Euler's step reproduces only the first two terms, since $y^{'} (t_{n}) = f (t_{n}, y_{n})$ . The first term it omits is $\frac{h^{2}}{2} y^{″}$ , so the error committed in one step is $O (h^{2})$ , the local error.

Why the global order is one less. Reaching a fixed $T$ takes $N = (T - t_{0}) / h$ steps, so $N$ is proportional to $1 / h$ . Accumulating $N$ local errors of size $O (h^{2})$ gives a total of order

\frac{1}{h} \cdot h^{2} = h,

which is first order. The same trade explains the others: improved Euler's averaging cancels the $h^{2}$ term so its local error is $O (h^{3})$ and its global error $O (h^{2})$ ; RK4's weighted slopes match the expansion through $h^{4}$ , giving local $O (h^{5})$ and global $O (h^{4})$ .

This is why the weights in RK4 are what they are. The coefficients $\frac{1}{6} (1, 2, 2, 1)$ are chosen precisely so that the expansion of the weighted average agrees with the Taylor series to fourth order. They are not arbitrary, and a different weighting would lose the cancellation. The Taylor expansion the argument rests on is the one the series unit develops.

Verified against computation. For $y^{'} = y$ , $y (0) = 1$ on $[0, 1]$ , where $y (1) = e = 2.7182818285$ , the error ratio between successive halvings approaches $2^{p}$ :

$n$	Euler ratio	Improved Euler ratio	RK4 ratio
2	$1.53$	$2.81$	$10.63$
4	$1.69$	$3.32$	$13.01$
8	$1.82$	$3.64$	$14.42$
16	$1.90$	$3.81$	$15.19$
32	$1.95$	$3.91$	—
64	$1.97$	—	—
128	$1.99$	—	—

Predicted limits $2$ , $4$ , $16$ ✓.

Two things the theorem does not say.

It does not promise the ratio at any particular $h$ . Every column above approaches its limit from below. At $n = 2$ RK4's ratio is $10.6$ , not $16$ ; the asymptotic claim needs small $h$ , and how small depends on the constant $C$ and on the size of the derivatives of $f$ .

It does not say the error can be driven to zero. The bound $C h^{p}$ describes truncation error only, meaning the discretisation. Floating-point round-off accumulates in the opposite direction, growing as the number of steps grows, so total error eventually stops improving and then worsens. The theorem's hypotheses also require $f$ smooth along the whole path, which fails if the solution escapes to infinity.

Corollary (choosing a step). If a method of order $p$ gives error $E$ at step $h$ , then attaining a tolerance $ε$ requires about

h_{new} \approx h {(\frac{ε}{E})}^{1 / p} .

For $p = 4$ , a thousandfold improvement needs only about a $5.6$ -fold reduction in $h$ ; for $p = 1$ it needs a thousandfold. That gap is the argument for higher-order methods.

Warning

A smaller step is not always a better answer

Refinement is the natural response to an inaccurate result, and the order table encourages it. The claim that a smaller $h$ is always better fails in two distinct ways.

---

Instability: too large a step can produce nonsense, not merely inaccuracy.

Take $y^{'} = - 20 y$ with $y (0) = 1$ on $[0, 1]$ . The exact solution $y = e^{- 20 t}$ decays fast and monotonically to $y (1) = 2.061153622439 \times 10^{- 9}$ , positive at every point.

Applying Euler with $n = 5$ , that is $h = 0.2$ :

y_{n + 1} = y_{n} + h (- 20 y_{n}) = (1 - 4) y_{n} = - 3 y_{n} .

The computed values are $1, - 3, 9, - 27, 81, - 243$ , and the answer returned for $y (1)$ is

- 2.43 \times 10^{2} .

The true value is about $2 \times 10^{- 9}$ . The output is not slightly wrong; it has the wrong sign, the wrong magnitude by eleven orders, and grows where the solution decays.

Why, exactly. Euler on $y^{'} = λ y$ gives $y_{n + 1} = (1 + h λ) y_{n}$ , so the computed values form a geometric sequence with ratio $1 + h λ$ . Decay requires

| 1 + h λ | < 1 ⟺ h < \frac{2}{| λ |} = \frac{2}{20} = 0.1 .

At $h = 0.2$ the ratio is $- 3$ , whose magnitude exceeds 1, hence growth, and the negative sign explains the alternation. The computed values behave as the formula dictates; the formula has simply stopped tracking the equation.

Crossing the limit:

$n$	$h$	$1 + h λ$	Euler's $y (1)$	Verdict
5	$0.2$	$- 3$	$- 2.430000 \times 10^{2}$	unstable
10	$0.1$	$- 1$	$+ 1.000000$	exactly on the boundary — no decay
20	$0.05$	$0$	$0$	stable
40	$0.025$	$0.5$	$9.094947 \times 10^{- 13}$	stable
100	$0.01$	$0.8$	$2.037036 \times 10^{- 10}$	stable

Note $h = 0.1$ exactly: the ratio is $- 1$ , so values alternate between $1$ and $- 1$ forever and the magnitude never decays. The condition is a strict inequality for a reason.

The practical warning is that the required $h$ depends on how fast the solution decays, not on how accurate one wants to be. An equation with $λ = - 1000$ forces $h < 0.002$ for Euler, however modest the accuracy sought. Equations like this are called stiff, and they are why implicit methods, which have no such limit, exist.

---

Round-off: too small a step stops helping, then hurts.

Truncation error falls as $h^{p}$ , but every step also commits a small floating-point rounding error, and halving $h$ doubles the number of steps over which those accumulate. Total error behaves roughly as

\underset{truncation, falling}{\underset{⏟}{C h^{p}}} + \underset{round-off, rising}{\underset{⏟}{\frac{δ}{h}}},

with $δ$ set by the machine precision. The sum has a minimum at some finite $h$ : below it, more work gives a worse answer. In double precision this floor is usually far below any step one would choose, but with a low-order method in single precision it is reachable.

---

What to do instead. Refine, but verify rather than assume. Compare answers at $h$ and $h / 2$ : if they agree to the digits needed, those digits are probably sound; if they disagree wildly, or the sequence oscillates, suspect instability before suspecting the equation. And check the output against what the equation says. A solution of $y^{'} = - 20 y$ starting positive must stay positive.

Figure

Euler, improved Euler and RK4 at equal cost

same cost, three accuracies: order beats step count

The same budget of 32 evaluations of $f$ , spent three ways on $y^{'} = y$ , $y (0) = 1$ over $[0, 1]$ . Euler takes 32 steps, improved Euler 16, RK4 8.

The coarsest method by step count is the closest by result. Euler ends at $2.67699$ , improved Euler at $2.71659$ , RK4 at $2.7182768$ against the exact $e = 2.7182818$ — errors of $4.1 \times 10^{- 2}$ , $1.7 \times 10^{- 3}$ and $5.0 \times 10^{- 6}$ for identical cost.

Read the shapes, not just the endpoints. Euler's polyline falls below the curve everywhere and the deficit compounds, because each step uses a slope measured only at its left end on a function that is increasing. The higher-order methods sample within the step, so their error per step is a higher power of $h$ and the coarser grid is more than repaid.

Contrast

Same work, three methods

Comparing methods per step is misleading, because a step of RK4 costs four evaluations of $f$ and a step of Euler costs one. Evaluating $f$ is usually the expensive part, so the fair comparison holds the number of evaluations fixed.

The test. $y^{'} = y$ with $y (0) = 1$ on $[0, 1]$ , exact answer $e = 2.7182818285$ . Budget: 32 evaluations of $f$ . That buys 32 Euler steps, 16 improved Euler steps, or 8 RK4 steps.

Method	Steps	$h$	Result	Error
Euler	32	$0.03125$	$2.67699013$	$4.13 \times 10^{- 2}$
Improved Euler	16	$0.0625$	$2.7165935225$	$1.69 \times 10^{- 3}$
RK4	8	$0.125$	$2.718276844417$	$4.98 \times 10^{- 6}$

For identical cost, RK4 is about 8300 times more accurate than Euler and 340 times more accurate than improved Euler. The coarser step of the higher-order method is more than repaid by the exponent.

The gap widens with the budget. Doubling to 64 evaluations:

Method	Steps	Error at 32 evals	Error at 64 evals	Improvement
Euler	32 → 64	$4.13 \times 10^{- 2}$	$2.09 \times 10^{- 2}$	$\times 2$
Improved Euler	16 → 32	$1.69 \times 10^{- 3}$	$4.32 \times 10^{- 4}$	$\times 3.9$
RK4	8 → 16	$4.98 \times 10^{- 6}$	$3.28 \times 10^{- 7}$	$\times 15.2$

Each doubling of effort buys Euler one factor of 2 and RK4 a factor of 15. The advantage compounds, so the more accuracy is wanted, the more decisively the higher-order method wins.

A single RK4 step, costing four evaluations, gives error $9.9 \times 10^{- 3}$ , already better than Euler with 128 evaluations ( $1.05 \times 10^{- 2}$ ). That is a 32-fold cost saving on the crudest possible setting.

When the ranking does not hold.

When $f$ is not smooth. The order theorem assumes enough derivatives exist. Across a kink in $f$ , RK4's four samples straddle the discontinuity and its advantage collapses toward first order.
When stability binds. For a stiff equation the step is capped by the stability limit, not by accuracy. If $h$ must be tiny anyway, the extra accuracy per step is wasted, and an implicit method, which has no such cap, beats all three.
When only a rough answer is wanted. For one or two digits, Euler's simplicity can be worth more than RK4's accuracy, especially by hand.
When $f$ is cheap and coding time is not. Euler is three lines.

Higher order is not a free improvement, since it costs evaluations per step, but the exponent beats the constant factor as soon as any real accuracy is required. That is why RK4 is the default in step 2 of the procedure, and why the exceptions are about smoothness and stability rather than about accuracy.

Application

Numerical solution in applied settings

The exact methods in this module cover separable, linear, constant-coefficient and linear-system equations. That list sounds broad and is not. Most equations arising from a real situation fall outside it, and the reason is usually the same: something in the model is nonlinear.

Orbital motion. Newton's law of gravitation gives an inverse-square force, so the equations of motion carry $1 / r^{2}$ with $r$ itself depending on position. Two bodies can be solved exactly, that is Kepler's achievement. Three cannot, and no formula exists for the general case. Every trajectory flown to another planet is computed by stepping, and the step size is chosen against exactly the trade-off this unit describes.

Epidemics. The standard model has susceptible and infected populations with a transmission term proportional to the product $S I$ . That product is nonlinear, so nothing separates and no integrating factor applies, although the system is only two equations and every term has a plain meaning. Case-number projections come from numerical solution.

Weather and climate. The governing equations are nonlinear partial differential equations; discretising space turns them into an enormous system of ordinary differential equations in time, stepped forward numerically. Here the stability limit is not an academic caution. It dictates the time step given the grid spacing, and getting it wrong produces a forecast that blows up rather than one that is merely inaccurate.

Circuits and control. A resistor–inductor–capacitor circuit with linear components is solvable by the characteristic equation of the second-order unit. Add a diode or a transistor, whose current–voltage relation is exponential, and the equation becomes nonlinear. Circuit simulators step it numerically, which is what a SPICE simulation is.

Chemical kinetics. Reaction rates depend on products of concentrations, so kinetics is nonlinear almost by construction. These systems are also frequently stiff: a fast reaction and a slow one in the same mixture give widely separated decay rates, and the stability limit of an explicit method is set by the fastest, forcing tiny steps even when only the slow behaviour is of interest. This is the practical reason implicit methods were developed.

---

What changes when the answer is numerical.

	Exact solution	Numerical solution
Output	a function of $t$ and the parameters	numbers at chosen points, for one set of inputs
Changing a parameter	substitute and read off	run the computation again
Long-run behaviour	read from the formula's structure	inferred from a finite run, or from theory
Error	none	present, and requiring estimation

The second row is the real cost. A formula shows how the answer depends on its inputs, and a table of numbers does not, which is why a parameter study means many runs, and why the eigenvalue reading of the previous unit is so valuable: it extracts qualitative behaviour without solving at all.

Which is why step 1 of the procedure says to look for an exact solution first. Not from nostalgia, but because a formula answers questions a simulation cannot. When none exists, the normal case, the methods here supply numbers, and the order and stability theory says how far to trust them.

And why the earlier units still matter. Their equations serve as the test cases against which methods are checked: $y^{'} = y$ has the known answer $e$ , and the error tables in this unit were built by comparing against it. A method is trusted on equations that cannot be solved because it was verified on equations that can.

Numerical Solution of Differential Equations

When there is no formula to find

Why one slope per step is not enough

Reading the formulas: where each slope is evaluated

Approximating a solution numerically

Four Euler steps, then the same problem improved

An equation with no closed-form solution

Order of accuracy, and where the exponents come from

A smaller step is not always a better answer

Euler, improved Euler and RK4 at equal cost

Same work, three methods

Numerical solution in applied settings

Settings

Focus timer

Today

Last seven days

Where the time went

Recent sessions

Settings

Notifications

Your history