A decaying quantity follows y′=−y, so its exact solution stays positive and tends smoothly to zero. Yet explicit Euler with a time step of 1.5 produces 1, −0.5, +0.25, −0.125, …: the computed result alternates between negative and positive while its magnitude shrinks. With a step of 2.5 the numerical magnitude grows without bound. The continuous physical model did not change—only the time step did.
This article derives the forward Euler update and its stability interval, separates stability from positivity and accuracy, calculates representative test cases, and explains why stiff systems and real nonnegative state variables require more than an arbitrary ‘small enough’ timestep.
Find the question you need
All chapters remain visible. Choose a mathematical derivation, worked example, model check or practice route.
Derive the Euler stability rule
Chapters 1–3. Start here.
Separate stability and accuracy
Chapters 4–6. Start here.
Check stiffness and alternatives
Chapters 7–10. Start here.
Practise and validate
Chapters 11–12. Start here.
Complete chapter index · Practice · Scientific Computing subject root
Chapter index
Derive the Euler stability rule · Chapters 1–3
Separate stability and accuracy · Chapters 4–6
Check stiffness and alternatives · Chapters 7–10
Practise and validate · Chapters 11–12
Consider the differential equation dy/dt=−ky for a nonnegative quantity y and constant rate k>0. It can model an idealised radioactive inventory, a simple first-order chemical decay or a linear relaxation process over a suitable interval. The exact solution is y(t)=y₀e^(−kt) for initial amount y₀.
When y₀>0, the exact solution stays positive at every finite t and decreases monotonically toward zero. For k=1 per unit time and y₀=1, the value at t=2 is e^(−2)≈0.135335 and at t=3 it is e^(−3)≈0.049787.
These exact numbers give an independent standard for assessing a numerical method. A simulation that reports growing oscillations for this simple decay model is exhibiting a numerical pathology, not discovering a new physical process.
The test also isolates mathematics from measurement noise or uncertain parameters. Real systems can have changing k, sources and nonlinear interactions, but a method should first reproduce a simple case with a known analytic solution.
Euler’s method approximates the derivative over a step h by (yₙ₊₁−yₙ)/h and uses the derivative at the beginning of that step. Thus yₙ₊₁≈yₙ+h f(tₙ,yₙ), where f is the differential-equation right-hand side.
For the decay equation f(t,y)=−ky, substitute to obtain yₙ₊₁=(1−kh)yₙ. Each step multiplies the current estimate by the same factor q=1−kh. After n equal steps, yₙ=y₀(1−kh)ⁿ.
This update is easy to implement and often useful when h is sufficiently well chosen. But a multiplication factor becomes physically meaningful only after the step and model parameters are examined together. The dimensionless product kh, not h alone, controls the behaviour.
The units are a useful checkpoint: k has units inverse time and h has units time, so kh is dimensionless. If a student uses k measured per second and h recorded in minutes without conversion, the numerical stability decision can be wrong by a factor of sixty.
A numerical solution of the unforced decay test should have its perturbations shrink. For the repeated multiplier q=1−kh, the magnitude after n steps is |yₙ|=|y₀||q|ⁿ. It tends to zero when |1−kh|<1.
Solving −1<1−kh<1 gives 0<kh<2, or 0<h<2/k for positive k. This is the absolute-stability interval for explicit Euler applied to this linear negative-real test equation.
At kh=2, q=−1 and the magnitude does not decay: the numerical solution alternates indefinitely. At kh>2, |q|>1 and the magnitude grows, sometimes dramatically, even though the exact differential equation continues to decay.
This bound is a property of the numerical method on the test equation, not a universal time-step rule for every nonlinear physical system. It provides an important first diagnostic and illustrates how numerical dynamics can contradict the model being solved.
Two conditions are easy to conflate. For a positive starting amount, Euler estimates stay nonnegative and monotonically decrease if 0≤q≤1, equivalent to 0≤kh≤1. At kh=1, the next estimate becomes zero and stays there.
For 1<kh<2, q is negative but has magnitude less than one. The numerical trajectory alternates sign while approaching zero. It is stable in the absolute-error amplification sense for this linear model, but it violates the positivity expected of a decaying inventory or concentration.
At kh=1.5, q=−0.5. Starting at one, Euler produces 1, −0.5, +0.25, −0.125 and so forth. This is not a valid positive inventory at intermediate times, even though the oscillation amplitude decays.
An application may therefore require a stricter step-size restriction than stability alone. Preserving physical invariants such as nonnegativity, conservation or monotonicity is a separate numerical-design objective.
Fix k=1 and y₀=1. With h=0.5, q=0.5; over four steps to t=2, Euler produces y₄=0.5⁴=0.0625. The exact result at t=2 is e^(−2)≈0.135335. The numerical value is positive and stable but not highly accurate.
With h=1, q=0. Euler gives zero after the first step. This remains stable in the narrow multiplier sense and nonnegative, but prematurely wipes out a smooth exponentially decaying quantity.
With h=1.5, q=−0.5, producing an alternating sequence whose magnitude shrinks. With h=2, q=−1, so the values are 1,−1,1,−1,… and do not approach zero. With h=2.5, q=−1.5 and the magnitude grows like 1.5ⁿ.
These five cases use one differential equation and one starting value. The change in qualitative behaviour belongs entirely to the discretisation method and timestep, not to any new source term or changed physical decay constant.
| k h | Euler multiplier 1−kh | Numerical behaviour | Nonnegative? |
|---|---|---|---|
| 0.5 | 0.5 | Stable monotone decay | Yes |
| 1.0 | 0 | Becomes zero in one step | Yes |
| 1.5 | −0.5 | Stable alternating decay | No |
| 2.0 | −1 | Persistent oscillation | No |
| 2.5 | −1.5 | Growing oscillations | No |
At h=0.5 the explicit Euler result at t=2 is 0.0625, whereas the exact value is about 0.135335; the absolute error is roughly 0.072835. Reducing h to 0.25 gives eight steps and y₈=0.75⁸≈0.100113, reducing the error to about 0.035223.
Euler’s method has first-order global accuracy for a sufficiently smooth problem under suitable conditions: halving a small step tends to reduce the leading global error by about a factor of two. The exact ratio depends on the regime, and an unstable step does not become trustworthy through wishful convergence statements.
Even an absolutely stable step near kh=2 may make a poor approximation because the numerical solution can oscillate or decay at a wholly wrong rate. It is not enough to check that the plotted values eventually shrink.
An accuracy assessment compares against a reference solution, a verified higher-order or smaller-step computation, or a carefully controlled error estimate. A stability assessment asks whether perturbations grow. Both are necessary, and neither automatically supplies the other.
Over one step h, the exact solution multiplies by e^(−kh). Euler uses 1−kh. The Taylor expansion e^(−kh)=1−kh+(kh)²/2−(kh)³/6+… shows why Euler becomes accurate as kh becomes small.
For h=0.5 and k=1, the exact one-step factor is e^(−0.5)≈0.60653, while Euler uses 0.5. After four steps, the differences compound: 0.60653⁴=e^(−2)≈0.135335, versus 0.5⁴=0.0625.
The difference is a time-discretisation error. It is not the same mechanism as floating-point cancellation in a numerical derivative or summation order loss, although these can also affect a computer implementation.
This distinction matters when debugging. If results are wrong because h is too large for Euler’s approximation, switching from binary64 to arbitrary precision will calculate the wrong discrete recurrence more accurately. The method and step must be addressed.
Now imagine two independent components, u′=−u and v′=−100v. The slow component has rate one, while the fast component has rate 100 per time unit. Both are physically stable and decay exponentially in the exact solution.
Explicit Euler uses multipliers 1−h for u and 1−100h for v. To be absolutely stable for both, the faster mode imposes h<2/100=0.02. If positivity is also required for both, the stricter bound is h≤0.01.
A step h=0.1 might look modest compared with the slow process, but it gives a fast-mode multiplier of 1−10=−9. The simulated fast component grows explosively in alternating sign, even as the real fast component dies away.
This separation of timescales is the essence of a stiffness issue: stability of an explicit integrator can force steps far smaller than the timescale on which the user wishes to observe the slow solution. A more suitable numerical method may be needed.
Backward, or implicit, Euler evaluates the derivative at the new value: yₙ₊₁=yₙ−hk yₙ₊₁. Rearranging gives yₙ₊₁=yₙ/(1+kh), so its multiplier is 1/(1+kh). For k,h>0 it is positive and lies between zero and one.
Thus the backward Euler test solution is stable and positive for any positive h on this linear decay equation. At h=2.5 with k=1, the one-step estimate is 1/3.5≈0.285714, while the exact value at t=2.5 is e^(−2.5)≈0.082085.
The implicit method avoids the explicit scheme’s divergent oscillation, but the large-step value can be inaccurate. ‘Unconditionally stable’ in this test-equation sense does not mean ‘correct at any step size’.
Implicit methods can require solving algebraic or nonlinear equations at each step for general systems. They have additional computational costs and error controls, so the method choice should be guided by stiffness, conservation and required accuracy.
A concentration, stored population or inventory should not become negative merely because a numerical update overshot a stable decay state. A negative Euler value at kh=1.5 is an algorithmic artefact, not a prediction of the original model.
Clipping the result to zero after every step can prevent a negative display but changes the numerical method and can introduce its own mass-balance or accuracy problems. It is not generally a substitute for choosing a positivity-preserving method or appropriate timestep.
For complex reaction networks, a variable may interact with conservation constraints and nonlinear rates. An independently stable linear test is a beginning, not proof that the complete numerical model respects all physical requirements.
Document any positivity correction, limiter or projection operation. A computed trace should say whether values came directly from an integration rule or from postprocessing that altered them.
Exercise A: A decay equation has k=4 per minute. What is explicit Euler’s absolute-stability interval? 0<h<0.5 minute, since h<2/4. To retain nonnegative monotone values, use the stricter bound h≤0.25 minute.
Exercise B: For k=1 and h=1.8, the multiplier is −0.8. Is it absolutely stable? Yes, since |−0.8|<1. Does it respect nonnegative decay? No, because consecutive estimates alternate sign.
Exercise C: With k=10 and h=0.25, the multiplier is 1−2.5=−1.5. The numerical magnitude grows, contradicting the exact exponential decay. Reducing h below 0.2 is necessary for absolute stability on this test equation.
Exercise D: An implicit Euler run uses h=5 and reports a positive, decaying sequence. Has accuracy been proved? No. Check the exact or reference solution; stability alone does not quantify its approximation error.
Begin with the mathematical model, units, initial conditions and exact properties such as positivity, monotonic decay and conserved quantities. Use a test equation with a known solution to validate the code’s update formula and timestep interpretation.
Second, examine the amplification factors or stability properties of the chosen method for the relevant modes. In an explicit method, the fastest mode may determine a much smaller permissible h than the slow phenomenon of primary interest.
Third, perform a step-refinement study and compare with a trustworthy reference. Check physical invariants separately from numerical stability, and report whether positivity clipping, adaptive timestep control or implicit solves were used.
Continue with The Gold Standard of Scientific Computing, the finite-difference step-size problem and compensated summation. Different numerical algorithms fail for different reasons, so the correct remedy must match the actual error mechanism.
Applied numerical stability lab
Reproduce the multiplier classification by hand
Write q=1−kh and construct a table for kh=0.5,1,1.5,2 and 2.5. The multipliers are respectively 0.5,0,−0.5,−1 and −1.5. The categories are monotone positive decay, immediate zero, alternating decay, persistent oscillation and alternating divergence.
This compact table is more informative than trying a series of arbitrary time-step numbers. The dimensionless product kh carries the entire stability condition for the idealised single-rate test.
Ask a learner to identify whether each trace respects the physical sign and whether it converges numerically to zero. Those answers show that positivity and asymptotic stability are separate tests.
Use two mode rates to recognise stiffness
For u′=−u and v′=−100v, choose h=0.01. Explicit Euler gives multipliers 0.99 for u and 0 for v, so the scheme is stable and nonnegative for these two components. With h=0.015, the fast multiplier becomes −0.5: still decaying in absolute size but alternating sign.
At h=0.025 the fast multiplier becomes −1.5 and diverges, even though the slow multiplier 0.975 is benign. This demonstrates why a timestep selected from the slow dynamics may fail a fast component.
A student should use the maximum relevant k when deriving a common explicit-Euler stability bound for independent decays, then check whether that bound is small enough to make the calculation impractical.
Find an exact-solution check over the same time interval
For k=1 and t=2, compute the exact result e^(−2)≈0.135335. With h=0.5 use four Euler steps, yielding 0.0625. With h=0.25 use eight steps, yielding approximately 0.100113. Both are stable and positive but differ noticeably in accuracy.
Now reduce h further and track convergence toward the exact value. The expectation of first-order global convergence should emerge once the step is sufficiently small and other arithmetic errors are controlled.
This is a useful experimental distinction: compare estimates at the same final physical time. Comparing one Euler step of h=0.5 with one step of h=0.25 would mean comparing two different times, not measuring time-step accuracy.
Interpret what an implicit solver really changes
Backward Euler solves a different discrete update equation. On the decay test its multiplier 1/(1+kh) remains positive for all positive kh, making it valuable for stiff decay dynamics.
However, its large-step trajectory still deviates from the exact exponential. Implicit stability controls numerical growth, whereas model fit, time resolution and global error remain separate concerns.
Use this comparison to teach that a stable computer simulation is not automatically an accurate scientific model. Reporting the integrator and tolerance settings is part of explaining how the numerical answer was produced.
Frequently asked questions
Can a stable numerical solution be negative?
Yes. For 1<kh<2, explicit Euler’s magnitude decays but successive values alternate in sign, violating nonnegative-decay behaviour.
Does backward Euler eliminate all error?
No. It is stable on this linear decay test for any positive step, but its large-step approximation can still be inaccurate.
Why can the fastest mode control the timestep?
For multiple decay rates, the largest positive k gives the most restrictive explicit Euler bound h<2/k for absolute stability.
Sources and connected reading
All examples are original idealised decay systems and exact or calculated discretisation results. Real engineering simulations need model-specific checks.
- University of Cambridge — Euler’s method and stability on exponential decay
- The Gold Standard of Scientific Computing — subject root
- The Gold Standard of Scientific Computing — finite differences and step size
Advanced stability lab: convergence, complex modes and adaptive steps
A step-refinement table that shows first-order convergence
Hold the final physical time at t=2 and the initial value at y₀=1 for k=1. With h=0.5, use four explicit Euler steps and obtain 0.062500. With h=0.25, use eight steps and obtain 0.100113. With h=0.125, use sixteen steps and obtain about 0.118067.
The exact value e^(−2)≈0.135335. The corresponding absolute errors are approximately 0.072835, 0.035222 and 0.017268. Halving h nearly halves the error once the method is in its expected convergence regime, consistent with a first-order global method.
At h=0.05 with forty steps, Euler gives approximately 0.128512 and error 0.006823. The improvement is real, but the value remains visibly below the exact one. Stability of all four tested steps did not make any finite-step result exactly correct.
This check must compare the same final t. Computing four steps for every h would end at different times and would not isolate time-step accuracy. An educational table should include h, step count, final time, numerical value and reference error together.
Extend the stability idea to oscillatory complex modes
The simple real decay test is part of a larger analysis. For a linear test equation y′=λy with complex λ, explicit Euler multiplies by 1+hλ at each step. Absolute stability requires |1+hλ|<1. In the complex hλ plane this is the open disk of radius one centred at −1.
Suppose λ=−1+2i per time unit. The exact solution decays in magnitude because the real part of λ is negative, while its imaginary part represents oscillation. With h=0.2, Euler’s multiplier is 0.8+0.4i and its magnitude is √0.8≈0.8944, so this mode is numerically damped.
Increase h to 0.5 and the multiplier becomes 0.5+1.0i, whose magnitude is √1.25≈1.118. Its numerical amplitude grows even though the exact solution’s amplitude still decays. A timestep acceptable for one real decay constant need not be acceptable for a different oscillatory mode.
This extension shows why numerical stability is connected to a method’s stability region and the eigenvalues of a linearised system. In complex physical simulations, a simple positive decay-rate checklist is necessary background but not the whole stability analysis.
Separate adaptive accuracy control from positivity guarantees
An adaptive ODE solver can adjust h based on an estimated local error. Smaller steps are used where the model changes rapidly and larger steps where the solution varies smoothly. This can be efficient, but error control and invariant preservation are distinct requirements.
A solver might satisfy a relative and absolute error tolerance yet briefly produce a small negative concentration if its numerical method does not preserve nonnegativity. Depending on the model, that may be harmless for a display, or it may invalidate a subsequent reaction-rate expression.
Likewise a method may preserve positive values but deliver a poor approximation with a large h. Our backward Euler example at h=2.5 stays positive and stable but is far from e^(−2.5). A full validation plan checks both error estimates and physical constraints.
When software clips a negative value to zero or renormalises a mass fraction, that postprocessing should be stated. Such corrections can change conservation laws and should not silently replace a more appropriate integrator or a well-justified step selection.
A self-contained classroom algorithm check
Write the exact model as y(t)=exp(−kt) for y₀=1. Then write the forward Euler recurrence y←y(1−kh). Use k=1 and test h=0.5, 1, 1.5, 2 and 2.5 for a handful of steps. Label each result as positive decay, alternating decay, persistent oscillation or growing oscillation.
Next calculate the stability bound 0<kh<2 and the stricter nonnegative monotone condition 0<kh≤1. Ask students to predict the result category before running the loop. This sequence connects algebra to the computational evidence and prevents the plot from becoming the only explanation.
Then compare one forward Euler estimate with the exact solution at the same time and repeat after halving h. The student should see that the improved stable estimate still has nonzero approximation error. They should identify stability, accuracy and positivity as three different claims.
Finally repeat with two independent decay rates, k₁=1 and k₂=100, on one shared grid. The fast mode controls explicit stability. The exercise demonstrates why a simulation can look well behaved in one plotted variable while another hidden variable has already become numerically unstable.
A decision rule for selecting a timestep
A practical timestep decision should start with the fastest significant decay rate, not the slowest plotted curve. For independent real decay modes, choose h below 2/k_max for absolute stability, and below or equal to 1/k_max if the method must keep every nonnegative component nonnegative in the simple explicit update. Then perform an accuracy-refinement study; satisfying the inequality alone is not enough.
If h would need to be extremely small to handle a fast mode that is not itself of scientific interest, consider an appropriate implicit or specialised stiff integrator. Record its error tolerances and any event handling. The goal is to preserve the physical model and obtain reliable results efficiently, rather than to make an unstable graph appear smooth.
Time units remain a basic audit item: a k stated per second and h stated in minutes must be converted before applying the product kh. A correct stability derivation with inconsistent units can still produce a computationally disastrous implementation.