8
$\begingroup$

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.

Improve this question
$\endgroup$
2
  • $\begingroup$ Most common methods would be unsuitable because they have errors that behave as $h^{s}$, so requiring $h^s\sim 10^{-M}$ leads to number of steps that is exponential in $M$: $1/h\sim 10^{M/s}$. Presumably some methods like spectral methods would allow you to do this if you can show the error is exponentially small in $1/h$. Since each step is a finite number of rational operations on numbers of length $O(M)$, not counting black-box evaluations, each taking time polynomial in $M$, presumably it is the step size that matters the most. $\endgroup$
    – Kirill
    May 22 '15 at 5:51
  • $\begingroup$ @Kirill: In a spectral approach (or using other traditional ODE methods), I think the OP would need to know leading constants for the asymptotics to get a certificate of accuracy. These constants would come from analysis, or computation via interval arithmetic. $\endgroup$
    – Geoff Oxberry
    May 22 '15 at 6:15
9
$\begingroup$

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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ I am tempted to retract my answer because @geoffoxberry's is just so much better than mine... $\endgroup$
    – Wolfgang Bangerth
    May 22 '15 at 6:38
2
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.