Testing ODE Functions Against Analytic Solutions

Why Compare an ODE Implementation With an Analytic Solution?

A numerical solver can return a smooth, plausible trajectory even when the function supplied to it represents the wrong model. A decay term with the wrong sign, a missing input, or an incorrect parameter can all produce numerical output without an integration failure. When a model has a known analytic solution, that solution supplies an independent benchmark: the implemented ODE and its numerical solution should agree with the exact trajectory within a defined tolerance.

This is a focused check of both the right-hand-side function and its use in numerical integration. It does not establish that a model is biologically appropriate, nor does it replace checks of units, state indexing, or conservation. Those checks help identify errors; an analytic comparison tests whether a simple implemented model produces its known trajectory.

Two Analytic Benchmarks

For exponential decay, let \(x(t)\) be an amount or concentration, \(k>0\) a first-order loss-rate constant, and \(x_0\) the value at time \(t=0\). The model and its solution are:

\frac{dx}{dt}=-kx,\qquad x(0)=x_0, \qquad x(t)=x_0 e^{-kt}.

The solution applies when the state changes only through this first-order loss, with constant \(k\), and there are no doses or other time-varying inputs over the interval being tested. If \(x\) has units of amount per volume, \(k\) has units of inverse time, so \(kt\) is dimensionless and the solution has the same units as \(x_0\).

For constant-input turnover, let \(R(t)\) be a state, \(k_{\mathrm{in}}\) a constant input rate, and \(k_{\mathrm{out}}>0\) a first-order loss-rate constant. With initial value \(R_0\), the model and solution are:

\frac{dR}{dt}=k_{\mathrm{in}}-k_{\mathrm{out}}R,\qquad R(0)=R_0, \qquad R(t)=\frac{k_{\mathrm{in}}}{k_{\mathrm{out}}} +\left(R_0-\frac{k_{\mathrm{in}}}{k_{\mathrm{out}}}\right)e^{-k_{\mathrm{out}}t}.

Here, \(k_{\mathrm{in}}\) has units of state per time and \(k_{\mathrm{out}}\) has units of inverse time. The benchmark assumes both rates are constant over the interval, the initial condition is \(R_0\), and no additional processes or discontinuous inputs act on \(R\). The analytic expression gives the exact solution of this specified mathematical model, not a claim that the model is a complete description of a biological system.

Definition: An analytic benchmark test compares a numerical solution of an implemented ODE with a known exact solution of the same ODE, evaluated at the same times, initial conditions, and parameter values. The comparison is meaningful only when the implemented equations and the analytic benchmark represent the same model.

How the Comparison Works

For each chosen time \(t_i\), compute the analytic value \(x(t_i)\) or \(R(t_i)\), obtain the numerical value from the solver at that same time, and calculate their discrepancy. A direct measure is absolute error, which retains the units of the state. For values that vary greatly in scale, a combined absolute-relative rule is often more informative:

\left|y_{\mathrm{num}}(t_i)-y_{\mathrm{exact}}(t_i)\right| \leq \mathrm{atol}+\mathrm{rtol}\left|y_{\mathrm{exact}}(t_i)\right|.

Here, \(\mathrm{atol}\) is an absolute tolerance in the state’s units and \(\mathrm{rtol}\) is a dimensionless relative tolerance. The values should be chosen before interpreting the comparison, with regard to the state scale and the precision required for the intended calculation. This acceptance rule is a test criterion; it does not imply that the solver’s internal error estimate guarantees the same bound at every requested output time.

Use multiple times, including the initial time and intermediate points, rather than checking only the final output. The initial-time comparison checks that the implemented simulation starts from the stated initial condition. Intermediate comparisons can reveal a wrong rate or trajectory shape that an endpoint-only check misses. In a system with several states, compare each state separately in its own units.

Deriving the Turnover Solution

The turnover expression follows by solving the linear ODE with an integrating factor. Rewrite the equation as \(dR/dt+k_{\mathrm{out}}R=k_{\mathrm{in}}\), then multiply both sides by \(e^{k_{\mathrm{out}}t}\). The left side becomes the derivative of a product:

