# Optimal ODE method for fixed number of RHS evaluations

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$).

• The usual correspondence of fixed step Runge-Kutta methods with Newton-Cotes methods applies to the case of the RK method being applied to the IVP $y^\prime=f(x)$; e.g., applying the classical fourth-order method to that IVP is equivalent to performing Simpson's rule on $f(x)$. Commented May 21, 2013 at 13:45
• @J.M.: I'm aware of that. I only intended to use the quadrature rules as an example of characterizing the accuracy of a numerical method for a certain set of inputs when the number of function evaluations is limited. Aside from that I'm interested in "true" ODEs, i.e. those which do not reduce to standard integration. Commented May 21, 2013 at 14:32
• This is getting more interesting. Now the number $N$ by itself doesn't mean anything. What might be helpful is to know $\lambda N/T$, where $T$ is the length of the interval of integration and $\lambda$ is the Lipschitz constant of $f$ with respect to $x$. This will tell us how stiff the problem really is. Assuming it's stiff, a likely candidate is the 2nd order BDF method. Commented May 22, 2013 at 17:17
• @DavidKetcheson: I'm more interested in the general approach for choosing a suitable method for a given problem rather than in the optimal method for a specific problem. We have a larger number of models that vary greatly in stiffness and timing requirements. Commented May 23, 2013 at 7:02
• You say that $f$ is very expensive to evaluate. Can you compute a Jacobian at all? What about some approximation that can correct the principle stiffness? You're not in good shape if your problem is very stiff and you have no way to correct it. Commented May 23, 2013 at 22:38

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).

• Thank you. The paper by Hosea and Shampine is indeed very interesting. Do you know of similar results for stiff problems? I'm aware that one usually uses implicit methods for those, but these have no a priori bound on the number of RHS evaluations, so they are of little use in my case. Commented May 22, 2013 at 6:19
• I don't know of anything like this for stiff problems, but I suspect something exists. Like you say, the question is more subtle when using implicit methods. One approach might be to use Rosenbrock methods, which handle stiff problems well but have a fixed number of RHS evaluations. Commented May 22, 2013 at 6:42

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:

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).

• It can happen that using a variable step (or even variable order) method might be more efficient than doggedly sticking to a fixed step method. One might for instance consider using an extrapolative method like Bulirsch-Stoer: do a few evaluations at some steps, and then build (ostensibly) more accurate estimates from the results of those steps. Commented May 21, 2013 at 14:47
• True. As a matter of fact many of the optimal methods are equivalent in some sense to a variable step version of another time-stepper. Runge-Kutta-Chebshev for example can be seen as forward Euler applied with the variable time-steps being Chebyshev points. Commented May 21, 2013 at 16:22
• @J.M.: Exactly. But is there a way to judge the accuracy of these approaches w.r.t. the number of RHS evaluations, aside from numerical experiments (which would be very involved, given the high number of possible approaches)? Commented May 22, 2013 at 7:16
• @Florian, not in general. You've heard of Lorenz's equations, I presume? Commented May 22, 2013 at 7:17
• @J.M.: Yes :) That's why I mentioned the quadrature example, where accuracy is measured w.r.t. a subset (polynomials) of the original problem space. I'd be happy with results that only work for a certain subset of problems. Commented May 22, 2013 at 7:24

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.

• +1. Whenever somebody asks about efficient ODE solvers, I just assume they're interested in huge systems of ODEs coming from PDE semi-discretizations or large n-body problems. Commented May 22, 2013 at 6:44
• Can you please explain how this relates to my question? I don't see the connection, since I'm interested in the case where evaluating f(x) is not free but rather so expensive that the number of evaluations is limited. Commented May 22, 2013 at 7:02
• @DavidKetcheson: This is not the case here. It's rather that we have very strict timing requirements (hard real-time) on weak hardware (embedded devices). The ODE systems themselves are comparatively small. Commented May 22, 2013 at 7:11
• @FlorianBrucker I think I may have misunderstood your question. It's not that you don't care about function evaluations but that you say the number of function evaluations is bounded. I guess then the question is: Is $N$ large or small? Because the development of ODE solvers has been guided by minimizing $N$ keeping the error constant for decades, and at least for asymptotically large $N$ we know the best methods: high order methods (e.g. RK45) with adaptive time stepping. For these, it is also clear that the larger $N$ the smaller the error. Commented May 22, 2013 at 13:48
• @FlorianBrucker: What I think isn't quite clear is whether the same methods that we know work well for asymptotically large $N$ are also the best for small $N$, say $N<1000$. In other words, you need to tell us something about your budget in terms of how many function evaluations you can afford. Commented May 22, 2013 at 13:50

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.

• You're correct that there's more to execution time and accuracy than the number of RHS evaluations. However, for the other factors (e.g. the linear algebra) there is a lot of material on execution speed vs. problem size. Hence my focus on $N$ in this question. Commented May 24, 2013 at 5:51