Why Numerical Stability Matters in Dosing
A deployed dosing model can receive valid data and execute the intended equations yet still produce an unreliable result because of finite numerical precision, rounding, or an unsuitable solver tolerance. The practical concern is not whether two computers print identical digits. It is whether numerical variation can change a dose recommendation, its clinical interpretation, or the system’s ability to detect an invalid result.
Numerical stability concerns how computational errors affect the result produced by an algorithm. This is distinct from conditioning: conditioning describes how sensitive the mathematical problem is to small changes in its inputs, while stability describes how well the chosen computation controls the errors it introduces. A well-conditioned problem can be computed poorly, and a stable algorithm cannot remove sensitivity inherent in an ill-conditioned problem.
For a dosing application, numerical adequacy is judged in context. A difference of \(0.01\) mg/L may be immaterial when predictions are far from a decision threshold, but important if the recommendation changes at that threshold. Numerical checks should therefore relate computational discrepancies to units, model behavior, and the distance to a decision boundary.
Conditioning, Precision, and Small Differences
A simple local sensitivity measure for a differentiable scalar function \(f(x)\) is its relative condition number. It is defined only when both \(x\ne 0\) and \(f(x)\ne 0\), because relative input and output changes require division by those quantities.
For a small input change \(\Delta x\), differentiability gives \(f(x+\Delta x)-f(x)\approx f'(x)\Delta x\). Dividing by \(f(x)\) and rewriting \(\Delta x=(\Delta x/x)x\) gives the relative-change approximation
Thus, for small perturbations, the magnitude of the relative output change is approximately \(\kappa_f(x)\) times the magnitude of the relative input change. This is a local approximation, not a guarantee for large changes or near a discontinuity. If \(x=0\) or \(f(x)=0\), relative change is not an appropriate measure there; an absolute sensitivity measure, such as \(|f'(x)|\), may be more useful.
Worked Example: Input Sensitivity in an Exponential Prediction
Suppose a hypothetical one-compartment prediction is \(f(k)=\exp(-kt)\), where \(k=1.00\ \text{h}^{-1}\) and \(t=5.00\ \text{h}\). Assess the local effect of a 0.1% relative change in \(k\), holding time fixed.
Here \(f(k)=\exp(-kt)\) and \(f'(k)=-t\exp(-kt)\). Since \(k\ne0\) and \(f(k)\ne0\), the relative condition number is
A 0.1% relative perturbation in \(k\) has magnitude \(0.001\). The local approximation therefore predicts a relative change in \(f\) of about \(5.00(0.001)=0.00500\), or 0.500%. For a direct check, increase \(k\) by 0.1% to \(1.001\ \text{h}^{-1}\). The ratio of the new prediction to the old one is
The prediction decreases by about 0.499%, close to the local approximation. This quantifies sensitivity of the mathematical prediction; it does not by itself show that a software implementation is stable or that the input change is clinically plausible.
Finite precision can be especially troublesome when an algorithm subtracts nearly equal numbers. This is called cancellation: the subtraction may remove the leading, accurately represented digits and leave a small result dominated by rounding. For example, near 1, single-precision floating-point spacing is approximately \(1.19\times10^{-7}\). A value \(1+4\times10^{-8}\) can round to the same representable number as 1. Subtracting the stored values then gives zero, although the exact difference is \(4\times10^{-8}\). Whether that matters depends on how the difference is used; a near-zero residual or decision margin can be more consequential than a large predicted concentration.
Worked Example: A Rounded Margin Changes a Branch
Consider a hypothetical dosing algorithm that increases a dose if a predicted trough is below \(5.000\ \text{mg/L}\), and otherwise retains the dose. A full-precision calculation gives \(4.9996\ \text{mg/L}\). Show what happens if the value is rounded to three decimal places before the comparison.
Using the unrounded value, \(4.9996<5.000\), so the algorithm selects the increase branch. If the value is first rounded to three decimal places, it becomes \(5.000\ \text{mg/L}\). The comparison is now \(5.000<5.000\), which is false, so the algorithm selects the retain branch.
The numerical discrepancy is only \(0.0004\ \text{mg/L}\), but it changes the branch because the prediction is close to the cutoff. This example does not establish which clinical action is appropriate; it shows why internal rounding and boundary behavior must be checked against the actual decision rule. The displayed precision may be lower than the internal precision, but the internal value and comparison rule must be specified consistently.
Solver Tolerances and Model Boundaries
Many dosing models solve differential equations numerically. An adaptive solver estimates local error and adjusts step size to meet specified absolute and relative tolerances. A common error scaling for a state \(y_i\) is based on \(\text{atol}_i+\text{rtol}|y_i|\), where \(\text{atol}_i\) has the units of \(y_i\) and \(\text{rtol}\) is dimensionless. This scaling helps accommodate states of different magnitudes. It is not, by itself, a guarantee that the final prediction error is below that amount: local errors accumulate and the solver’s error estimate is not necessarily the true error.
Tolerances should be assessed in the units and scale of each modeled state. A single absolute tolerance may be unsuitable when one state is near zero and another is large, or when states have different units. In pharmacokinetic models, dosing events can also create discontinuities in the input or state. The solver must handle these events at the specified times, rather than stepping across them as if the model were smooth. Output should be evaluated at the intended event-relative time, with a clearly defined convention for a sample exactly at a dose time.
Worked Example: Interpreting a Solver Error Scale
For a hypothetical one-compartment model, \(dC/dt=-kC\), with \(C(0)=10.0\ \text{mg/L}\) and \(k=0.200\ \text{h}^{-1}\), the analytical concentration at \(t=5.00\ \text{h}\) is
Suppose a solver uses \(\text{atol}=0.001\ \text{mg/L}\) and \(\text{rtol}=0.001\). At a state near \(3.679\ \text{mg/L}\), the error scale is approximately
If the solver’s estimated local error for a proposed step is \(0.0040\ \text{mg/L}\), its normalized estimate is \(0.0040/0.004679\approx0.855\), below 1; under this common criterion the step can be accepted. If the estimate is \(0.0060\ \text{mg/L}\), the ratio is \(0.0060/0.004679\approx1.28\), so the step would be rejected and typically retried with a smaller step. These calculations explain the local acceptance scale, not a guaranteed error in the final concentration. To assess final-output adequacy, tighten tolerances and compare predictions with an independent calculation or a suitable analytical solution.
Checks for a Dosing Application
The checks below complement the reference test cases covered in the earlier material on computational verification. They focus on whether a result remains numerically credible when precision, solver settings, or inputs near a boundary are varied.
Set allowable differences in clinically meaningful units, such as mg/L for concentration or mg for dose. Relate the limits to the decision rule and its margins; do not use a universal tolerance without regard to the application.
Compare the production computation with a higher-precision calculation where feasible, and tighten solver tolerances in stages. Record changes in both continuous outputs and discrete recommendation branches.
Test values just below, at, and just above thresholds; valid input limits; zero or near-zero states; and dose or sample times. Confirm that branch rules, units, and event-time conventions are explicit.
Where the model requires it, verify nonnegative concentrations, finite outputs, expected monotonicity, and conservation or balance relationships. Treat a violated invariant as a diagnostic signal rather than silently clipping the result.
Compare the numerical difference with the distance to the nearest decision boundary. If plausible numerical variation can change a recommendation, flag or recompute the result according to a prespecified rule rather than implying false precision.
Tighter tolerances are a useful diagnostic, not a complete validation strategy. If predictions change materially as tolerances tighten, investigate solver configuration, discontinuities, scaling, or model implementation. Agreement between two calculations is stronger evidence when they use meaningfully different numerical methods; agreement alone is not proof that both correctly represent the intended model.
Approaches, Failure Modes, and Limitations
| Approach | Strength | Limitation or failure mode |
|---|---|---|
| More internal precision | Can reduce representation and intermediate-rounding error. | Does not fix an ill-conditioned model, incorrect equations, or inappropriate event handling. |
| Tighter solver tolerances | Can reveal whether an ODE result is sensitive to numerical integration settings. | May increase computation and still does not guarantee a clinically adequate final error bound. |
| Boundary and invariant tests | Expose discontinuous decisions, invalid values, and violations of expected behavior. | Require model-specific expectations and do not cover every possible input. |
| Independent numerical method or analytical comparison | Can help distinguish implementation or solver problems from ordinary rounding differences. | Methods can share assumptions or errors; comparison is not a substitute for checking the model specification. |
A common mistake is to assume that more decimal places mean more accuracy. Printing ten digits cannot recover information lost by rounding inputs or intermediate values, and it cannot correct a biased or misspecified model. The correction is to preserve adequate internal precision, avoid unnecessary rounding before a decision, test the sensitivity of outputs to computational settings, and report only precision justified by the model and data.
Another failure is treating a solver’s tolerance as a clinical error bound. Tolerances govern an algorithm’s local error-control procedure; they are not automatically a bound on the final prediction or on dose error. Numerical checks also cannot resolve uncertainty in the patient’s measurements, model parameters, or structural assumptions. Those uncertainties require their own assessment.
Check Your Understanding
Answer each question, then compare your reasoning with the brief solutions.
- A calculation reports a value of \(0.0000\) after subtracting two nearly equal numbers. Why might this be misleading, and what check would you perform?
- For \(f(x)=\exp(-2x)\) at \(x=3\), calculate the relative condition number and estimate the relative output change from a 0.2% relative input change.
- A solver reports a local error estimate of \(0.003\ \text{mg/L}\) at a state of \(2.0\ \text{mg/L}\), with \(\text{atol}=0.001\ \text{mg/L}\) and \(\text{rtol}=0.001\). Calculate the normalized error estimate using the stated scale and say whether it is below 1.
- A prediction is \(6.004\ \text{mg/L}\), and the algorithm changes branch at \(6.000\ \text{mg/L}\). Why should the numerical difference be assessed in relation to the boundary rather than only by its number of decimal places?
Solutions. (1) Cancellation and finite precision may have erased a small nonzero difference. Recompute at higher precision or reformulate the calculation, then compare with an independent method if possible. (2) Since \(f'(x)=-2\exp(-2x)\), \(\kappa_f(x)=|x f'(x)/f(x)|=2x=6\). The approximate relative output change has magnitude \(6(0.002)=0.012\), or 1.2%. (3) The scale is \(0.001+0.001(2.0)=0.003\ \text{mg/L}\); the normalized estimate is \(0.003/0.003=1\), not below 1. (4) The prediction is only \(0.004\ \text{mg/L}\) above the branch point, so numerical variation of that order could change the recommendation. Check stability under precision and tolerance changes and define what the application does when the result is too close to the boundary.