Numerical approximation when exact solutions are impossible
After completing this section, you should be able to:
Many differential equations lack closed-form analytical solutions. When we cannot solve $y' = f(x, y)$ explicitly, we turn to numerical methods that approximate the solution step by step, following the direction field at each point. Euler's method is the simplest and most intuitive approach to this problem.
Given the initial value problem
we want to find approximate values of $y$ at points $x_1, x_2, \ldots$ ahead of $x_0$.
At each current point $(x_n, y_n)$, the slope of the solution curve is $f(x_n, y_n)$. Instead of following the actual curve, Euler's method simply follows the tangent line for a small step of size $h$:
where $x_{n+1} = x_n + h$ and $h$ is the step size (also called mesh width).
Starting from $(x_0, y_0)$, compute successive approximations via: $$x_{n+1} = x_n + h$$ $$y_{n+1} = y_n + h \cdot f(x_n, y_n)$$ for $n = 0, 1, 2, \ldots$ until reaching the target $x$ value.
At each point $(x_n, y_n)$, we draw the tangent line to the solution curve (which has slope $f(x_n, y_n)$). We then move horizontally a distance $h$, and vertically an amount $h \cdot f(x_n, y_n)$ to find the next approximation point. This is like walking along the direction field with a fixed step size.
Local Truncation Error (LTE): The error introduced in a single step, assuming all previous steps are exact. For Euler's method, LTE $= O(h^2)$.
Global Truncation Error (GTE): The total accumulated error over all steps from $x_0$ to $x_f$. For Euler's method, GTE $= O(h)$.
This makes Euler a first-order method. Doubling the number of steps (halving $h$) roughly halves the global error, but at the cost of doing twice the work.
Halving the step size $h$ improves accuracy significantly, but it doubles the computational cost. In practice, we often need to balance accuracy against computational feasibility.
A better strategy is to use a two-stage corrector method. First, we predict using Euler's method, then we correct using the average slope between the current point and the predicted next point.
Predictor step: $$\tilde{y}_{n+1} = y_n + h \cdot f(x_n, y_n)$$
Corrector step: $$y_{n+1} = y_n + \frac{h}{2} \left[ f(x_n, y_n) + f(x_{n+1}, \tilde{y}_{n+1}) \right]$$
The improved Euler method (also called Heun's method) has:
For certain differential equations (especially stiff equations), if the step size $h$ is too large, the numerical solution can become unstable or diverge. Checking that your computed values remain reasonable is an important sanity check. If in doubt, use a smaller step size.
To test stability, compute the solution with step sizes $h$ and $h/2$. If the two results agree to your desired precision, you can trust the smaller-$h$ result. This is the basis of adaptive step-size control in modern codes.
This is the exponential growth equation. The exact solution is $y(x) = e^x$. We'll take 5 steps with $h = 0.1$ to reach $x = 0.5$.
Setup: $f(x, y) = y$, $(x_0, y_0) = (0, 1)$, $h = 0.1$
Computations:
| $n$ | $x_n$ | $y_n$ (Euler) | $y_{\text{exact}}$ | Error |
|---|---|---|---|---|
| 0 | 0.0 | 1.0000 | 1.0000 | 0.0000 |
| 1 | 0.1 | 1.1000 | 1.1052 | 0.0052 |
| 2 | 0.2 | 1.2100 | 1.2214 | 0.0114 |
| 3 | 0.3 | 1.3310 | 1.3499 | 0.0189 |
| 4 | 0.4 | 1.4641 | 1.4918 | 0.0277 |
| 5 | 0.5 | 1.6105 | 1.6487 | 0.0382 |
Step-by-step:
Step 0 to 1: $m_0 = f(0, 1) = 1 \Rightarrow y_1 = 1 + 0.1(1) = 1.1$
Step 1 to 2: $m_1 = f(0.1, 1.1) = 1.1 \Rightarrow y_2 = 1.1 + 0.1(1.1) = 1.21$
Continue similarly for steps 3, 4, and 5.
Result: $y(0.5) \approx 1.6105$ (Euler) vs. $e^{0.5} \approx 1.6487$ (exact). The error grows because we're accumulating local errors.
This is a linear first-order ODE. The exact solution is $y(x) = 2e^x - x - 1$. We'll take 5 steps with $h = 0.2$ from $x = 0$ to $x = 1$.
Setup: $f(x, y) = x + y$, $(x_0, y_0) = (0, 1)$, $h = 0.2$
| $n$ | $x_n$ | $f(x_n, y_n)$ | $y_n$ (Euler) | $y_{\text{exact}}$ | Error |
|---|---|---|---|---|---|
| 0 | 0.0 | 1.0 | 1.0000 | 1.0000 | 0.0000 |
| 1 | 0.2 | 1.4000 | 1.2000 | 1.2428 | 0.0428 |
| 2 | 0.4 | 1.8800 | 1.4800 | 1.5836 | 0.1036 |
| 3 | 0.6 | 2.4560 | 1.8560 | 2.0442 | 0.1882 |
| 4 | 0.8 | 3.1472 | 2.3472 | 2.6511 | 0.3039 |
| 5 | 1.0 | — | 2.9766 | 3.4366 | 0.4600 |
Result: $y(1) \approx 2.9766$ (Euler) vs. $2e - 2 \approx 3.4366$ (exact). The error is substantial with this step size, showing why finer grids or better methods are needed.
Same problem as Example 1, but now using the improved Euler (Heun's method) to show the dramatic improvement in accuracy.
Setup: $f(x, y) = y$, $(x_0, y_0) = (0, 1)$, $h = 0.1$
| $n$ | $x_n$ | $\tilde{y}_{n+1}$ (Pred.) | $y_n$ (Improved) | $y_{\text{exact}}$ | Error |
|---|---|---|---|---|---|
| 0 | 0.0 | — | 1.0000 | 1.0000 | 0.0000 |
| 1 | 0.1 | 1.1000 | 1.1050 | 1.1052 | 0.0002 |
| 2 | 0.2 | 1.2155 | 1.2214 | 1.2214 | 0.0000 |
| 3 | 0.3 | 1.3436 | 1.3499 | 1.3499 | 0.0000 |
| 4 | 0.4 | 1.4849 | 1.4918 | 1.4918 | 0.0000 |
| 5 | 0.5 | 1.6410 | 1.6487 | 1.6487 | 0.0000 |
Detail for Step 1:
Predictor: $\tilde{y}_1 = 1 + 0.1 \cdot 1 = 1.1$
Corrector: $y_1 = 1 + \frac{0.1}{2}[1 + 1.1] = 1 + 0.05(2.1) = 1.105$
Observation: The improved Euler method is dramatically more accurate! The errors are on the order of $10^{-4}$ or better, compared to $10^{-2}$ for basic Euler. This is the power of second-order methods.
This equation has the exact solution $y(x) = e^{-x^2}$ (a Gaussian curve). We'll use Euler's method to see how it tracks this smooth, decaying function.
Setup: $f(x, y) = -2xy$, $(x_0, y_0) = (0, 1)$, $h = 0.1$
| $n$ | $x_n$ | $m_n = -2x_n y_n$ | $y_n$ (Euler) | $y_{\text{exact}}$ | Error |
|---|---|---|---|---|---|
| 0 | 0.0 | 0.0 | 1.0000 | 1.0000 | 0.0000 |
| 1 | 0.1 | -0.2000 | 1.0000 | 0.9900 | 0.0100 |
| 2 | 0.2 | -0.3920 | 0.9800 | 0.9608 | 0.0192 |
| 3 | 0.3 | -0.5645 | 0.9408 | 0.9139 | 0.0269 |
| 4 | 0.4 | -0.7075 | 0.8844 | 0.8521 | 0.0322 |
| 5 | 0.5 | — | 0.8136 | 0.7788 | 0.0348 |
Result: The method captures the general decay, but as $x$ increases, the slope becomes steeper (in absolute value), and Euler's method underestimates the decay. Using a smaller step size would help.
To illustrate the convergence of Euler's method as we decrease the step size, we solve the same problem with different values of $h$. The exact answer is $e \approx 2.71828$.
| Step Size $h$ | Number of Steps | $y_{\text{approx}}(1)$ | Error | Error Ratio |
|---|---|---|---|---|
| 0.5 | 2 | 2.2500 | 0.4683 | — |
| 0.25 | 4 | 2.4414 | 0.2769 | 1.69 |
| 0.1 | 10 | 2.5937 | 0.1246 | 2.22 |
| 0.05 | 20 | 2.6533 | 0.0650 | 1.92 |
Observation: As $h$ decreases, the error roughly decreases proportionally to $h$ (the error ratio is close to $h_{\text{old}} / h_{\text{new}}$). This confirms that Euler's method is first-order: error $\propto h$.
For example, from $h = 0.5$ to $h = 0.25$ (halving $h$), the error ratio is approximately $2$. From $h = 0.25$ to $h = 0.1$, the ratio is about $2.2$. This linear convergence is characteristic of first-order methods.
These are written in the style of past exam questions on numerical solutions of first-order IVPs. Carry at least six decimal places while iterating — rounding early is the single most common source of lost marks. Click a problem to reveal the full solution.
Step 1: Identify the pieces
$$f(x,y) = x + y,\qquad (x_0,y_0) = (0,1),\qquad h = 0.1,\qquad y_{n+1} = y_n + h\,f(x_n,y_n)$$Step 2: Find the exact solution for comparison
This is linear: $y' - y = x$, with integrating factor $e^{-x}$, giving $\left(e^{-x}y\right)' = xe^{-x}$ and hence
$$y(x) = 2e^{x} - x - 1$$Check: $y(0) = 2 - 0 - 1 = 1$ ✓ and $y' = 2e^x - 1 = x + (2e^x - x - 1) = x + y$ ✓
Step 3: Iterate five steps
$$y_1 = 1 + 0.1(0 + 1) = 1.1,\qquad y_2 = 1.1 + 0.1(0.1 + 1.1) = 1.22$$ $$y_3 = 1.22 + 0.1(0.2 + 1.22) = 1.362,\qquad y_4 = 1.362 + 0.1(0.3 + 1.362) = 1.5282$$ $$y_5 = 1.5282 + 0.1(0.4 + 1.5282) = 1.72102$$| $n$ | $x_n$ | $y_n$ (Euler) | $y(x_n)$ exact | Error |
|---|---|---|---|---|
| 0 | 0.0 | 1.000000 | 1.000000 | 0.000000 |
| 1 | 0.1 | 1.100000 | 1.110342 | 0.010342 |
| 2 | 0.2 | 1.220000 | 1.242806 | 0.022806 |
| 3 | 0.3 | 1.362000 | 1.399718 | 0.037718 |
| 4 | 0.4 | 1.528200 | 1.583649 | 0.055449 |
| 5 | 0.5 | 1.721020 | 1.797443 | 0.076423 |
Step 4: Confirm the first-order behaviour
Repeating with $h = 0.05$ (ten steps) gives $y(0.5) \approx 1.757789$, an error of $0.039653$. The error ratio is
$$\frac{0.076423}{0.039653} \approx 1.93 \approx 2$$exactly what an $O(h)$ global error predicts.
$$\boxed{y(0.5) \approx 1.721020,\qquad \text{exact } 1.797443,\qquad \text{error } \approx 0.0764}$$Step 1: Write down the predictor–corrector pair
$$k_1 = f(x_n, y_n),\qquad y^{*}_{n+1} = y_n + h\,k_1,\qquad k_2 = f(x_n + h,\; y^{*}_{n+1})$$ $$y_{n+1} = y_n + \frac{h}{2}\left(k_1 + k_2\right)$$Here $f(x,y) = 2xy$, $h = 0.25$, and the exact solution is $y = e^{x^2}$ (separable: $dy/y = 2x\,dx$).
Step 2: First step, $x_0 = 0 \to x_1 = 0.25$
$$k_1 = 2(0)(1) = 0,\qquad y^{*}_1 = 1 + 0.25(0) = 1,\qquad k_2 = 2(0.25)(1) = 0.5$$ $$y_1 = 1 + \frac{0.25}{2}(0 + 0.5) = 1 + 0.0625 = 1.0625$$Exact: $e^{0.0625} = 1.064494$, so the error after one step is only $0.001994$.
Step 3: Second step, $x_1 = 0.25 \to x_2 = 0.5$
$$k_1 = 2(0.25)(1.0625) = 0.53125$$ $$y^{*}_2 = 1.0625 + 0.25(0.53125) = 1.1953125$$ $$k_2 = 2(0.5)(1.1953125) = 1.1953125$$ $$y_2 = 1.0625 + \frac{0.25}{2}(0.53125 + 1.1953125) = 1.0625 + 0.2158203125 = 1.2783203125$$Step 4: Compare
| Method ($h = 0.25$) | $y(0.5)$ | Error |
|---|---|---|
| Euler | 1.125000 | 0.159025 |
| Improved Euler | 1.278320 | 0.005705 |
| Exact $e^{0.25}$ | 1.284025 | — |
The corrector step costs one extra evaluation of $f$ per step and buys a factor of about $28$ in accuracy here.
$$\boxed{y(0.5) \approx 1.278320 \quad (\text{exact } 1.284025,\ \text{error } \approx 5.7\times10^{-3})}$$Step 1: Reduce one Euler step to a multiplication
With $f(x,y) = \lambda y$ and $\lambda = -20$,
$$y_{n+1} = y_n + h(\lambda y_n) = (1 + h\lambda)y_n \quad\Longrightarrow\quad y_n = (1 + h\lambda)^n y_0$$The amplification factor is $G = 1 + h\lambda = 1 - 20h$.
Step 2: Impose $|G| < 1$
$$|1 - 20h| < 1 \quad\Longleftrightarrow\quad -1 < 1 - 20h < 1 \quad\Longleftrightarrow\quad 0 < h < 0.1$$The exact solution $y = e^{-20x}$ decays, so any $h \ge 0.1$ produces a numerical solution that does not decay — no matter how small the local truncation error is.
Step 3: $h = 0.2$ — unstable, $G = -3$
| $n$ | $x_n$ | $y_n$ | Exact $e^{-20x_n}$ |
|---|---|---|---|
| 1 | 0.2 | $-3$ | 0.018316 |
| 2 | 0.4 | $9$ | 0.000335 |
| 3 | 0.6 | $-27$ | 0.000006 |
| 4 | 0.8 | $81$ | 0.000000 |
The signs alternate and the magnitude triples every step — the classic signature of an unstable explicit step on a stiff problem.
Step 4: $h = 0.02$ — stable, $G = 0.6$
| $n$ | $x_n$ | $y_n$ | Exact $e^{-20x_n}$ |
|---|---|---|---|
| 1 | 0.02 | 0.600000 | 0.670320 |
| 2 | 0.04 | 0.360000 | 0.449329 |
| 3 | 0.06 | 0.216000 | 0.301194 |
| 4 | 0.08 | 0.129600 | 0.201897 |
Now the numerical solution decays monotonically like the true one. It is still not especially accurate — that is a separate question from stability.
$$\boxed{\text{Euler is stable for }0 < h < \tfrac{1}{10};\quad h = 0.2 \Rightarrow y_n = (-3)^n \to \pm\infty}$$Test your understanding with these multiple-choice questions. Click a choice to check your answer.
$$y_{n+1} = y_n + h \cdot f(x_n, y_n)$$ $$x_{n+1} = x_n + h$$
Order: 1st (error $\propto h$)
$$\tilde{y}_{n+1} = y_n + h f(x_n, y_n)$$ $$y_{n+1} = y_n + \frac{h}{2}[f(x_n, y_n) + f(x_{n+1}, \tilde{y}_{n+1})]$$
Order: 2nd (error $\propto h^2$)
Euler: LTE $O(h^2)$, GTE $O(h)$
Improved: LTE $O(h^3)$, GTE $O(h^2)$
Better accuracy requires smaller $h$ or higher-order methods.
• Smaller $h$ = more accurate but slower
• Test with $h$ and $h/2$ to verify stability
• Improved Euler gives $\approx 4\times$ better accuracy for same work