Numerical Solution of Differential Equations

Stepping methods that approximate a solution when no formula exists, carried out by hand for Euler and improved Euler, and separately judged: the order of accuracy that says how each responds to a smaller step, the stability limit that can make a smaller step necessary rather than merely better, and what a computed answer does and does not establish.

Definition

A numerical method for the initial value problem

y^{'} = f (t, y), y (t_{0}) = y_{0}

produces approximations $y_{1}, y_{2}, \dots$ to $y (t_{1}), y (t_{2}), \dots$ at points $t_{n} = t_{0} + n h$ , where $h$ is the step size. It computes values, not a formula.

Euler's method. The differential equation gives the slope at any point, so follow it for one step:

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

This is the tangent line to the solution through $(t_{n}, y_{n})$ , evaluated at $t_{n + 1}$ .

Improved Euler (Heun's method). The slope at the start of a step is generally not the slope across it, so average the slope at both ends, using an Euler step to predict the far end:

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}) .

Classical Runge–Kutta (RK4). Four slope evaluations, weighted toward the midpoint:

\begin{aligned} k_{1} & = f (t_{n}, y_{n}), & k_{2} & = f (t_{n} + \frac{h}{2}, y_{n} + \frac{h}{2} k_{1}), \\ k_{3} & = f (t_{n} + \frac{h}{2}, y_{n} + \frac{h}{2} k_{2}), & k_{4} & = f (t_{n} + h, y_{n} + h k_{3}), \end{aligned}

y_{n + 1} = y_{n} + \frac{h}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4}) .

Order of accuracy. A method has order $p$ when the error at a fixed final time behaves as $O (h^{p})$ . Halving $h$ then divides the error by $2^{p}$ :

Method	Slope evaluations per step	Order $p$	Error when $h$ is halved
Euler	1	1	halved
Improved Euler	2	2	quartered
RK4	4	4	divided by 16

Order describes the trend as $h \to 0$ , not the error at any particular $h$ .

Stability. Accuracy is not the only constraint. Applying Euler to $y^{'} = λ y$ gives $y_{n + 1} = (1 + h λ) y_{n}$ , so the computed values stay bounded only when

| 1 + h λ | < 1 .

For $λ < 0$ this forces $h < 2 / | λ |$ . A step above that limit produces growing, sign-alternating output no matter how smooth the true solution is. A failure of the method, not of the problem.

Assumptions and scope

These methods approximate. They return numbers at chosen points, never a formula, and no amount of refinement produces one.
Order describes asymptotic behaviour as $h \to 0$ . Observed error ratios approach $2$ , $4$ and $16$ rather than equalling them, and at coarse steps they can be well short.
Stability is a separate requirement from accuracy. Euler on $y^{'} = λ y$ with $λ < 0$ needs $h < 2 / | λ |$ regardless of the accuracy sought.
Halving the step doubles the number of steps, so error per step and total work move in opposite directions; the comparison between methods is per unit of work, not per step.
Round-off sets a floor. Below a certain step size, accumulated floating-point error grows while truncation error shrinks, so the total stops improving and eventually worsens.
The methods assume $f$ is defined along the whole path taken. A solution that escapes to infinity in finite time, as $y^{'} = y^{2}$ does, will produce numbers past the blow-up that mean nothing.

Worked material

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.

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.

Common errors

Common misconception

A smaller step size always gives a better numerical answer, so accuracy is limited only by how long one is willing to compute.

Common misconception

The step size in a numerical method is chosen from the accuracy wanted, so a problem needing only one or two digits may always be integrated with a coarse step.

Numerical Solution of Differential Equations

Assumptions and scope

Worked material

An equation with no closed-form solution

Same work, three methods

Common errors

Related units

Requires

Connected

Learn this topic

Used in

Sources

Settings

Focus timer

Today

Last seven days

Where the time went

Recent sessions

Settings

Notifications

Your history