# Step-size selection for an Trapezoidal Method ODE solver (ode23t)

I was reading the documentation of the MatLab ODE solver ode23t, and I've seen that the trapezoidal rule is used. Moreover, the error is estimated by differentiating a cubic polynomial, and that's still fine to me.

My big question is about the step size selection: most of ODE solver usually are embedded Runge-Kutta methods and the step-size selection is done (modulus lots of improvements made by MatLab developers) in a way that is presented in classical texts like "Hairer-Wanner-Norseet (I and II volume)".

What can one do, in order to change the step-size, if only trapezoidal rule is used? I'm really curious about that.

EDIT: Thanks @Chris

Still a doubt

I've read the SPICE book (the interesting pages are from 308 to 312), and the algorithm at page 310 is also clear: it uses the number of iteration required to solve the non-linear system in order to change the step-size. But... as stated above, MatLab estimates the error by differentiating a cubic interpolant. What are the points (and the derivatives) used to build this interpolant? Your code seems to use the solution just computed $$u_{n+1}$$, and the $$u_{n}, u_{n-1}$$, but still can't see where are the derivatives.

• Some solvers use an adaptive step-size that can change over time. – R zu Oct 29 '18 at 17:46

The only documentation I know about for the implementation of ode23t is in the paper which documents the implementation of ode23tb, the TRBDF2 method in MATLAB.

As usually implemented, the trapezoidal rule is not strictly a one-step method because the truncation error estimate makes use of two previously computed solution values. The method is not efficient for very stiff problems because it is not strongly stable.

and it cites The SPICE book. Two things here. First of all, the Trapezoidal rule, also known as Crank-Nicholson, is not L-stable and is generally not suited for highly stiff equations (but this lack of dampening can give undesirable outcomes in some models). This is exacerbated by the instability of the error estimate, which is from a truncation error estimate. Basically, you save some previous points and use derivative approximations in order to match the leading truncation error terms, and then use this as the input to standard P or PI adaptive algorithms.

For an open source (MIT licensed) implementation of this algorithm, you can take a look at the Trapezoid method from DifferentialEquations.jl. The implementation can be found in the source code here and it was derived from the description in the SPICE book reference given by Shampine. Note that due to licensing issues this has not been directly checked against the MATLAB source code, only directly against these sources given by Shampine.

## Responding to Edits

Continuing the conversation is kind of frowned upon in SO, but I'll do it since it makes sense. For part one, the linked method isn't implicit Euler in the first step. Instead, the adaptivity doesn't exist in the first step since it uses two steps to get a second derivative approximation. However, being a 1-step method, the Trapezoidal rule is still used even when it takes the non-adaptive step. This keeps it always second order (though you can have a single first order step and still be second order).

But... as stated above, MatLab estimates the error by differentiating a cubic interpolant.

Where is it stated that MATLAB does this? MATLAB uses a cubic interpolant error for dde23, but I have never seen it stated that it does this on ode23t. In the comments I show the paper trail starting from MATLAB's documentation to Shampine's paper and to the SPICE book. In the SPICE book they derive the truncation error estimate which uses two steps to estimate the third derivative term which is the leading truncation error term. To the best of my knowledge, and not looking at the MATLAB source since I am an open source developer and do not want licensing issues, that is how it's done according to the sources given. If someone has the time, they can double change and confirm or deny this by looking at the MATLAB source.

To answer the question of how you could come up with an error estimate using the interpolant, let me show you that Hermite won't work here. I am just going to point to the live code since I know it's correct and it's easier than TeXing it up. The Hermite interpolant is defined here:

with ($$y_0$$,k[1]) and ($$y_1$$,k[2]) as the values at t and t+dt respectively, calculating at the point t+$$\theta$$dt. To get the derivative of the interpolant, you differentiate by $$\theta$$ and get
Notice that when $$\theta=1$$, the derivative approximated by the cubic Hermite polynomial is $$k[2]$$, which is the derivative estimate at $$t+dt$$ calculated from $$f$$ itself. So it doesn't have a natural error estimate at the end points of the interval against the derivative since it's made to perfectly hit those derivative values. When this estimate is used it's called a residual correction method, discussed here, where this derivative estimate is used against extra $$f$$ evaluations in the middle of the interval to give a continuous derivative error estimate. By using specific points you can estimate the maximum difference, and this done on RK4 is how ddesd works. RK4 is simple, and an open source implementation of the two-point maximal residual error estimate for it can be found here. But while MATLAB does use this in a different integrator, Shampine writes about it, and you technically can apply this to the Trapezoid method, from what I can tell this is not what's done in ode23t.
• I wouldn't have guessed that the Hosea and Shampine paper documents MATLAB's ode23tb implementation specifically, rather than the TR-BDF2 method in general. Do you have a reference for that statement? – Rahul Oct 29 '18 at 7:43
• @ChrisRackauckas it's almost everything clear, except one thing: we you say "in the SPICE book they derive the truncation error estimate which uses two steps to estimate the third derivative term [...] " , how can it use just two steps? Also in the source code you linked (here) it seems that four values of the solution $u(t)$ are required: u, uprev, uprev2,uprev3. What am I missing? – VoB Oct 30 '18 at 19:03