Forward Euler Method - Introduction to Scientific Computing

The Initial Value Problem¶

We seek to approximate the solution of

\frac{du}{dt} = f(t, u), \quad u(t_0) = u_0

(1)

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:

\frac{u(t_{n+1}) - u(t_n)}{h} \approx f(t_n, u(t_n))

(2)

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

Geometric Interpretation¶

Forward Euler follows the tangent line:

At $(t_n, u_n)$ , the slope is $f(t_n, u_n)$
Follow this slope for time $h$
Arrive at $u_{n+1} = u_n + h \cdot \text{slope}$

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

Implementation¶

def forward_euler(f, u0, t0, T, h):
    """
    Solve u' = f(t, u) from t0 to T with step size h.
    """
    t = t0
    u = u0
    times = [t]
    values = [u]

    while t < T:
        u = u + h * f(t, u)
        t = t + h
        times.append(t)
        values.append(u)

    return np.array(times), np.array(values)

Local Truncation Error¶

The local truncation error measures how well the exact solution satisfies the numerical scheme.

For forward Euler, Taylor expand the exact solution:

u(t_{n+1}) = u(t_n) + h u'(t_n) + \frac{h^2}{2}u''(t_n) + O(h^3)

(5)

Since $u'(t_n) = f(t_n, u(t_n))$ :

\tau_n = \frac{h}{2}u''(t_n) + O(h^2) = O(h)

(6)

Global Error and Convergence¶

The global error $E_n = u_n - u(t_n)$ accumulates local errors over all steps.

Error Propagation: Building Intuition¶

To understand how local errors accumulate, consider the linear test problem $u' = \lambda u + g(t)$ . The exact solution satisfies:

u(t_{n+1}) = (1 + h\lambda)u(t_n) + hg(t_n) + h\tau_n

(9)

The numerical solution satisfies:

u_{n+1} = (1 + h\lambda)u_n + hg(t_n)

(10)

Subtracting gives the error recurrence:

E_{n+1} = (1 + h\lambda)E_n - h\tau_n

(11)

Unrolling this recurrence:

E_n = (1 + h\lambda)^n E_0 - h\sum_{j=0}^{n-1}(1 + h\lambda)^{n-1-j}\tau_j

(12)

This is the discrete Duhamel principle: the global error at step $n$ equals the initial error propagated forward, plus all local truncation errors propagated by the appropriate powers of the amplification factor.

Remark 2 (Connection to Neumann Series)

This structure is identical to the Neumann series from iterative methods! Compare:

Context	Recurrence	Solution
Iterative methods	$x_{k+1} = Ax_k + b$	$x_n = A^n x_0 + \sum_{j=0}^{n-1} A^{n-1-j}b$
ODE error	$E_{n+1} = RE_n - h\tau_n$	$E_n = R^n E_0 - h\sum_{j=0}^{n-1} R^{n-1-j}\tau_j$

where $R = 1 + h\lambda$ is the amplification factor.

Both are discrete convolutions with a geometric kernel. The convergence/stability condition is the same: $|R| < 1$ (ODE stability) parallels $\|A\| < 1$ (Neumann series convergence).

This also discretizes the continuous variation of constants formula:

u(t) = e^{\lambda t}u_0 + \int_0^t e^{\lambda(t-s)}g(s)\,ds

(13)

The exponential $e^{\lambda t}$ becomes the power $R^n = (1+h\lambda)^n$ , and the integral becomes a sum.

Convergence Theorem¶

Proof 1

We prove this for forward Euler on a nonlinear problem; the general case is similar.

Step 1: Error recurrence. The numerical method gives $u_{n+1} = u_n + hf(t_n, u_n)$ . The exact solution satisfies $u(t_{n+1}) = u(t_n) + hf(t_n, u(t_n)) + h\tau_n$ .

Subtracting:

E_{n+1} = E_n + h[f(t_n, u_n) - f(t_n, u(t_n))] - h\tau_n

(16)

Step 2: Lipschitz bound. Since $f$ is Lipschitz with constant $L$ :

|f(t_n, u_n) - f(t_n, u(t_n))| \leq L|E_n|

(17)

Thus:

|E_{n+1}| \leq |E_n| + hL|E_n| + h|\tau_n| = (1 + hL)|E_n| + h|\tau_n|

(18)

Step 3: Unroll the recurrence. Let $\tau_{\max} = \max_j |\tau_j|$ . Then:

