Reliable Differential-Equation Computation, differential equations lesson

A solver can return a smooth curve and a tidy table of numbers for a problem nobody meant to pose. Each check in this lesson looks at a different part of the problem, the equation, the data, the interval, the error trend or the units, and none of them replaces the others. Throughout, $t$ is the independent variable and $y' = dy/dt$ . The running example is the dimensionless problem $y' = y^2$ with $y(0) = 1$ , whose solution $y = 1/(1 - t)$ exists for $t < 1$ , and the lesson uses it to show what each check can and cannot establish.

Key terms

Equation residual

Definition For a differential equation $y' = f(t, y)$ and a candidate function $u$ on an interval, the equation residual is

r(t) = u'(t) - f\bigl(t, u(t)\bigr).

Here $u$ is continuous and differentiable except possibly at finitely many points, and $r$ is computed wherever $u'$ exists. A solution has $r(t) = 0$ throughout its interval, and an approximation with residual $r$ solves the nearby equation $u' = f(t, u) + r(t)$ exactly. For output known only at grid points, $u$ is an interpolant, so the residual also reflects the interpolation, and differentiating noisy data can magnify errors. A small residual is not by itself a bound on the error. Suppose $\lvert r(t) \rvert \le \delta$ on $[t_0, T]$ , the initial error is $\varepsilon_0 = \lvert u(t_0) - y(t_0) \rvert$ , and $\lvert f(t, y_1) - f(t, y_2) \rvert \le L\lvert y_1 - y_2 \rvert$ with $L > 0$ on a region containing the graphs of $u$ and of the solution $y$ . Then the error $E(t) = \lvert u(t) - y(t) \rvert$ satisfies

\begin{aligned} E(t) &\le \varepsilon_0 e^{L(t - t_0)} \\[0.4em] &\quad + \frac{\delta}{L}\bigl(e^{L(t - t_0)} - 1\bigr) \end{aligned}

for $t_0 \le t \le T$ . The growing factor $e^{L(t - t_0)}$ shows how a sensitive problem can turn a small residual into a large error.

In plain words Put the computed curve back into the equation and see how badly it fails. Zero failure means a solution. A small failure means the curve solves an equation close to the intended one, and how far apart the two solutions drift depends on how sensitive the problem is.

Example For $u = 1/(1 - t)$ with $t < 1$ ,

u' = \frac{1}{(1 - t)^2} = u^2,

so the residual is zero. Two Euler steps with $h = 1/4$ give $Y_1 = 5/4$ and $Y_2 = 105/64$ , and joining the points by straight lines gives a polygon $p$ . On the first step $p(t) = 1 + t$ , so its residual $1 - (1 + t)^2$ falls from $0$ to $-9/16$ at $t = 1/4$ . On the second step the slope is $25/16$ , and the residual reaches $25/16 - (105/64)^2 = -4625/4096 \approx -1.13$ at $t = 1/2$ . For $0 \le t \le 1/2$ both graphs stay in the strip $1 \le y \le 2$ , where $\lvert y_1^2 - y_2^2 \rvert \le 4\lvert y_1 - y_2 \rvert$ , so the bound with $\varepsilon_0 = 0$ , $\delta \approx 1.13$ and $L = 4$ guarantees an error of at most about $1.8$ at $t = 1/2$ . The actual error is $2 - 105/64 = 23/64 \approx 0.36$ . The bound is valid but loose, and the large residual already warns that the steps are too coarse.

Initial and boundary checks

Definition A solution of an initial-value problem must satisfy the differential equation on an interval containing $t_0$ and also the initial condition $y(t_0) = y_0$ . A solution of a boundary-value problem must satisfy every boundary condition. These checks are separate from the residual, because the equation alone has a whole family of solutions, and a zero residual cannot tell which member was intended. For computed output, the same check confirms that the data actually entered are the data stated, including the starting time and the units.

In plain words Passing the equation shows that the curve belongs to the right family. Passing the data shows that it is the right member of that family. A computation needs both.

Example The function $v(t) = 1/(2 - t)$ has zero residual for $y' = y^2$ on $t < 2$ , since $v' = 1/(2 - t)^2 = v^2$ , but $v(0) = 1/2$ . It solves the equation and fails the problem with $y(0) = 1$ . In fact the zero function and every function $1/(c - t)$ solve the equation, and only the initial value picks out $c = 1$ . In the same way, take the boundary-value problem $y'' = -2$ on $0 \le x \le 1$ , with $x$ as the independent variable and $y(0) = y(1) = 0$ . The function $2x - x^2$ satisfies the equation and the left condition but has the value $1$ at $x = 1$ , so only the boundary check exposes it.

Event detection

