Can stiff ODE methods have an interpolant of the same order?

Are there implicit methods that generate a continuous extension (i.e. dense output) at the same order as the integration? Many non-stiff solvers have this property, but I don’t know of any 4th or 5th order implicit solvers that have the same order interpolant.

  1. If the local error is O(h^{p+1}) then you can still achieve the same global accuracy after O(1/h) steps with an O(h^{p}) interpolant [1]. Maybe you just don’t need it?
  2. On the other hand, there is a classic result showing the equivalence between implicit methods and collocation. So the p-order solver is equivalently estimating some p-order collocation polynomial, with nodes from the left column of the Butcher tableau. Some polynomial exists!

Is it just too inefficient to estimate the function values at the nodes?


Yes, there are the general linear methods of John Butcher. I’ve written code for methods of orders 2 thru 5 in Python. I’ve not gotten around to porting them the Julia yet.


1 Like

Not all methods are collocation methods. RadauIIA5 and BDF have collocation polynomials, but not every SDIRK is guaranteed to have one.

But the bigger deal is that most of the special interpolations for stiff solvers use some kind of ki = W \ bi k’s for building the polynomial since the k’s give a smooth interpolation of the smooth manifold. Directly using the evaluation points can be less stable. That then reduces the number of points you have access to similarly to the explicit RK methods, and in general it’s then much easier to construct one with error one order less.

1 Like

Thank you both, the two answers here resolve my question. They exist, but are not easy to construct. Perhaps I should be more impressed that explicit RK methods like Verner/Feagin can have a same-order interpolant without too many extra stages.

I would speculate that implicit schemes could be efficiently extended with calculations between stages by re-using factorized Jacobians, to produce a higher order interpolant. Maybe the benefits for delay and integro-differential equations don’t outweigh the extra work per step.

And you couldn’t do it lazily because otherwise you’d have to store the factorized Jacobian (or recompute it)

1 Like