Forward Euler Method

The Initial Value Problem¶

We seek to approximate the solution of

on a grid of times $t_0 < t_1 < t_2 < \cdots < t_N = T$ with step size $h = t_{n+1} - t_n$.

Let $u_n$ denote our approximation to $u(t_n)$.

Derivation¶

Replace $du/dt$ with a forward difference:

This gives the forward Euler (or explicit Euler) method:

Geometric Interpretation¶

Forward Euler follows the tangent line:

1. At $(t_n, u_n)$, the slope is $f(t_n, u_n)$

2. Follow this slope for time $h$

3. Arrive at $u_{n+1} = u_n + h \cdot \text{slope}$

Each step follows a tangent to a (possibly different) solution curve.

Implementation¶

See the Euler’s Method notebook for an interactive implementation with visualizations of convergence, stability, and stiffness.

Sensitivity of the Initial Value Problem¶

Before analyzing the algorithm’s accuracy, we should ask: how sensitive is the problem itself?

The colored curves in the figure above are exact solutions of the same ODE with different initial conditions. How quickly do nearby solutions diverge or converge?

Lipschitz Continuity¶

Condition Number of the IVP¶

Suppose $u(t)$ and $v(t)$ both solve $u' = f(t,u)$ but with different initial conditions $u(0) = u_0$ and $v(0) = v_0$. Their difference satisfies

If $f$ is Lipschitz continuous in $u$ with constant $L$ (Definition 2), then $\frac{d}{dt}|u - v| \leq L\,|u - v|$. By Grönwall’s inequality:

The factor $e^{LT}$ is the condition number of the IVP. It is a property of the problem, not the algorithm:

• If $L > 0$: nearby solution curves diverge. Small perturbations grow exponentially, and $e^{LT}$ can be very large. The problem is ill-conditioned over long times.

• If $L < 0$: nearby solution curves converge. Perturbations are damped, and $e^{LT} < 1$. The problem is well-conditioned.