Definition An event is a time $t^*$ at which a chosen scalar event function $g(t, y)$ of the solution is zero, such as the first time a state reaches a threshold. Event detection checks the sign of $g$ along the computed solution after each step. When $g$ changes sign between $t_n$ and $t_{n+1}$ , the solver locates the root inside the step from its interpolant, by bisection or a secant-type method, to a set tolerance. It records the direction of the crossing, increasing when $g$ goes from negative to positive, and it may stop the integration there. Reporting the first grid time past the threshold instead gives an event time that can be late by up to one step, and a sign test at grid points alone can miss two crossings that fall inside one step.

In plain words Often the answer wanted is a moment, not a whole curve. The solver watches a quantity, notices when it changes sign during a step, and then zooms in on the crossing instead of rounding to the next grid point.

Example For $y = 1/(1 - t)$ and the event function $g = y - 4$ , the equation $1/(1 - t) = 4$ gives $1 - t = 1/4$ , so $t^* = 3/4$ . There $y' = y^2 = 16 > 0$ , so $g$ crosses zero in the increasing direction. Suppose accurate values are known at $t = 0.7$ and $t = 0.8$ , namely $10/3$ and $5$ . The first grid time past the level is $0.8$ , which is $0.05$ late. Linear interpolation between the two points gives $t \approx 0.74$ , off by $0.01$ . A cubic interpolant that also uses the slopes $y^2$ at both points gives $t \approx 0.7504$ , off by about $4 \times 10^{-4}$ .

Events versus singularities

Definition An event lies inside the interval of existence, and the solution continues past it. A singularity, in the sense used here, is a finite endpoint of the maximal interval of existence, where the solution becomes unbounded or comes arbitrarily close to the edge of the region on which the equation is defined. No classical solution of the problem continues through it. A third kind of limit is a validity boundary of the model, beyond which the equation may still have a solution but no longer describes the system. A stopping event, a validity boundary and a singularity are different reasons to end a computation, and a report should say which one occurred.

In plain words Reaching a level is something that happens to a solution on its way. A singularity is the end of the road for that solution. A solver that keeps printing numbers past a singularity is not showing a continuation.

Example The solution $y = 1/(1 - t)$ has maximal interval $(-\infty, 1)$ , and both $y$ and $y' = y^2$ tend to $\infty$ as $t \to 1^-$ . The event at $t = 3/4$ is an ordinary point of this interval, while $t = 1$ is its end. Fixed-step methods do not notice the end. Euler with $h = 1/4$ produces the finite values $3.65$ at $t = 1$ and $6.99$ at $t = 1.25$ , so its polygon first passes the level $4$ after the singularity. None of the values at or beyond $t = 1$ describe the solution. The formula $1/(1 - t)$ is also defined for $t > 1$ , where it is negative, but that branch solves a different problem, such as $y(2) = -1$ , and no finite value at $t = 1$ joins it to the data at $t = 0$ . An adaptive solver usually reveals the singularity by shrinking its steps without limit near $t = 1$ and stopping with a warning, which is a signal to analyze the problem, not a failure to hide.

Definition A refinement study repeats a computation with smaller fixed steps, or with tighter tolerances for an adaptive solver, and compares the quantity of interest, such as the state at a fixed time or an event time. For a method of order $p$ on a smooth problem, the error at a fixed time behaves like $Ch^p$ for small $h$ , with $C$ independent of $h$ , so halving the step should divide the error, and the differences between successive runs, by about $2^p$ . The observed ratio approaches $2^p$ only once $h$ is small enough, and the difference between the two finest runs, divided by $2^p - 1$ , then estimates the error of the finer one. An independent check compares the result with evidence that could fail in a different way, such as an exact solution of a special case, a conserved quantity, a different method or solver, or a physical limit. Refinement shows that the digits of the computed problem have settled, not that the computed problem is the intended one.

In plain words Run the computation again with smaller steps and see whether the answer stops moving at the rate the method promises. Then test it against something that does not share its assumptions, because a wrong equation converges just as nicely as a right one.

Example Euler’s method for $y' = y^2$ with $y(0) = 1$ gives these values at $t = 1/2$ , where the exact value is $2$ .

Step $h$	Euler value	Error
$1/4$	$1.6406$	$0.3594$
$1/8$	$1.7662$	$0.2338$
$1/16$	$1.8611$	$0.1389$
$1/32$	$1.9230$	$0.0770$
$1/64$	$1.9593$	$0.0407$
$1/128$	$1.9790$	$0.0210$

The successive error ratios are $1.54$ , $1.68$ , $1.80$ , $1.89$ and $1.94$ . They approach $2$ , as they should for a first-order method, but only slowly, because the derivatives of the solution grow quickly toward the blowup at $t = 1$ and keep the higher-order error terms significant. With $p = 1$ , the difference $1.9790 - 1.9593 = 0.0197$ between the two finest runs is a fair estimate of the remaining error $0.0210$ . Here the exact solution serves as the independent check.

Unit consistency