|E_n| \leq (1 + hL)^n|E_0| + h\tau_{\max}\sum_{j=0}^{n-1}(1 + hL)^j

(19)

Step 4: Bound the geometric sum. Using $(1 + hL)^k \leq e^{khL}$ and $\sum_{j=0}^{n-1}(1+hL)^j \leq \frac{(1+hL)^n - 1}{hL}$ :

|E_n| \leq e^{nhL}|E_0| + \frac{e^{nhL} - 1}{L}\tau_{\max}

(20)

Step 5: Substitute $nh = t_n - t_0$ and $\tau_{\max} = O(h^p)$ . With $E_0 = 0$ (exact initial condition):

|E_n| \leq \frac{e^{L(t_n - t_0)} - 1}{L} \cdot Ch^p = O(h^p)

(21)

Proof 2

This follows directly from Theorem 1. If the method has consistency order $p$ , then $\|\tau\| \leq Ch^p$ . The theorem gives:

\|E_n\| \leq \frac{C}{L_\Psi}\left(e^{L_\Psi T} - 1\right) h^p = O(h^p)

(23)

The key point: the constant $\frac{C}{L_\Psi}(e^{L_\Psi T} - 1)$ depends on the problem (through $L_\Psi$ and $T$ ) but not on $h$ . Thus consistency of order $p$ implies convergence of order $p$ .

Remark 3 (The Forward/Backward Error Framework for ODEs)

The convergence theorem is the ODE version of our fundamental error relationship:

\text{Forward error} \leq \text{Condition number} \times \text{Backward error}

(24)

In the ODE setting:

Backward error = Local truncation error $\tau_n = O(h^p)$
Forward error = Global error $\|E_n\| = O(h^p)$
Condition number = $e^{L(t_n - t_0)}$ (Lipschitz amplification over time)

The Lipschitz constant $L$ measures how sensitive the ODE is to perturbations. The factor $e^{LT}$ is the condition number for time evolution—a property of the problem, not the algorithm.

A consistent method (small backward error per step) achieves convergence (controlled forward error) as long as the problem isn’t too ill-conditioned (moderate $L$ ).

Stability Analysis¶

Convergence (as $h \to 0$ ) doesn’t guarantee useful results for fixed $h > 0$ . We need absolute stability.

The convergence theorem tells us that errors are bounded as $h \to 0$ . But in practice, we use a fixed step size $h > 0$ . The question becomes: does the error recurrence

E_{n+1} = (1 + h\lambda)E_n - h\tau_n

(25)

produce bounded errors for our chosen $h$ ? This depends on the amplification factor $R = 1 + h\lambda$ . If $|R| > 1$ , errors grow exponentially with each step—the method is unstable. Stability is the condition that ensures our backward errors (local truncation errors) don’t get amplified catastrophically as we march forward in time.

The Test Equation¶

Apply forward Euler to $u' = \lambda u$ with $\text{Re}(\lambda) < 0$ (decaying solutions):

u_{n+1} = (1 + h\lambda)u_n = R(z)u_n, \quad z = h\lambda

(26)

where $R(z) = 1 + z$ is the stability function.

Stability Region¶

For the numerical solution to remain bounded, we need $|R(z)| = |1 + z| \leq 1$ .

For real $\lambda < 0$ : stability requires $-2 \leq h\lambda \leq 0$ , so:

h \leq \frac{2}{|\lambda|}

(28)

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

Remark 4 (The Exponential Factor)

The factor $e^{L(t_n - t_0)}$ in the convergence bound is unavoidable—it reflects the problem’s sensitivity (condition number), not the algorithm. The Lipschitz constant $L$ measures how quickly nearby solution curves diverge or converge:

If $L = 0$ (e.g., $u' = g(t)$ ): solution curves are parallel, errors don’t grow
If $L > 0$ : solution curves diverge, errors can grow exponentially
If $L < 0$ : solution curves converge, errors are damped

Limitations¶

Forward Euler’s conditional stability becomes problematic for stiff problems—where large eigenvalues force tiny step sizes even when the solution varies slowly.

For the test equation with $\lambda = -1000$ :

Stability requires $h < 0.002$
But the solution $e^{-1000t}$ decays to negligible levels by $t = 0.01$

We’re forced to take thousands of tiny steps tracking a transient we don’t care about. This motivates implicit methods like backward Euler, which we develop next.