e^{k_{\mathrm{out}}t}\frac{dR}{dt} +k_{\mathrm{out}}e^{k_{\mathrm{out}}t}R = \frac{d}{dt}\left(e^{k_{\mathrm{out}}t}R(t)\right) = k_{\mathrm{in}}e^{k_{\mathrm{out}}t}.

Integrate from \(0\) to \(t\) and apply \(R(0)=R_0\):

e^{k_{\mathrm{out}}t}R(t)-R_0 = \frac{k_{\mathrm{in}}}{k_{\mathrm{out}}} \left(e^{k_{\mathrm{out}}t}-1\right).

Dividing by \(e^{k_{\mathrm{out}}t}\) and rearranging gives \(R(t)=k_{\mathrm{in}}/k_{\mathrm{out}}+(R_0-k_{\mathrm{in}}/k_{\mathrm{out}})e^{-k_{\mathrm{out}}t}\), as stated. This derivation also clarifies why the input must be constant for this particular formula: the integral on the right uses a constant \(k_{\mathrm{in}}\).

Worked Tests

Worked Example: Testing Exponential Decay at Several Times

Suppose a hypothetical state \(x\), measured in mg/L, is implemented as \(dx/dt=-kx\), with \(x_0=12\) mg/L and \(k=0.35\ \mathrm{h}^{-1}\). The solver is asked for values at \(t=0\), \(1\), and \(4\) h. The constant-rate, no-input assumptions match the analytic decay model, and the initial value and rate are specified in consistent units.

The reference values are \(x(t)=12e^{-0.35t}\) mg/L. At \(t=1\) h:

x(1)=12e^{-0.35(1)} =12(0.7046881) =8.456257\ \mathrm{mg/L}.

At \(t=4\) h:

x(4)=12e^{-0.35(4)} =12(0.2465970) =2.959164\ \mathrm{mg/L}.

Suppose hypothetical solver output is 12.0000, 8.4560, and 2.9600 mg/L at those times. The initial value agrees exactly. The absolute discrepancies at 1 and 4 h are \(|8.4560-8.456257|=0.000257\) mg/L and \(|2.9600-2.959164|=0.000836\) mg/L, respectively. These values would pass a preselected absolute tolerance of 0.001 mg/L at both times. That result supports agreement for these tested conditions; it does not show that all parameter values or other equations in a larger model are correct.

Worked Example: Benchmarking Constant-Input Turnover

Consider a hypothetical concentration \(R\) in nM, with constant input \(k_{\mathrm{in}}=6\ \mathrm{nM/h}\), loss rate \(k_{\mathrm{out}}=0.2\ \mathrm{h}^{-1}\), and initial value \(R_0=20\) nM. The inputs are constant over the test interval, and the implemented ODE is \(dR/dt=6-0.2R\). The analytic benchmark has \(k_{\mathrm{in}}/k_{\mathrm{out}}=30\) nM:

R(t)=30+(20-30)e^{-0.2t} =30-10e^{-0.2t}\ \mathrm{nM}.

At \(t=0\), the reference is \(R(0)=30-10=20\) nM. At \(t=3\) h:

R(3)=30-10e^{-0.6} =30-10(0.5488116361) =24.511883639\ \mathrm{nM}.

At \(t=10\) h, \(R(10)=30-10e^{-2}=28.646647168\) nM. Suppose the numerical output at 3 h is 24.5200 nM. Its absolute discrepancy is \(|24.5200-24.511883639|=0.008116361\) nM, or 0.0081 nM to four decimal places. If the stated absolute tolerance is 0.01 nM, this result passes at 3 h. The comparison at 10 h must also be made before concluding that the tested trajectory passes at all selected times.

Worked Example: Detecting a Missing Input Term

Use the turnover model and parameters from the preceding example, but suppose a coding error omits the input and the implemented function returns \(dR/dt=-0.2R\). A solver may still integrate this incorrect function without warning. The correct analytic benchmark at 3 h remains \(24.511883639\) nM, because the intended model still includes the 6 nM/h input.

The erroneous function instead generates pure decay from \(R_0=20\) nM:

R_{\mathrm{bug}}(3)=20e^{-0.2(3)} =20(0.5488116361) =10.976232722\ \mathrm{nM}.