Definition An equation is unit consistent when all terms that are added, subtracted or equated have the same dimension once their units are converted, and every argument of an exponential, a logarithm or a trigonometric function is dimensionless. A derivative $dC/dt$ carries the units of $C$ divided by those of $t$ . Parameters get their units from this requirement, and every numerical value must be entered in one consistent system of units. Error tolerances carry units too. A solver typically accepts a step when the estimated error in each component $y_i$ is at most $\text{atol}_i + \text{rtol}\,\lvert y_i \rvert$ , so an absolute tolerance has the units of its component and has to be chosen for each component separately.

In plain words Every term of a model has to measure the same kind of quantity in the same units, or the equation adds things that cannot be added. A unit check costs almost nothing and catches mistakes that no amount of refinement can find.

Example Suppose a concentration $C$ in $\mathrm{mol/L}$ obeys $dC/dt = kC^2$ with $t$ in minutes. The left side has units $\mathrm{mol/(L \cdot min)}$ and $C^2$ has units $\mathrm{mol^2/L^2}$ , so $k$ must have units $\mathrm{L/(mol \cdot min)}$ , which is inverse concentration per time. With $C(0) = C_0$ , the scaled variables $y = C/C_0$ and $\tau = kC_0t$ are dimensionless, and

\begin{gathered} \frac{dC}{dt} = kC_0^2\,\frac{dy}{d\tau}, \\[0.4em] kC^2 = kC_0^2\,y^2, \end{gathered}

so the model becomes $dy/d\tau = y^2$ with $y(0) = 1$ , the running example with $\tau$ in place of $t$ , as in the lesson on units and nondimensionalization. For $k = 0.25\ \mathrm{L/(mol \cdot min)}$ and $C_0 = 2\ \mathrm{mol/L}$ , the time scale is $1/(kC_0) = 2\ \mathrm{min}$ . The singularity at $\tau = 1$ is then at $t = 2\ \mathrm{min}$ , and the event $y = 4$ , which is $C = 8\ \mathrm{mol/L}$ , comes at $t = 1.5\ \mathrm{min}$ .

Check yourself

A program reports $u(t) = 2/(2 - t)$ as the solution of $y' = y^2$ with $y(0) = 1$ . Which checks does it pass, and which does it fail?

Show the answer

It passes the initial check, since $u(0) = 2/2 = 1$ . For the residual,

\begin{aligned} u' - u^2 &= \frac{2}{(2 - t)^2} - \frac{4}{(2 - t)^2} \\[0.4em] &= -\frac{2}{(2 - t)^2}, \end{aligned}

which is never zero and equals $-1/2$ at $t = 0$ , so it fails the equation. The function $1/(2 - t)$ fails the other way, with zero residual and the wrong initial value. Only $1/(1 - t)$ passes both checks.

For the running example, find the time when $y$ reaches $10$ and the direction of the crossing. A fixed-step Euler run with $h = 0.1$ reports $y(1.1) \approx 9.89$ and $y(1.2) \approx 19.7$ and places the crossing of $10$ between $t = 1.1$ and $t = 1.2$ . What is wrong with that report?

Show the answer

The equation $1/(1 - t) = 10$ gives $t = 0.9$ , and there $y' = y^2 = 100 > 0$ , so the crossing is increasing. The Euler report places the event after $t = 1$ , where the solution has already become unbounded. The maximal interval of the problem is $(-\infty, 1)$ , so the values at $t = 1.1$ and $t = 1.2$ belong to no solution of this problem, and the method has simply stepped over the singularity. Euler with this step is also far too coarse, since at $t = 0.9$ it gives only about $4.29$ instead of $10$ .

In the concentration model of the unit card, the rate constant is really $k = 0.25\ \mathrm{L/(mol \cdot s)}$ , but the program measures $t$ in minutes and uses the number $0.25$ . Runs with two very different tolerances agree that $C$ reaches $8\ \mathrm{mol/L}$ at $t = 1.5\ \mathrm{min}$ . What does the agreement show, and what is the correct event time?

Show the answer

The agreement shows only that the program solved its own problem accurately. That problem has the wrong rate constant, and refinement cannot detect it. In minutes the constant is $k = 0.25 \cdot 60 = 15\ \mathrm{L/(mol \cdot min)}$ , so $kC_0 = 30$ per minute and the time scale is $1/30\ \mathrm{min} = 2\ \mathrm{s}$ . The event $C = 8\ \mathrm{mol/L}$ is still at $\tau = 3/4$ , which is

\begin{aligned} t &= \frac{3/4}{30}\ \mathrm{min} \\[0.4em] &= 0.025\ \mathrm{min} = 1.5\ \mathrm{s}. \end{aligned}

The true solution blows up at $t = 2\ \mathrm{s}$ , so the reported time of $1.5\ \mathrm{min}$ lies far beyond its interval of existence. A unit check of $k$ catches the mistake at once.

Reliable Differential-Equation Computation

Key terms

Equation residual

Initial and boundary checks

Event detection

Events versus singularities

Convergence under refinement

Unit consistency

Check yourself

How this lesson connects

Builds on