2. Time integration

Can a decaying quantity oscillate numerically?

A positive quantity obeying $x'=-kx$ , with $k>0$ , never changes sign. Yet a numerical approximation can alternate between positive and negative values. This lesson explains why before running the same decay model with an adaptive method.

The equation describes a normalized temperature difference in a simple cooling model and appears in first-order circuit dynamics. Its physical derivation is in the ODE lesson. Here we need only $x(0)=1$ and the exact solution $x(t)=e^{-kt}$ .

Replace a derivative by a step

Approximate the derivative at $t_n$ by a forward difference:

\frac{x_{n+1}-x_n}{\Delta t}=-kx_n, \qquad x_{n+1}=(1-k\Delta t)x_n.

This is forward Euler. Its amplification factor is $g=1-z$ , where $z=k\Delta t$ . The exact factor is $e^{-z}$ . Taylor expansion gives $e^{-z}=1-z+z^2/2+O(z^3)$ : the difference per step starts at second order, but accumulated error over a fixed interval is generally first order.

For positive $z$ , decay in magnitude requires $|1-z|<1$ , or $0<z<2$ . Preserving nonnegative values further requires $z\leq1$ . Stability and positivity are different requirements.

$z$	Forward Euler factor	Predicted behavior
$0.5$	$0.5$	Positive decay, faster than exact decay
$1.5$	$-0.5$	Alternating signs with decaying magnitude
$2$	$-1$	Persistent sign alternation
$2.5$	$-1.5$	Growing oscillation

These entries follow from the recurrence, not from an Eqiora run. Starting from 1 with $z=2.5$ , the first three steps are $-1.5$ , $2.25$ , and $-3.375$ . The differential equation remains a decay equation throughout this failure.

Move the slope to the new time

Backward Euler instead uses $-kx_{n+1}$ :

\frac{x_{n+1}-x_n}{\Delta t}=-kx_{n+1}, \qquad x_{n+1}=\frac{x_n}{1+z}.

Its factor lies between zero and one for every positive $z$ . A very large step is therefore stable for this problem, but need not be accurate. At $z=1$ , one step gives $1/2$ , compared with $e^{-1}\approx0.367879$ .

For $k=1\,\mathrm{s}^{-1}$ and $T=1\,\mathrm{s}$ , taking $N$ equal steps gives the following formulas, which you can evaluate without implementing a solver:

x_N^{\mathrm{FE}}=(1-1/N)^N, \qquad x_N^{\mathrm{BE}}=(1+1/N)^{-N}.

Both tend to $e^{-1}$ . Their errors approach zero from opposite sides. Expand the logarithms to see the leading error proportional to $1/N$ .

Express the model once

// A minimal implicit ODE: x decays at the rate supplied by the user.
model decay() {
  // `x` is dimensionless and starts at 1. `rate` has inverse-time units.
  state x: 1;
  initial { x = 1; }
  parameter rate: 1 / s = 1;

  // Eqiora writes the evolution law as a residual equal to zero.
  relation flow {
    derivative(x) + rate * x = 0;
  }
}

There is no time-step size in this source. The rate and initial value describe the mathematical problem; the method describes how to approximate it.

Run the adaptive method

Use the files and environment from Get started:

uv run python run.py

"""Run the decay model shared by Get started and the ODE lesson."""

from __future__ import annotations

import argparse
from pathlib import Path

import eqiora


OUTPUT_TIMES_S = (0.25, 0.5, 1.0)


def solve(source_path: Path) -> tuple[tuple[float, float], ...]:
    """Return the requested (time in seconds, dimensionless value) samples."""

    model = eqiora.compile(path=source_path)
    field = model.field(model.field_ids[0])
    plan = eqiora.resolve(
        model,
        temporal=eqiora.time.Tsitouras45(
            initial_step_s=0.01,
            relative_tolerance=1.0e-9,
            absolute_tolerances={(field, 0, 0, False): 1.0e-11},
        ),
    )
    result = eqiora.run(
        plan,
        state=eqiora.State.initial(plan),
        until_s=OUTPUT_TIMES_S[-1],
        output_times_s=OUTPUT_TIMES_S,
    )
    series = result.series(field)
    times = series.time.numpy()
    values = series.values.numpy()
    return tuple(
        (float(time_s), float(value))
        for time_s, value in zip(times, values, strict=True)
    )


def main() -> None:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("model", nargs="?", type=Path, default=Path("decay.eqi"))
    args = parser.parse_args()
    for time_s, value in solve(args.model):
        print(f"t={time_s:.2f}, x={value:.10f}")


if __name__ == "__main__":
    main()

This program uses Tsitouras45, an explicit embedded Runge–Kutta method. The pair supplies two accuracy estimates from related stage evaluations; their difference guides time-step adaptation. See Tsitouras’s original paper for the method’s construction.

initial_step_s starts the adaptive process; it does not fix every subsequent step. output_times_s requests observations and is also independent of the internal steps. Requesting more output samples does not make the integration more accurate.

Compare the three values with $e^{-t}$ . Then hold the model, interval, and observation times fixed and change the requested relative tolerance from $10^{-9}$ to $10^{-6}$ , together with the absolute tolerance from $10^{-11}$ to $10^{-8}$ . Record the actual errors; a local error-control parameter is not an exact bound on every reported global error.

The program wraps the compile, resolve, and run calls in a solve function, which you can reuse for repeated comparisons; the equation itself stays in .eqi. Open the script.

Recognize stiffness without confusing it with instability

If a system contains both $e^{-t}$ and $e^{-1000t}$ , the fast component may force a small stable explicit step long after it becomes physically small. That mismatch between time scales is a reason to investigate stiff methods. One scalar run at a single rate says little about a coupled stiff system. Note also that the Euler formulas above are worked by hand; the script uses Tsitouras45, which has its own stability limits.

Exercises

Compute both Euler endpoint formulas for $N=4,8,16$ . Divide successive absolute errors and compare with the predicted factor of two.
For $k=20\,\mathrm{s}^{-1}$ , derive the forward Euler stability and nonnegativity bounds on $\Delta t$ .
Add an output time at $0.75\,\mathrm{s}$ . Explain why this is a different experiment from tightening the integration tolerances.
Repeat the adaptive run with rate 2. State the new exact comparison before running; identify the time scale that changed.
Why can a stable backward Euler answer still be unsuitable for estimating a short transient’s peak?

Previous: Errors and residuals · Book map · Next: Weak forms and finite elements