4. Ordinary differential equations

State makes history matter

An ordinary differential equation relates quantities to derivatives with respect to one independent variable, here time. A first-order initial-value problem needs both an evolution relation and an initial state. Without the initial value, the relation describes a family of trajectories rather than one trajectory.

Learning outcomes

After this chapter, you should be able to:

identify state, parameter, derivative relation, and initial value;
derive the closed form for constant-rate scalar decay;
keep the time method and output schedule outside the Model;
run the example with the installed Python package and select a series by its FieldRef; and
diagnose a physically wrong sign that remains dimensionally consistent.

From a rate of cooling to a trajectory

Return to the bead in chapter 1. Its heat capacity $C$ stores energy, and its conductance $G$ transfers energy to a room at fixed $T_\infty$ . The balance $C\dot T+G(T-T_\infty)=0$ predicts a rate, but our thermometer asks for a temperature at a later time. Integrating the rate equation bridges those questions. Before calculating, predict three features: the bead stays between its initial temperature and room temperature, the difference decreases monotonically, and the cooling slows with time.

Normalize with $x=(T-T_\infty)/(T_0-T_\infty)$ and let $k=G/C$ . The resulting source contains just one state, one rate, and one initial value. It also describes other first-order relaxation systems after the physical quantities have been mapped correctly. A shared equation does not make their parameters interchangeable.

Derive the decay equation

Let dimensionless state $x(t)$ decay with constant rate $k>0$ :

\frac{\mathrm{d}x}{\mathrm{d}t} + kx = 0, \qquad x(0)=1.

Separating variables and integrating from $0$ to $t$ gives

\int_1^{x(t)}\frac{1}{x}\,\mathrm{d}x =-\int_0^t k\,\mathrm{d}t, \qquad \log x(t)=-kt,

so the solution to compare against is

x(t)=\exp(-kt).

The executable lesson uses $k=1\,\mathrm{s}^{-1}$ and observes $t=0.25$ , $0.5$ , and $1\,\mathrm{s}$ . They let us compare the trajectory at a quarter, half, and one characteristic time.

What the solution tells us before a run

The characteristic time is $\tau=1/k$ . At $t=\tau$ , the remaining fraction is $e^{-1}\approx0.367879$ , independent of the unit used to express time. The half-life is found from $e^{-kt_{1/2}}=1/2$ , hence $t_{1/2}=\log(2)/k$ . Repeated equal time intervals multiply $x$ by the same factor; they do not subtract equal amounts.

For $k=1\,\mathrm{s}^{-1}$ , the analytic values are:

Time (s)	$\exp(-kt)$	Interpretation
0.25	0.7788007831	about 78% remains
0.50	0.6065306597	about 61% remains
1.00	0.3678794412	one characteristic time

These numbers come from evaluating the analytic solution, not reading a solver trace. The curve satisfies both $x(0)=1$ and $\dot x=-kx$ ; checking the initial value matters because $2e^{-kt}$ satisfies the same differential equation.

Recover physical temperature using $T(t)=T_\infty+(T_0-T_\infty)e^{-kt}$ . With the illustrative bead’s $\tau=100\,\mathrm s$ , $T_0=340\,\mathrm K$ , and $T_\infty=300\,\mathrm K$ , the predicted temperature after 100 seconds is about $314.715\,\mathrm K$ . Changing to this rate requires changing rate to 0.01; the default script’s one-second window then covers only one hundredth of a characteristic time.

Model source

// 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;
  }
}

The initial value and rate belong to the Model. The adaptive method, its tolerances, the final time, and requested observations belong to the Realization or Run.

Execute the model

Follow Get started to install the published package and download the same decay.eqi and run.py used here. Set rate to 1, then run from their folder:

$ uv run python run.py

The script compiles, resolves, and runs the model:

"""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()

View the run script.

The program prints each time and computed value. With these settings, each value should be within $2\times10^{-8}$ of $\exp(-t)$ , far tighter than the rounded percentages in the table. If you change the model, the interval, or the integration settings, recompute the expected values and reconsider the tolerance.

Read the result

The result series, selected with the field handle, contains the three requested output times; the solver may take additional internal steps between them.

If the values differ from the closed form, inspect the equation sign, initial value, rate, and solver settings.

What can we reuse?

decay.eqi is already the complete mathematical description for this small problem. There is no need to hide a one-line evolution law behind a new component. Reuse its source and run script, then open their definitions when you need to change the law. In the circuits path, resistors and sources are larger reusable physical parts, and imports remove repetition after their equations have been derived.

Change an assumption before changing a solver

Suppose the room temperature becomes $T_\infty(t)$ . The energy balance remains $C\dot T+G(T-T_\infty(t))=0$ , but the old normalization no longer produces unforced decay. For $\theta=T-T_\infty(t)$ ,

\dot\theta+k\theta=-\dot T_\infty(t).

This follows by differentiating the temperature difference. It explains why a varying ambient temperature needs a forcing term: cooling a difference and cooling an absolute temperature are different statements when the reference moves. Work this extension on paper before editing source.

By contrast, recording a fourth output time leaves the differential equation unchanged. Tightening an integrator tolerance changes how accurately the trajectory is approximated. The numerical-simulation path examines that distinction through time-step error and convergence.

Deliberate failure: the wrong sign still compiles

Change the residual to:

relation flow {
  derivative(x) - rate * x = 0;
}

Both terms still have inverse-time dimension, so compilation can succeed. But the relation now implies $x(t)=\exp(+kt)$ : growth instead of decay, so the comparison with the decay solution fails. Fix the sign by revisiting the physical assumption; loosening the tolerance would only hide the error.

Exercises

Derive the solution for an arbitrary positive initial value $x(0)=x_0$ .
Predict the three observed values when rate is $2\,\mathrm{s}^{-1}$ , then state which part of the source and which comparison expression must change.
Run the wrong-sign mutant. Explain why dimensional checking accepts it and the closed-form comparison rejects it.
Add one requested output time without changing the Model. Explain why this changes a Run request rather than mathematical meaning.
List two reasons the three passing samples do not establish a convergence rate or suitability for stiff systems.

Exercise guidance

Exercise 1 gives $x=x_0e^{-kt}$ . For exercise 2, replace $t$ with $2t$ in the comparison; the source’s rate and the comparison must describe the same model. For exercise 5, several samples on one smooth trajectory do not reveal how error changes as resolution changes. A stiff system introduces competing time scales that this scalar example does not contain.

Inverse question. Given two positive observations $x(t_1)$ and $x(t_2)$ , derive $k=\log(x(t_1)/x(t_2))/(t_2-t_1)$ . Ask what happens if measurement noise makes the second value larger than the first. The inverse-problems path takes up the distinction between estimating a parameter and validating the assumed law.

Continue with stochastic evolution

The stochastic-process notebook authors Ornstein–Uhlenbeck and geometric Brownian motion from scratch, then executes individual paths and computes typed ensemble moments. It separates sampling uncertainty from time-discretization bias and compares coupled GBM paths with the exact Ito solution. Every equation and execution step appears in the notebook; it requires the development checkout installed in its kernel.

References

Arthur Mattuck, Haynes Miller, Jeremy Orloff, and John Lewis, Differential Equations, MIT OpenCourseWare 18.03SC (2011), Basic DE’s and separable equations, especially “The Most Important DE” and “Separation of Variables.” The exponential and its initial-value interpretation underlie this chapter’s comparison.

Previous: Algebraic relations and networks · Back to the series map · Next: Fields and spatial domains