• If $L = 0$ (e.g., $u' = g(t)$): solution curves are parallel. Errors neither grow nor shrink.

No numerical method can do better than the problem allows. We will see this same factor reappear in the convergence bound (Theorem 1).

Local Truncation Error¶

If we take a single step of forward Euler starting from the exact solution, how far off are we? The local truncation error measures this per-step accuracy: the error committed in one step, assuming we started with the exact answer.

Consistency Order¶

This is the same notion of “order of accuracy” we saw for finite differences: the forward difference $\frac{f(x+h) - f(x)}{h}$ approximates $f'(x)$ with error $O(h)$. Forward Euler is a forward difference applied to $du/dt$, so it inherits that same first-order accuracy.

Global Error and Convergence¶

The local truncation error tells us how accurate a single step is. But we take many steps, so the real question is: do these per-step errors accumulate controllably, or do they blow up?

The global error is the cumulative effect of all one-step errors. Two things determine its size:

1. Error accumulation: Recall that the one-step error is $\mathcal{L}_n = h\tau_n = O(h^{p+1})$. Over $N = T/h$ steps, these pile up: $N \times O(h^{p+1}) = O(h^p)$.

2. Stability: Local errors don’t just add up; each one gets amplified by subsequent steps. If the amplification factor is bounded, errors accumulate gently ($O(h^p)$). If it exceeds 1, errors grow exponentially and the method is useless.

Convergence Theorem¶

Stability Analysis¶

The convergence theorem guarantees that errors vanish as $h \to 0$. But in practice we use a fixed step size. Can errors grow from step to step?

The error propagation example (Example 2) already answered this for the linear problem $u' = \lambda u + g(t)$: the error at each step is multiplied by the amplification factor $R = 1 + h\lambda$. If $|R| > 1$, errors grow exponentially and the method is useless.

This is why the scalar test equation $u' = \lambda u$ plays a central role. Any smooth ODE is locally linear (by Taylor expansion), so the behavior of a method on $u' = \lambda u$ reveals its stability properties in general. Applying forward Euler gives

so the numerical solution is $u_n = (1 + h\lambda)^n u_0$. For the solution to remain bounded, we need $|1 + h\lambda| \leq 1$. Writing $z = h\lambda$, this defines the stability region.

Definition 6 (Stability Region)

The stability region of forward Euler is:

$$\mathcal{S} = \{z \in \mathbb{C} : |1 + z| \leq 1\}$$

(31)

This is a disk of radius 1 centered at $(-1, 0)$ in the complex plane.

For real $\lambda < 0$, stability requires $-2 \leq h\lambda \leq 0$, giving the step-size restriction

Forward Euler is conditionally stable: the step size is restricted by the problem’s eigenvalues.

Stiffness¶

When $|\lambda|$ is large, the stability constraint $h \leq 2/|\lambda|$ forces extremely small step sizes even when the solution itself varies slowly. Stiffness is not about the solution being complicated; it is about the equation pulling toward a slow manifold so aggressively that explicit methods cannot keep up.

Example 3 (Stiffness in Action (Prothero-Robinson))

Consider $u' = \lambda(u - \cos t) - \sin t$ with $u(0) = 1$. The exact solution is $u(t) = \cos t$ regardless of $\lambda$: the same smooth, slowly varying function whether $\lambda = -1$ or $\lambda = -10^6$.

Since $u(t) = \cos t$, we have $|u''(t)| = |\cos t| \leq 1$.

Accuracy step size. By Proposition 1, the one-step error is $\mathcal{L}_n = h\tau_n = \frac{h^2}{2}u''(t_n) + O(h^3)$. For the Prothero-Robinson problem $|u''(t)| = |\cos t| \leq 1$, so

$$|\mathcal{L}_n| \leq \frac{h^2}{2}$$

(34)

To keep the per-step error below a tolerance $\delta$, we need $h \leq \sqrt{2\delta}$. This depends only on the smoothness of the solution, not on $\lambda$.

Over $N = T/h$ steps, the one-step errors accumulate to a global error of roughly $|E_N| \approx N \cdot |\mathcal{L}| = \frac{T}{h} \cdot \frac{h^2}{2} = \frac{Th}{2}$. For a global tolerance $\varepsilon = 10^{-2}$ over $T = 2$:

$$h_{\text{accuracy}} \approx \frac{2\varepsilon}{T} = \frac{2 \times 10^{-2}}{2} = 10^{-2}$$

(35)

This does not involve $\lambda$ at all.

Stability step size. From the stability analysis above, forward Euler requires $h < 2/|\lambda|$:

$\lambda$	$h_{\text{stability}}$	$h_{\text{accuracy}}$	Ratio	Steps needed
-10	0.2	0.01	0.05 (not stiff)	200
-1000	0.002	0.01	5 (stiff)	1000
-106	$2 \times 10^{-6}$	0.01	5000	106

The solution is identical in every row; only the stability constraint changes. With $\lambda = -10^6$, forward Euler needs a million steps to integrate to $T = 2$, all wasted on stability rather than accuracy. An implicit method (backward Euler) can take $h = 0.01$ and finish in 200 steps.

Why the convergence theorem cannot see stiffness¶

The convergence theorem (Theorem 1) uses the Lipschitz constant $L = |\lambda|$ as its measure of problem sensitivity. For stiff problems this gives a catastrophically pessimistic bound.

What About Rounding Errors?¶

The convergence theorem (Theorem 1) assumes exact arithmetic. On a computer, each step incurs a rounding error $\mu_n$ in addition to the truncation error.

The bound has two competing terms as $h \to 0$:

• Truncation error $Ch^p$: decreases with smaller $h$.

• Rounding error $\max \|\mu_l\|/h$: increases with smaller $h$, since we divide the per-step rounding error (roughly machine epsilon) by an ever-smaller step size.

This creates three regimes:

1. Large $h$ — the discretization is too coarse to resolve the solution; errors decrease as $h$ shrinks.

2. Moderate $h$ — the “sweet spot” where truncation error dominates and convergence at rate $O(h^p)$ is observed.

3. Very small $h$ — rounding errors dominate and the total error increases as $h$ shrinks further.

There is an optimal step size $h_{\text{opt}}$ that balances the two contributions. Below $h_{\text{opt}}$, making $h$ smaller makes things worse.

Limitations¶

Low order. Forward Euler is only $O(h)$ accurate. Just as central differences improved on forward differences by averaging, we can build higher-order ODE methods by evaluating $f$ at cleverly chosen intermediate points and combining the results. This is the idea behind Runge-Kutta methods; see the optional section on Runge-Kutta and adaptive methods.

Fixed step size. For a nonlinear ODE, the effective $\lambda$ changes along the solution. A step size that is stable in one region may be unstable in another, and a step size that is accurate during rapid transients may be wasteful during slow evolution. With a fixed $h$, there is no way to adapt. This motivates adaptive time-stepping, where we estimate the local truncation error at each step and adjust $h$ automatically; see the optional section on Runge-Kutta and adaptive methods.