Make Euler stable on a stiff equation

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

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

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.

Open in the workbench Download .vivOutputs below are from the site build’s run; Run all runs it again in your browser.

Look inside this program

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

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$ ).

Your task

Write these 2 functions:

decay_step(y, h): 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}).
system_step(y: Arr, h): 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.

Make these tests pass 4 of 4 fail

one step of 0.05 multiplies y by 1 / (1 + 2.5) = 0.285714 assert_near: |-1.5 - 0.285714| = 1.78571 > 1e-15 (fails)
twelve steps of 0.05 leave 2.9593e-7: small and positive, never negative assert_gt: -1.5 <= 0 (and 6 more) (fails)
a step 25 times the forward-Euler limit, h = 1, gives 1/51 = 0.0196 assert_near: |-49 - 0.0196078| = 49.0196 > 1e-15 (fails)
ten steps of 0.1 on the system end at numpy's (0.38554, 0.38554): the fast mode is gone assert_near: |9.04382e+19 - 0.385543| = 9.04382e+19 > 1e-13 (and 2 more) (fails)

Printed

slow_rate: -1

fast_rate: -1000

y: 129.746

lambda_h: 2.5

Your task

Solutions 0