Optimal ODE method for fixed number of RHS evaluations

Question

In practice, the runtime of numerically solving an IVP $$ \dot{x}(t) = f(t, x(t)) \quad \text{ for } t \in [t_0, t_1] $$ $$ x(t_0) = x_0 $$ is often dominated by the duration of evaluating the right-hand side (RHS) $f$. Let us therefore assume that all other operations are instant (i.e. without computational cost). If the overall runtime for solving the IVP is limited then this is equivalent to limiting the number of evaluations of $f$ to some $N \in \mathbb{N}$.

We are only interested in the final value $x(t_1)$.

I'm looking for theoretical and practical results that help me choose the best ODE method in such a setting.

If, for example, $N = 2$ then we could solve the IVP using two explicit Euler steps of width $(t_1 - t_0)/2$ or one step of width $t_1 - t_0$ using the midpoint method. It is not immediately clear to me which one is preferrable. For larger $N$, one can of course also think about multi-step methods, iterated Runge-Kutta schemes, etc.

What I'm looking for are results similar to the ones that exist, for example, for quadrature rules: We can pick $n$ weights $\{w_i\}$ and associated points $\{x_i\}$ such that the quadrature rule $\sum_{i = 1}^n w_i g(x_i)$ is exact for all polynomials $g$ such that $deg(g) \le 2n - 1$.

Hence I'm looking for upper or lower bounds on the global accuracy of ODE methods, given a limited number of allowed evaluations of the RHS $f$. It's OK if the bounds only hold for some classes of RHS or pose additional constraints on the solution $x$ (just like the result for the quadrature rule which only holds for polynomials up to a certain degree).

EDIT: Some background information: This is for hard real-time applications, i.e. the result $x(t_1)$ must be available before a known deadline. Hence the limit on the number of RHS evaluations $N$ as the dominating cost factor. Typically our problems are stiff and comparatively small.

EDIT2: Unfortunately I don't have the precise timing requirements, but it is safe to assume that $N$ will be rather small (definitely <100, propably closer to 10). Given the real-time requirement we have to find a tradeoff between the accuracy of the models (with better models leading to longer execution times of the RHS and hence to a lower $N$) and the accuracy of the ODE method (with better methods requiring higher values of $N$).

Answer 1

I think a key reference to answer your question is this paper by Hosea and Shampine. Now I will give some background.

In general, the step size that you may use when numerically integrating an IVP can be restricted by stability or accuracy. If you want to choose the best solver in terms of stability, you need to consider the region of absolute stability. For a one-step method this is

$$S = \{ z \in {\mathbb C} : |P(z)|\le 1\}.$$

Here $P(z)$ is the stability function of the method; see e.g. the text of Hairer et. al. A necessary condition for stability is that $\lambda h \in S$ where $\lambda$ ranges over the eigenvalues of the jacobian of $f$ and $h$ is the step size. This is not always a sufficient condition for nonlinear problems but it is usually a good rule of thumb and is used in practice.

For an extensive treatment of the problem of finding (explicit) methods that allow large stable step sizes, see this paper of mine on stability polynomials and this one on optimization of Runge-Kutta methods for compressible fluid simulations.

Stability is the relevant concern if you find that the largest stable step size already gives you enough accuracy. On the other hand, the step size may instead be restricted by your accuracy requirements. What is typically done is local error control. The solution is computed using two methods, and their difference is used as an estimate of the error in the less accurate one. The step size is chosen adaptively to achieve as closely as possible the prescribed tolerance.

Two theoretical measures are key to predicting accuracy efficiency. The first is the order of accuracy of the method, which describes the rate at which the error approaches zero when the step size is decreased. The second is the accuracy efficiency index (see the paper of Hosea and Shampine linked in the first sentence above) which takes into account the constants appearing in the error terms and allows comparisons between methods of the same order.

Accuracy and stability efficiency of a wide range of methods can be computed in a simple and automated way using NodePy (disclaimer: NodePy is developed by me).

Answer 2

There aren't many results in this direction because it is more difficult than just fixing accuracy, since stability considerations can often require you to pick time-steps that are smaller than you would need for the accuracy you want. So results are split between the stiff and non-stiff cases. The former case the time-steps and RHS evaluations requirements are generally not governed by accuracy, and in the latter case they are.

I'm going to focus on explicit methods, since the implicit case is far less obvious how many RHS evaluations you will need to use.. that entirely depends on how you decide to solve the resulting system.

For non-stiff systems:

There are stage limitations for explicit Runge-Kutta methods for which say how many stages (RHS evaluations) are needed to achieve a certain order of accuracy. After fourth order the number of stages exceeds the order of accuracy, and the disparity continues to grow. Butcher's big ODE book: http://books.google.com/books/about/Numerical_Methods_for_Ordinary_Different.html?id=opd2NkBmMxsC

does a good job explaining some of these 'non-existence' proofs.

Your quadrature rule example leads either to a multistep type method such as Adams-bashforth, or to what are now called spectral-deferred-correction methods. For adams-bashforth you only need one RHS evaluation per step, but since the stability regions are so small in general for these methods you generally end up doing the same amount of work in terms of RHS evaluations as a Runge-Kutta method with the same order.

Here is a paper on spectral deferred correction:

https://www.google.com/search?q=spectral+deferred+correction&aq=f&oq=spectral+deferred+correction&aqs=chrome.0.57j0l2j62.3336j0&sourceid=chrome&ie=UTF-8

I am unsure how these integration methods play out against standard explicit methods, they often require a lot more memory to save solution states at quadrature nodes and so I have never used them myself.

For stiff systems:

there are 'optimized' time-steppers, but precise theoretical results concerning how good these can get are unfortunately limited to some simple cases (and even those turned out to not be a trivial bit of work). The three standard results say that for Runge-Kutta methods with $S$ stages: The most negative real axis it can include in its stability region is an interval of length $2S^2$, the most imaginary axis it can contain is an interval of length $S-1$, and the largest circle tangent to the imaginary axis it can contain has radius $S$ (all of these are mutually exclusive as well).

Answer 3

Reid and David have already answered the technical questions, but I'd like to provide some context to it anyway. Today, for most ODEs you can think of, you can achieve essentially perfect accuracy (say, $10^{-14}$) with most of the good ODE solver packages in essentially no time if evaluating $f(x)$ is free.

There are of course exceptions (very large systems, very stiff systems) but a common sentiment in the community is that the question of designing ODE solvers for "standard" systems is a solved one. Consequently, I think that the question you pose is not a very interesting one -- it addresses a component of ODE solver design that is no longer of importance. This may also explain the lack of literature on the subject.

Answer 4

My reasoning is the following. Since your problems are stiff and not big then most probably you should use some implicit method. If you use implicit method then you solve at least one linear system on each step. It means that if the Jacobian is not structured then each step will take $O(\mathit{dim}^3)$ or $O(\mathit{dim}^2)$ flops on each step (depending on how frequently you LU-decompose the Jacobian), aside from evaluating the RHS.

So the first point is to make sure if your RHS is really more expensive than the underlying linear algebra.

The second point: it is known from literature that solvers based on "expensive" methods (i. e. explicit RK methods) sometimes perform faster than "cheaper" ones (explicit multistep methods).

Summing up, I think that you should not only consider the RHS evaluations count.