Notebooks · posts/forward-euler-blows-up-past-h-equals-2-over-lambda

Forward and backward Euler on either side of h = 2/λ

Adam Zettel · September 25, 2026 · written with an agent

The exact solution of this equation is $cos t$ , which a step of 0.1 would follow to a few digits. With $λ = 1000$ , forward Euler is within $4 \times 1 0^{- 7}$ of it at $t = 2$ with $h = 0.0019$ , and off by $2.7 \times 1 0^{33}$ at $h = 0.0021$ . A disturbance is multiplied by $1 - hλ$ each step, and the measured factor matches $∣1 - hλ ∣$ to $2 \times 1 0^{- 14}$ at 36 step sizes. Backward Euler's factor is $1/ (1 + hλ)$ , below 1 at every step size, and at $h = 0.1$ it is still within $2 \times 1 0^{- 5}$ of $cos 2$ .

y^{'} = - λ (y - cos t) - sin t, ∣1 - hλ ∣ < 1 ⟺ h < 2/ λ

An equation is stiff when it has a component that decays much faster than the solution you want changes. Here the solution is $cos t$ , and any other starting value is pulled onto it at the rate $λ = 1000$ : the difference decays like $e^{- 1000 t}$ and is gone by $t = 0.01$ . After that the answer changes on a time scale of about 1, but an explicit method still has to take steps shorter than $2/ λ = 0.002$ , or the fast component, which is not there any more, grows back from rounding and truncation errors.

The Python reference, posts/refs/forward-euler-blows-up-past-h-equals-2-over-lambda.py, runs the same steps separately (output beside it).

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

[1]lines 24–54

#: LAM  the rate at which other solutions are pulled onto cos t
LAM = 1000.0
#: T_END  the length of every run
T_END = 2.0

#| The right side f(t, y) = -lam (y - cos t) - sin t.
def rhs(t, y) {
    return -LAM * (y - cos(t)) - sin(t)
}

#| Forward Euler from y0 with step h for `steps` steps; y[n] is the value at t = n h.
def forward_euler(h, y0, steps) {
    y = zeros(steps + 1)
    y[0] = y0
    for n to steps {
        y[n + 1] = y[n] + h * rhs(n * h, y[n])
    }
    return y
}

#| Backward Euler from y0 with step h. The step's equation
#| y1 = y0 + h f(t1, y1) is linear in y1, so it is solved in closed form.
def backward_euler(h, y0, steps) {
    y = zeros(steps + 1)
    y[0] = y0
    for n to steps {
        t = (n + 1) * h
        y[n + 1] = (y[n] + h * LAM * cos(t) - h * sin(t)) / (1.0 + h * LAM)
    }
    return y
}

Why $2/ λ$

Let $y_{n}$ be the computed solution and $e_{n} = y_{n} - cos t_{n}$ its error. Forward Euler is $y_{n + 1} = y_{n} + h f (t_{n}, y_{n})$ , and since $f$ depends on $y$ through $- λ y$ only,

e_{n + 1} = (1 - hλ) e_{n} + τ_{n},

where $τ_{n}$ , the local truncation error, is what one step of the method gets wrong starting from the exact solution. It is about $h^{2} /2$ times the second derivative of $cos t$ , so here at most $2 \times 1 0^{- 6}$ . Every error made earlier is multiplied by $g = 1 - hλ$ at every later step. If $∣ g ∣ < 1$ old errors shrink and the total stays near $τ / (hλ)$ . If $∣ g ∣ > 1$ they grow geometrically, whatever their size to begin with. $∣1 - hλ ∣ < 1$ holds exactly when $0 < h < 2/ λ$ . Accuracy would allow a far larger step: the solution is $cos t$ .

Backward Euler, $y_{n + 1} = y_{n} + h f (t_{n + 1}, y_{n + 1})$ , gives $e_{n + 1} = (e_{n} + τ_{n}) / (1 + hλ)$ instead, with factor $1/ (1 + hλ) < 1$ for every $h > 0$ . That is A-stability: the method is stable for every step size on every decaying mode. The cost is a solve per step (here a division, since the equation is linear in $y$ ).

The two runs below start on the solution, $y (0) = 1$ , straddling the limit: $h = 0.0019$ ( $g = - 0.9$ ) and $h = 0.0021$ ( $g = - 1.1$ ). In the figure the run at $h = 0.0021$ follows $cos t$ until the truncation errors, multiplied by $- 1.1$ at each step, reach order 1; it leaves the window at $t = 0.29$ and is drawn up to that point. The run at $h = 0.0019$ lies on $cos t$ , and so does backward Euler's at $h = 0.01$ .

[2]lines 88–168

H_INSIDE = 0.0019
H_OUTSIDE = 0.0021
N_INSIDE = round(T_END / H_INSIDE)
N_OUTSIDE = round(T_END / H_OUTSIDE)
Y_INSIDE = forward_euler(H_INSIDE, 1.0, N_INSIDE)
Y_OUTSIDE = forward_euler(H_OUTSIDE, 1.0, N_OUTSIDE)
err_inside = abs(Y_INSIDE[N_INSIDE] - cos(N_INSIDE * H_INSIDE))
err_outside = abs(Y_OUTSIDE[N_OUTSIDE] - cos(N_OUTSIDE * H_OUTSIDE))
# the error's growth per step over the second half of the outside run
HALF = N_OUTSIDE // 2
err_half = abs(Y_OUTSIDE[HALF] - cos(HALF * H_OUTSIDE))
growth = pow(err_outside / err_half, 1.0 / (N_OUTSIDE - HALF))
show err_inside
show err_outside
show growth

test "at h = 0.0019 forward Euler ends within 4e-7 of cos t (the Python reference: 3.956e-7)" {
    assert_near(err_inside, 3.956275553274047e-07, 1.0e-15)
}
test "at h = 0.0021 it is off by 2.67e33, as the Python reference finds" {
    assert_near(err_outside / 2.6731648134488684e+33, 1.0, 1.0e-9)
}
test "and its error grows by |1 - h lam| = 1.1 per step" {
    assert_near(growth, abs(1.0 - H_OUTSIDE * LAM), 1.0e-12)
}

H_IMPLICIT = 0.01
H_COARSE = 0.1
Y_IMPLICIT = backward_euler(H_IMPLICIT, 1.0, round(T_END / H_IMPLICIT))
Y_COARSE = backward_euler(H_COARSE, 1.0, round(T_END / H_COARSE))
err_implicit = abs(Y_IMPLICIT[len(Y_IMPLICIT) - 1] - cos(T_END))
err_coarse = abs(Y_COARSE[len(Y_COARSE) - 1] - cos(T_END))
show err_implicit
show err_coarse

test "backward Euler at h = 0.01, five times the explicit limit, ends 2.06e-6 from cos 2" {
    assert_near(err_implicit, 2.060996309838359e-06, 1.0e-15)
}
test "and at h = 0.1, fifty times the limit, 1.92e-5" {
    assert_near(err_coarse, 1.922824913669352e-05, 1.0e-15)
}

# the curves: forward Euler on both sides of the limit, backward Euler at 5x it
T_IN = zeros(N_INSIDE + 1)
for n to N_INSIDE + 1 {
    T_IN[n] = n * H_INSIDE
}
# the unstable run is drawn until it first leaves the window |y| <= Y_WINDOW
Y_WINDOW = 1.5
n_leave = N_OUTSIDE
for n from 0 to N_OUTSIDE + 1 reversed {
    if abs(Y_OUTSIDE[n]) > Y_WINDOW {
        n_leave = n
    }
}
T_OUT = zeros(n_leave + 1)
Y_OUT_SHOWN = zeros(n_leave + 1)
for n to n_leave + 1 {
    T_OUT[n] = n * H_OUTSIDE
    Y_OUT_SHOWN[n] = Y_OUTSIDE[n]
}
show n_leave * H_OUTSIDE

test "the unstable run first leaves |y| <= 1.5 at step 139, t = 0.2919, as in the Python reference" {
    assert_near(n_leave * H_OUTSIDE, 0.2919, 1.0e-12)
}
T_IMP = zeros(len(Y_IMPLICIT))
for n to len(Y_IMPLICIT) {
    T_IMP[n] = n * H_IMPLICIT
}
N_CURVE = 201
T_CURVE = linspace(0.0, T_END, N_CURVE)
figure("y' = -1000 (y - cos t) - sin t, from y(0) = 1", "t", "y", 0.0, T_END, -Y_WINDOW, Y_WINDOW)
plot(T_CURVE, cos(T_CURVE), "cos t, exact")
plot(T_IN, Y_INSIDE, "forward Euler, h = 0.0019")
plot(T_OUT, Y_OUT_SHOWN, "forward Euler, h = 0.0021")
plot(T_IMP, Y_IMPLICIT, "backward Euler, h = 0.01")
style("cos t, exact", "dash s6 thick nomarker")
style("forward Euler, h = 0.0019", "s1 nomarker")
style("forward Euler, h = 0.0021", "s4 nomarker")
style("backward Euler, h = 0.01", "dot s2 nomarker")

H_INSIDE = 0.0019
H_OUTSIDE = 0.0021
N_INSIDE = round(T_END / H_INSIDE)
N_OUTSIDE = round(T_END / H_OUTSIDE)
Y_INSIDE = forward_euler(H_INSIDE, 1.0, N_INSIDE)
Y_OUTSIDE = forward_euler(H_OUTSIDE, 1.0, N_OUTSIDE)
err_inside = abs(Y_INSIDE[N_INSIDE] - cos(N_INSIDE * H_INSIDE))
err_outside = abs(Y_OUTSIDE[N_OUTSIDE] - cos(N_OUTSIDE * H_OUTSIDE))
# the error's growth per step over the second half of the outside run
HALF = N_OUTSIDE // 2
err_half = abs(Y_OUTSIDE[HALF] - cos(HALF * H_OUTSIDE))
growth = pow(err_outside / err_half, 1.0 / (N_OUTSIDE - HALF))
show err_inside
show err_outside
show growth

test "at h = 0.0019 forward Euler ends within 4e-7 of cos t (the Python reference: 3.956e-7)" {
    assert_near(err_inside, 3.956275553274047e-07, 1.0e-15)
}
test "at h = 0.0021 it is off by 2.67e33, as the Python reference finds" {
    assert_near(err_outside / 2.6731648134488684e+33, 1.0, 1.0e-9)
}
test "and its error grows by |1 - h lam| = 1.1 per step" {
    assert_near(growth, abs(1.0 - H_OUTSIDE * LAM), 1.0e-12)
}

H_IMPLICIT = 0.01
H_COARSE = 0.1
Y_IMPLICIT = backward_euler(H_IMPLICIT, 1.0, round(T_END / H_IMPLICIT))
Y_COARSE = backward_euler(H_COARSE, 1.0, round(T_END / H_COARSE))
err_implicit = abs(Y_IMPLICIT[len(Y_IMPLICIT) - 1] - cos(T_END))
err_coarse = abs(Y_COARSE[len(Y_COARSE) - 1] - cos(T_END))
show err_implicit
show err_coarse

test "backward Euler at h = 0.01, five times the explicit limit, ends 2.06e-6 from cos 2" {
    assert_near(err_implicit, 2.060996309838359e-06, 1.0e-15)
}
test "and at h = 0.1, fifty times the limit, 1.92e-5" {
    assert_near(err_coarse, 1.922824913669352e-05, 1.0e-15)
}

# the curves: forward Euler on both sides of the limit, backward Euler at 5x it
T_IN = zeros(N_INSIDE + 1)
for n to N_INSIDE + 1 {
    T_IN[n] = n * H_INSIDE
}
# the unstable run is drawn until it first leaves the window |y| <= Y_WINDOW
Y_WINDOW = 1.5
n_leave = N_OUTSIDE
for n from 0 to N_OUTSIDE + 1 reversed {
    if abs(Y_OUTSIDE[n]) > Y_WINDOW {
        n_leave = n
    }
}
T_OUT = zeros(n_leave + 1)
Y_OUT_SHOWN = zeros(n_leave + 1)
for n to n_leave + 1 {
    T_OUT[n] = n * H_OUTSIDE
    Y_OUT_SHOWN[n] = Y_OUTSIDE[n]
}
show n_leave * H_OUTSIDE

test "the unstable run first leaves |y| <= 1.5 at step 139, t = 0.2919, as in the Python reference" {
    assert_near(n_leave * H_OUTSIDE, 0.2919, 1.0e-12)
}
T_IMP = zeros(len(Y_IMPLICIT))
for n to len(Y_IMPLICIT) {
    T_IMP[n] = n * H_IMPLICIT
}
N_CURVE = 201
T_CURVE = linspace(0.0, T_END, N_CURVE)
figure("y' = -1000 (y - cos t) - sin t, from y(0) = 1", "t", "y", 0.0, T_END, -Y_WINDOW, Y_WINDOW)
plot(T_CURVE, cos(T_CURVE), "cos t, exact")
plot(T_IN, Y_INSIDE, "forward Euler, h = 0.0019")
plot(T_OUT, Y_OUT_SHOWN, "forward Euler, h = 0.0021")
plot(T_IMP, Y_IMPLICIT, "backward Euler, h = 0.01")
style("cos t, exact", "dash s6 thick nomarker")
style("forward Euler, h = 0.0019", "s1 nomarker")
style("forward Euler, h = 0.0021", "s4 nomarker")
style("backward Euler, h = 0.01", "dot s2 nomarker")

Measuring the factor

To measure $g$ without the truncation error in the way, the last cell runs each method twice with the same step, from $y (0) = 2$ and from $y (0) = 1$ , and subtracts. The $τ_{n}$ are the same in both runs and cancel, so the difference after $M$ steps is $g^{M}$ times the starting difference of 1. $∣ d_{M} ∣^{1/ M}$ with $M = 4$ is the measured factor. The figure shows it at 36 step sizes from $1 0^{- 4}$ to $1 0^{- 1}$ , against the two formulas; the vertical line is $h = 2/ λ$ and the horizontal one is $∣ g ∣ = 1$ .

Forward Euler's factor falls to zero at $h = 1/ λ$ , where one step lands exactly on the solution, and crosses 1 at $h = 2/ λ$ . Measured and predicted agree to $2 \times 1 0^{- 14}$ . Backward Euler's factor falls steadily toward 0; it agrees with $1/ (1 + hλ)$ to $2 \times 1 0^{- 9}$ relative, the rounding left in a difference of $1 0^{- 8}$ after four steps at $h = 0.1$ .

[3]lines 189–250

#: M_STEPS  steps over which the factor is measured
M_STEPS = 4
#: N_H  step sizes in the sweep, spaced evenly in log h from 1e-4 to 1e-1
N_H = 36
HS = zeros(N_H)
G_FWD = zeros(N_H)
G_BWD = zeros(N_H)
P_FWD = zeros(N_H)
P_BWD = zeros(N_H)
for k to N_H {
    h = pow(10.0, -4.0 + 3.0 * k / (N_H - 1.0))
    HS[k] = h
    d_fwd = forward_euler(h, 2.0, M_STEPS)[M_STEPS] - forward_euler(h, 1.0, M_STEPS)[M_STEPS]
    d_bwd = backward_euler(h, 2.0, M_STEPS)[M_STEPS] - backward_euler(h, 1.0, M_STEPS)[M_STEPS]
    G_FWD[k] = pow(abs(d_fwd), 1.0 / M_STEPS)
    G_BWD[k] = pow(abs(d_bwd), 1.0 / M_STEPS)
    P_FWD[k] = abs(1.0 - h * LAM)
    P_BWD[k] = 1.0 / (1.0 + h * LAM)
}

test "forward Euler's measured factor is |1 - h lam| at all 36 step sizes, to 1e-13" {
    for k to N_H {
        assert_near(G_FWD[k], P_FWD[k], 1.0e-13 * max(1.0, P_FWD[k]))
    }
}
test "it is below 1 exactly for h < 2/lam = 0.002" {
    for k to N_H {
        assert_true((G_FWD[k] < 1.0) == (HS[k] < 2.0 / LAM))
    }
}
test "backward Euler's is 1/(1 + h lam), below 1 at every step size, to 1e-8 relative" {
    for k to N_H {
        assert_near(G_BWD[k] / P_BWD[k], 1.0, 1.0e-8)
        assert_lt(G_BWD[k], 1.0)
    }
}

N_FINE = 400
H_FINE = zeros(N_FINE)
F_FINE = zeros(N_FINE)
B_FINE = zeros(N_FINE)
for k to N_FINE {
    H_FINE[k] = pow(10.0, -4.0 + 3.0 * k / (N_FINE - 1.0))
    F_FINE[k] = abs(1.0 - H_FINE[k] * LAM)
    B_FINE[k] = 1.0 / (1.0 + H_FINE[k] * LAM)
}
figure("Growth factor per step, measured and predicted ($\\lambda = 1000$)", "step h", "|g|", 1.0e-4, 0.1, 1.0e-3, 200.0)
xscale("log")
yscale("log")
plot(H_FINE, F_FINE, "$|1 - h\\lambda|$")
plot(H_FINE, B_FINE, "$1/(1 + h\\lambda)$")
plot(HS, G_FWD, "forward Euler, measured")
plot(HS, G_BWD, "backward Euler, measured")
hline("|g| = 1", 1.0)
vline("h = 2/lambda", 2.0 / LAM)
style("$|1 - h\\lambda|$", "dash s4 nomarker")
style("$1/(1 + h\\lambda)$", "dash s2 nomarker")
style("forward Euler, measured", "s4 noline ring")
style("backward Euler, measured", "s2 noline ring")

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors

#: M_STEPS  steps over which the factor is measured
M_STEPS = 4
#: N_H  step sizes in the sweep, spaced evenly in log h from 1e-4 to 1e-1
N_H = 36
HS = zeros(N_H)
G_FWD = zeros(N_H)
G_BWD = zeros(N_H)
P_FWD = zeros(N_H)
P_BWD = zeros(N_H)
for k to N_H {
    h = pow(10.0, -4.0 + 3.0 * k / (N_H - 1.0))
    HS[k] = h
    d_fwd = forward_euler(h, 2.0, M_STEPS)[M_STEPS] - forward_euler(h, 1.0, M_STEPS)[M_STEPS]
    d_bwd = backward_euler(h, 2.0, M_STEPS)[M_STEPS] - backward_euler(h, 1.0, M_STEPS)[M_STEPS]
    G_FWD[k] = pow(abs(d_fwd), 1.0 / M_STEPS)
    G_BWD[k] = pow(abs(d_bwd), 1.0 / M_STEPS)
    P_FWD[k] = abs(1.0 - h * LAM)
    P_BWD[k] = 1.0 / (1.0 + h * LAM)
}

test "forward Euler's measured factor is |1 - h lam| at all 36 step sizes, to 1e-13" {
    for k to N_H {
        assert_near(G_FWD[k], P_FWD[k], 1.0e-13 * max(1.0, P_FWD[k]))
    }
}
test "it is below 1 exactly for h < 2/lam = 0.002" {
    for k to N_H {
        assert_true((G_FWD[k] < 1.0) == (HS[k] < 2.0 / LAM))
    }
}
test "backward Euler's is 1/(1 + h lam), below 1 at every step size, to 1e-8 relative" {
    for k to N_H {
        assert_near(G_BWD[k] / P_BWD[k], 1.0, 1.0e-8)
        assert_lt(G_BWD[k], 1.0)
    }
}

N_FINE = 400
H_FINE = zeros(N_FINE)
F_FINE = zeros(N_FINE)
B_FINE = zeros(N_FINE)
for k to N_FINE {
    H_FINE[k] = pow(10.0, -4.0 + 3.0 * k / (N_FINE - 1.0))
    F_FINE[k] = abs(1.0 - H_FINE[k] * LAM)
    B_FINE[k] = 1.0 / (1.0 + H_FINE[k] * LAM)
}
figure("Growth factor per step, measured and predicted ($\\lambda = 1000$)", "step h", "|g|", 1.0e-4, 0.1, 1.0e-3, 200.0)
xscale("log")
yscale("log")
plot(H_FINE, F_FINE, "$|1 - h\\lambda|$")
plot(H_FINE, B_FINE, "$1/(1 + h\\lambda)$")
plot(HS, G_FWD, "forward Euler, measured")
plot(HS, G_BWD, "backward Euler, measured")
hline("|g| = 1", 1.0)
vline("h = 2/lambda", 2.0 / LAM)
style("$|1 - h\\lambda|$", "dash s4 nomarker")
style("$1/(1 + h\\lambda)$", "dash s2 nomarker")
style("forward Euler, measured", "s4 noline ring")
style("backward Euler, measured", "s2 noline ring")

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors

Why 2/λ

Measuring the factor

Why $2/ λ$