Notebooks · posts/make-euler-stable-on-a-stiff-equation

Make Euler stable on a stiff equation

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

The solution of $y^{'} = - 50 y$ decays to zero, but forward Euler with step $h = 0.05$ multiplies $y$ by $1 - 50 h = - 1.5$ every step, and after 12 steps it is 129.7 instead of $1 0^{- 13}$ . Rewrite decay_step and system_step as backward Euler, which takes the slope at the end of the step instead of the start, so that no step size makes them grow.

y^{'} = - 50 y, h = 0.05 : y_{n + 1} = (1 - 50 h) y_{n} = - 1.5 y_{n}

A problem is stiff when it contains a rate much faster than the one you want to follow. Forward Euler, $y_{n + 1} = y_{n} + h f (y_{n})$ , is stable on $y^{'} = - λ y$ only when $∣1 - λh ∣ \leq 1$ , that is $h \leq 2/ λ = 0.04$ here, however slowly the solution itself is changing. Backward Euler solves $y_{n + 1} = y_{n} + h f (y_{n + 1})$ for $y_{n + 1}$ . For $y^{'} = - λ y$ that is one division, and the factor per step is $1/ (1 + λh)$ , which is between 0 and 1 for every $h > 0$ .

The second function steps a system $y^{'} = A y$ whose matrix has eigenvalues $- 1$ and $- 1000$ . The slow part is what matters over $0 \leq t \leq 1$ ; the fast part dies out in a thousandth of that, but it is what makes forward Euler with $h = 0.1$ end at $9 \times 1 0^{19}$ . Backward Euler needs $(I - h A) y_{n + 1} = y_{n}$ : a 2 by 2 linear solve per step.

The expected numbers come from posts/refs/make-euler-stable-on-a-stiff-equation.py, which runs backward Euler in numpy. At $h = 0.1$ it is not accurate, only stable: the slow mode shrinks by $(1/1.1)^{10} = 0.3855$ where the exact answer is $e^{- 1} = 0.3679$ , and the last test checks the stable answer ( $0.3855$ ).

Look inside this program

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

This is a problem: its tests fail until the functions in it are written. Edit a cell and run it to try; the cells below it show whether the tests pass.

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 33–99

LAMBDA = 50.0
# the system's matrix: eigenvalues -1 (along (1, 1)) and -1000 (along (1, -1))
A_SYS = [[-500.5, 499.5], [499.5, -500.5]]

#| One step of $y' = -\lambda y$ with $\lambda$ = LAMBDA. As written it is forward
#| Euler. Make it backward Euler: $y_{n+1} = y_n + h f(y_{n+1})$.
def decay_step(y, h) {
    return y + h * (-LAMBDA * y)
}

#| One step of $y' = Ay$ for the 2 by 2 matrix A_SYS, returning the new y. As
#| written it is forward Euler. Make it backward Euler: solve $(I - hA)\,y_{n+1} = y_n$.
def system_step(y: Arr, h) {
    y_new = zeros(2)
    for i to 2 {
        y_new[i] = y[i] + h * (A_SYS[i, 0] * y[0] + A_SYS[i, 1] * y[1])
    }
    return y_new
}

# A_SYS (1, 1) = -(1, 1) and A_SYS (1, -1) = -1000 (1, -1): its two rates
slow_rate = A_SYS[0, 0] + A_SYS[0, 1]
fast_rate = A_SYS[0, 0] - A_SYS[0, 1]
show slow_rate
show fast_rate

# ten steps of 0.1 take the system from t = 0 to t = 1
H_SYS = 0.1
N_SYS = 10
# 12 steps of 0.05 take y' = -50 y from t = 0 to t = 0.6
H = 0.05
NSTEP = 12
y = 1.0
for n to NSTEP {
    y = decay_step(y, H)
}
show y
# lambda h, the number backward Euler divides by: 1 / (1 + lambda h) per step
lambda_h = LAMBDA * H
show lambda_h

test "one step of 0.05 multiplies y by 1 / (1 + 2.5) = 0.285714" {
    assert_near(decay_step(1.0, 0.05), 0.2857142857142857, 1.0e-15)
}
test "twelve steps of 0.05 leave 2.9593e-7: small and positive, never negative" {
    z = 1.0
    for n to NSTEP {
        z = decay_step(z, H)
        assert_gt(z, 0.0)
    }
    assert_near(z, 2.959262343537004e-07, 1.0e-20)
}
test "a step 25 times the forward-Euler limit, h = 1, gives 1/51 = 0.0196" {
    assert_near(decay_step(1.0, 1.0), 0.0196078431372549, 1.0e-15)
}
test "ten steps of 0.1 on the system end at numpy's (0.38554, 0.38554): the fast mode is gone" {
    w = [2.0, 0.0]
    for n to N_SYS {
        w = system_step(w, H_SYS)
    }
    assert_near(w[0], 0.3855432894295282, 1.0e-13)
    assert_near(w[1], 0.38554328942952826, 1.0e-13)
    assert_near(w[0] - w[1], 0.0, 1.0e-13)
}

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

LAMBDA = 50.0
# the system's matrix: eigenvalues -1 (along (1, 1)) and -1000 (along (1, -1))
A_SYS = [[-500.5, 499.5], [499.5, -500.5]]

#| One step of $y' = -\lambda y$ with $\lambda$ = LAMBDA. As written it is forward
#| Euler. Make it backward Euler: $y_{n+1} = y_n + h f(y_{n+1})$.
def decay_step(y, h) {
    return y + h * (-LAMBDA * y)
}

#| One step of $y' = Ay$ for the 2 by 2 matrix A_SYS, returning the new y. As
#| written it is forward Euler. Make it backward Euler: solve $(I - hA)\,y_{n+1} = y_n$.
def system_step(y: Arr, h) {
    y_new = zeros(2)
    for i to 2 {
        y_new[i] = y[i] + h * (A_SYS[i, 0] * y[0] + A_SYS[i, 1] * y[1])
    }
    return y_new
}

# A_SYS (1, 1) = -(1, 1) and A_SYS (1, -1) = -1000 (1, -1): its two rates
slow_rate = A_SYS[0, 0] + A_SYS[0, 1]
fast_rate = A_SYS[0, 0] - A_SYS[0, 1]
show slow_rate
show fast_rate

# ten steps of 0.1 take the system from t = 0 to t = 1
H_SYS = 0.1
N_SYS = 10
# 12 steps of 0.05 take y' = -50 y from t = 0 to t = 0.6
H = 0.05
NSTEP = 12
y = 1.0
for n to NSTEP {
    y = decay_step(y, H)
}
show y
# lambda h, the number backward Euler divides by: 1 / (1 + lambda h) per step
lambda_h = LAMBDA * H
show lambda_h

test "one step of 0.05 multiplies y by 1 / (1 + 2.5) = 0.285714" {
    assert_near(decay_step(1.0, 0.05), 0.2857142857142857, 1.0e-15)
}
test "twelve steps of 0.05 leave 2.9593e-7: small and positive, never negative" {
    z = 1.0
    for n to NSTEP {
        z = decay_step(z, H)
        assert_gt(z, 0.0)
    }
    assert_near(z, 2.959262343537004e-07, 1.0e-20)
}
test "a step 25 times the forward-Euler limit, h = 1, gives 1/51 = 0.0196" {
    assert_near(decay_step(1.0, 1.0), 0.0196078431372549, 1.0e-15)
}
test "ten steps of 0.1 on the system end at numpy's (0.38554, 0.38554): the fast mode is gone" {
    w = [2.0, 0.0]
    for n to N_SYS {
        w = system_step(w, H_SYS)
    }
    assert_near(w[0], 0.3855432894295282, 1.0e-13)
    assert_near(w[1], 0.38554328942952826, 1.0e-13)
    assert_near(w[0] - w[1], 0.0, 1.0e-13)
}

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