The discrepancy from the intended solution is \(|10.976232722-24.511883639|=13.535650917\) nM, far larger than a reasonable numerical tolerance for this example. Since the formula, units, initial value, and constant-rate assumptions match the intended benchmark, the large difference points to an equation or implementation mismatch rather than ordinary rounding. A separate evaluation of the right-hand side at \(R=20\) nM would also help localize the issue: the intended derivative is \(6-0.2(20)=2\ \mathrm{nM/h}\), whereas the buggy function returns \(-4\ \mathrm{nM/h}\).

Interpreting Discrepancies and Common Failure Modes

A passing comparison means the numerical outputs agree with the analytic solution to the chosen precision at the tested conditions and times. A failure does not, by itself, identify whether the cause is the ODE function, parameter values, initial conditions, units, or numerical integration. Inspect the implementation and the benchmark inputs together. In particular, confirm that time is expressed in the units used by the rate constants, the initial state is passed as intended, and the analytic formula corresponds to the actual input schedule.

The main alternatives answer different questions. Directly evaluating the right-hand-side function at a known state checks whether a derivative is calculated as expected, but does not test integration over time. Comparing against an analytic trajectory tests the combined function-and-solver result for a model with a known solution. For models without such a solution, comparisons across solver settings or with an independently implemented numerical method can provide useful evidence, but they are not exact analytic benchmarks.

Checking only the endpoint. An endpoint can agree even when intermediate values do not. Compare at the initial time and at multiple points along the trajectory.
Comparing different models. An analytic formula for constant input is not a valid reference for an implementation with a changing input. Match the model assumptions before judging errors.
Using inconsistent units. A rate in \(\mathrm{h}^{-1}\) must be paired with time in hours. A unit mismatch changes the trajectory even if the code runs.
Treating a pass as proof. A test covers only its tested equations, parameter values, initial conditions, and times. It cannot guarantee that a complex model is correct everywhere.

If discrepancies are unexpectedly large, first verify the analytic calculation and inputs, then check the implemented equation and state mapping, and finally inspect numerical integration settings and warnings. If discrepancies are small but exceed the prespecified test criterion, report that failure rather than silently loosening the tolerance after seeing the result. This makes the test interpretable and repeatable.

Specify the benchmark model.
Write its ODE, parameters, units, initial condition, and assumptions explicitly.

Calculate exact values.
Evaluate the analytic solution at the same times requested from the solver.

Compare against a stated rule.
Calculate discrepancies in state units and apply the chosen tolerance consistently.

Investigate failures.
Check inputs, units, initial values, right-hand-side calculations, and numerical output before concluding why the test failed.

Key takeaway: Analytic solutions for exponential decay and constant-input turnover provide exact, interpretable references for testing ODE implementations. Their value depends on matching the equation, assumptions, units, parameters, and initial conditions—not merely on obtaining a solver output.

Check Your Understanding

Use the stated model and show how you would compare a numerical result with its analytic benchmark.

For \(dx/dt=-0.25x\), \(x(0)=8\) mg/L, calculate the exact value at 2 h.
For \(dR/dt=4-0.5R\), \(R(0)=2\) units, calculate the exact value at 2 h.
A decay simulation matches its analytic solution at \(t=0\) but differs greatly at an intermediate time. Does the initial match establish that the implementation is correct? Explain.
For \(x(t)=10e^{-0.2t}\), suppose the numerical value at 3 h is 5.490. Calculate the absolute discrepancy and determine whether it passes \(\mathrm{atol}=0.002\) and \(\mathrm{rtol}=0.001\).

Solutions. (1) \(x(2)=8e^{-0.25(2)}=8e^{-0.5}=4.852245\) mg/L. (2) The constant-input solution is \(R(t)=4/0.5+(2-4/0.5)e^{-0.5t}=8-6e^{-0.5t}\), so \(R(2)=8-6e^{-1}=5.792723\) units. (3) No. Matching at \(t=0\) checks the initial value, but does not establish that the rate equation or later trajectory is correct; compare intermediate and later times as well. (4) The exact value is \(10e^{-0.6}=5.488116\), so the absolute discrepancy is \(|5.490-5.488116|=0.001884\). The allowed difference is \(0.002+0.001(5.488116)=0.007488\), so the result passes the stated rule.