How can I numerically solve an ODE to $N$ provably correct digits?

Question

Suppose we have an initial value problem of the form $$ \frac{\mathrm{d} \mathbf{x}}{\mathrm{d} t} = f(\mathbf{x}) \qquad \mathbf{x}(0) = \mathbf{x}_0 $$ where $\mathbf{x}_0 \in \mathbb{R}^n$ is known exactly (i.e. to unlimited precision) and we can efficiently evaluate $f: \mathbb{R}^n \to \mathbb{R}^n$ to any precision. That is to say, we have a black box which, given a vector $\mathbf{x} \in \mathbb{R}^n$ and an integer $M$, returns an approximation to $f(\mathbf{x})$ guaranteed to be correct to $M$ digits in time polynomial in $M$. I would like to know if there are any practical methods for obtaining an approximation to $\mathbf{x}(t_f)$ (where $t_f \in \mathbb{R}$ is some given final time) which is provably correct to $N$ digits.

Clearly, this will not be possible for just any function $f: \mathbb{R}^n \to \mathbb{R}^n$, since $f$ might have some crazy behavior that drastically alters the true solution but isn't picked up in any reasonable number of evaluations. Thus, I am also interested in knowing what kind of well-behavior conditions on $f$ (e.g. all partial derivatives exist and are bounded, small Lipschitz constant, etc.) would be necessary to do this.

Answer 1

I would like to know if there are any practical methods for obtaining an approximation to $\mathbf{x}(t_f)$ (where $t_f \in \mathbb{R}$ is some given final time) which is provably correct to $N$ [sic] digits.

That all depends on your opinion of the practicality of interval arithmetic. There are validated integrators available, such as the COSY code out of Martin Berz's group. You'd probably want to look at papers by Neumaier, Nedialkov, Berz & Makino, Chachuat, Stadtherr, and maybe a few other groups. Their papers tend to use the phrases "Taylor model", "validated integrator", and "interval arithmetic" among others.

Clearly, this will not be possible for just any function $f: \mathbb{R}^{n} \rightarrow \mathbb{R}^{n}$, since $f$ might have some crazy behavior that drastically alters the true solution but isn't picked up in any reasonable number of evaluations. Thus, I am also interested in knowing what kind of well-behavior conditions on $f$ (e.g. all partial derivatives exist and are bounded, small Lipschitz constant, etc.) would be necessary to do this.

The standard proofs (here, I am thinking of Dahlquist's proof, for instance) typically use Gronwall-type inequalities for error bounds, so in theory, you should just need a bounded Lipschitz constant over the domain in $\mathbb{R}^{n}$ you care about, which speaks to some of what WolfgangBangerth is talking about, although I don't know that a graduate level textbook like Hairer and Wanner would discuss specifically about the accuracy which which you can evaluate $f(\mathbf{x})$; Higham's book on accuracy and stability of numerical methods might discuss this point.

In practice, for the Taylor model methods I mentioned above, you typically need a combination of theoretical and practical conditions. The theory part is straightforward: if you want a $k$th order Taylor model, you need a function in $C^{k}$. The implementation part, which speaks to your issue about practicality, is harder.

From a user perspective, these integrators boil down to a two(-ish) things:

• Do I have the source code to evaluate $f$? (An object file isn't enough.)
• Can I augment this source code for compatibility with an automatic differentiation library, as well as an interval arithmetic library?

Assuming the function can be parsed or operator-overloaded (depending on the automatic differentiation technology used), and assuming you can generate a function that will calculate the interval extension of $f$ and its first $k$ derivatives, then you can implement a validated integration method with $k$th-order Taylor models.

As for typical ODE solving methods, to comment on Wolfgang's answer:

I don't think that you can get a certificate that the error is below a certain number, but you will get that the estimate is below your tolerance.

Any method that has an embedded error estimator has the information Wolfgang refers to. Usually, this means that the integration method really calculates two (or more solutions; e.g., DOP853 calculates 3 solutions) solutions and compares them via some norm. The assumption is that the higher order solution is more accurate, which may not actually be true, depending on the given problem, time step, initial conditions, etc. The solution returned by an implementation could be any of the candidate solutions computed. Taking the common Runge-Kutta 4(5) case as an example, one could return the 4th order solution or the 5th order solution; typical approaches use the Dormand-Prince formulas, which minimize the error in the 5th order solution, and return that, rather than the 4th order solution, because the 5th order solution is more likely to be more accurate. In addition to looking at stability issues, I think you should look at error control (Section II.4 of Hairer and Wanner); stability is necessary, but not sufficient for accuracy.

In contrast, a validated integrator will compute an interval containing the true solution. If the upper and lower bound on this interval agree to $M$ digits, then those are the first $M$ digits in the true solution, which is the certificate of accuracy, and sounds like what you want.

Answer 2

The first of your questions is something you can actually get from most of the canned ODE integrators since they all, in one way or the other, keep track of estimates of the error. I don't think that you can get a certificate that the error is below a certain number, but you will get that the estimate is below your tolerance.

The second of your questions is harder to answer: What is the relationship between the accuracy with which you can evaluate $f(\mathbf x)$ and how that affects the accuracy of your solution (assuming that you could integrate the ODE exactly). The ratio of these two accuracies is just the stability constant of your ODE. It will depend on things like the norm of $f$, the norm of $\nabla f$, and the length of the time period over which you integrate (and it will typically depend exponentially on this interval length). Most textbooks on numerical ODEs will have a section on stability.