Minimizing Residuals in ODE Integration Using Optimal Control

TL;DR

This work treats the residual (defect) of an ODE solution as a controllable perturbation and seeks interpolants through a given skeleton that minimize this residual. By formulating the interpolation as a multi-stage optimal-control problem, it derives necessary conditions of optimality for both $L^2$ and stage $L^{\infty}$ objectives, and analyzes normality to yield tractable two-point boundary-value problems. The authors provide analytic solutions for the Dahlquist test problem and a leaky-bucket variant, and perform numerical optimization for the Van der Pol equation, comparing the resulting residuals with MATLAB's deval across several solvers. The results show that the proposed residual-minimizing interpolants often yield substantially smaller residuals than standard polynomial interpolants, offering a robust diagnostic tool for validating numerical ODE solutions and understanding solver behavior, including order-reduction phenomena. This framework can enhance reproducibility and reliability in numerical ODE integration by quantifying how closely a skeleton fits a perturbed model via minimal residuals.

Abstract

Given the set of discrete solution points or nodes, called the skeleton, generated by an ODE solver, we study the problem of fitting a curve passing through the nodes in the skeleton minimizing a norm of the residual vector of the ODE. We reformulate this interpolation problem as a multi-stage optimal control problem and, for the minimization of two different norms, we apply the associated maximum principle to obtain the necessary conditions of optimality. We solve the problem analytically for the Dahlquist test problem and a variant of the leaky bucket problem, in terms of the given skeleton. We also consider the Van der Pol equation, for which we obtain interpolating curves with minimal residual norms by numerically solving a direct discretization of the problem through optimization software. With the skeletons obtained by various ODE solvers of MATLAB, we make comparisons between the residuals obtained by our approach and those obtained by the MATLAB function deval.

Figures (9)

• Figure 1: A kind of work-precision diagram for the solution of the simple harmonic oscillator $\ddot y + y = 0$ on $0 \le t \le 30$ with reasonable initial conditions, given to ode45 in the form of a first-order system as usual. On the horizontal axis we have the mean stepsize $h$. The code ode45 is a variable-stepsize method, and the stepsize at the tightest tolerance varied systematically between $9\times 10^{-3}$ and $1.2\times 10^{-2}$. Smaller mean $h$ implies that the code had to work harder to get to $t=30$. On the vertical axis, we have the maximum residual norm observed over the solution on $0 \le t \le 30$. We see that the expected order of the error, namely $O(h^5)$, plotted in red, is not attained; instead, the observed order of the residual error is very clearly $O(h^4)$ (plotted in blue).
• Figure 2: Example 1 -- (a) interpolant for $\dot{z}(t) = 3\,z(t)$, $0\le t\le1$, $z(0)=1$, and interpolant residuals computed using (b),(c),(d) Matlab's deval, (e),(f),(g) the $L^2$-residual in \ref{['simp_res_L2']} and (h),(i),(j) the stage $L^\infty$-residual in \ref{['simp_res_Linfty']}. The skeletons in each row were obtained using the Matlab packages ode15s, ode45 and ode113, as indicated on the graphs. We see incidentally that the polynomial interpolant for ode45 is indeed continuously differentiable, while those for ode15s and ode113 are not, although the residual (in (d)) for ode113 is zero as $s \to 1^-$ on each step.
• Figure 3: Example 2 -- (a) interpolant for $\dot{z}(t) = \sqrt{z(t)}$, $0\le t\le1$, $z(0)=1$, and interpolant residuals computed using (b),(c),(d) Matlab's deval, (e),(f),(g) the $L^2$-residual in \ref{['simp_res_L2']} and (h),(i),(j) the stage $L^\infty$-residual in \ref{['simp_res_Linfty']}. The skeletons in each row were obtained using the Matlab packages ode15s, ode45 and ode113, as indicated on the graphs. Figures (d), (g), and (j) are not very informative, because ode113 actually gets the exact solution to this problem: the minimal residual is zero, apart from roundoff error.
• Figure 4: Example 3 -- interpolant residual vector components computed for equation \ref{['ode:vdp']} and Matlab's ode15s, with $0\le t\le1$ and $x(0)=(-1,-3)$, using (a)--(b) Matlab's deval, (c)--(d) $L^2$-minimization and (e)--(f) stage $L^\infty$-minimization.
• Figure 5: Example 3 -- (a)--(b) control variable vector and (c)--(d) costate variable vector components, all for equation \ref{['ode:vdp']}, with $0\le t\le1$, $x(0)=(-1,-3)$, and Matlab's ode15s.

Theorems & Definitions (17)

• Lemma 1: Normality
• Remark 1: Optimal control
• Theorem 1: $L^2$-residual
• Remark 2
• Lemma 2: Normality
• Theorem 2: Stage $L^\infty$-residual
• Remark 3
• Lemma 3: Partial singularity
• Proposition 1: $L^2$-residual of \ref{['simpleeqn']}
• Remark 4