# How to select initial time step in adaptive time step ODE solver (TR-BDF2)

The Problem

I am currently reconstructing a TR-BDF2 scheme which contains the following two stages:

\begin{align} y_{n+\gamma} & = y_n + \gamma \frac{h}{2}\left( f_n + f_{n+\gamma} \right) \tag{1} \\ % y_{n+1} & = \frac{1}{\gamma(2-\gamma)}y_{n+\gamma} - \frac{(1-\gamma)^2}{\gamma(2-\gamma)}y_n + \frac{1-\gamma}{2-\gamma}hf_{n+1} \tag{2} \end{align}

From those, the local truncation error is derived as:

$$\begin{gather} e_l = 2k_\gamma \Delta t \left( \frac{1}{\gamma}f _n - \frac{1}{\gamma(1-\gamma)}f_{n+\gamma} + \frac{1}{1-\gamma} f_{n+1} \right), \ \text{where} \ k_\gamma = \frac{-3\gamma^2+4\gamma-2}{12(2-\gamma)} \tag{3} \end{gather}$$

Based on the above, a recommended method to calculate the next time step $$h$$ which I found in these lecture notes, would be via the below formula:

$$$$r = \frac{||e_l||}{||y_{n+1}||\epsilon_R+\epsilon_A} \tag{4} \\$$$$ where $$\epsilon_R$$ and $$\epsilon_A$$ are the user-set relative and absolute tolerances respectively.

1. if $$r\leq2$$ accept the solution $$y_{n+1}$$ and set $$h_{n+1}=h_n/r^{\frac{1}{p+1}}$$.
2. else redo the step by setting a new timestep $$h_{redo}=h_n/r^{\frac{1}{p+1}}$$.

where $$p=2$$.

The question

The above seems fine to me however my question is, what would be a rule of thumb in order to derive the initial time step that the method has to take?

Many numerical tips and theoretical explanations can be found in this book from Hairer and Wanner: https://www.springer.com/gp/book/9783540566700

In this book, a strategy is described, which uses a time step such that the relative variation of the solution during the first time step is below a certain threshold if you were using explicit Euler (omitting the time dependency of $$f$$ for simplicity): $$y_1 = y_0 + \Delta t f(y_0)$$ The variation is $$\Delta t f(y_0)$$. You can choose $$\Delta t$$ such that $$\Delta t |f(y_0)| < atol + rtol |y_0|$$, or a fraction of this value.

Some authors have suggested that the local error can be assumed to behave as $$e \approx C \Delta t^{p+1} \frac{d^{p+1} y}{dt^{p+1}}$$ for a method of order $$p$$. Therefore, a slightly more advanced solution is the following (page 169 of the book):

• Compute a first time step value $$\Delta t_0 = 10^{-2} \dfrac{atol + rtol||y_0||}{atol + rtol ||f(y_0)||}$$, which should already ensure that an explicit Euler step will only let the solution vary by ~1%.
• run one explicit Euler step: $$y_1 = y_0 + \Delta t_0 f(y_0)$$
• compute $$f(y_1)$$
• an estimate of the second time derivative of $$y$$ is then: $$d =\frac{||f(y_1) - f(y_0)||}{\Delta t_0}$$
• Therefore, to ensure we get a small enough error if we used explicit Euler, we can use $$\Delta t_1$$ such that: $$\Delta t_1^{2} ||d|| < 10^{-2}$$. Here the factor $$10^{-2}$$ is taken as a small enough factor to compensate the fact that the error constant $$C$$ is not given.

Additionaly, the authors suggest to choose the can take the minimum of $$(100\Delta t_0, \Delta_ 1)$$. They also suggest to extend this by considering the order $$p$$ of your method instead, and computing $$\Delta t_1$$ such that $$\Delta t_1^{p+1} ||d|| < 10^{-2}$$, i.e. they replace the (p+1)-th derivative of $$y$$ with its second derivative.

Otherwise, you can also estimate the dominant eigenvalue $$\lambda$$ of your system and chose $$\Delta t$$ such that $$|\lambda \Delta t| \ll 1$$ to ensure your method is used in a "zone" where it is precise. The eigenvalue can be estimated as: $$\lambda \approx \frac{||f(y_1) - f(y_0)||}{||y_1-y_0||}$$.

Personally, I often use a first step of $$10^{-8}$$ and let the time step adaptation algorithm increase $$\Delta t$$ from there, which should happen fairly quickly. Of course this may not be very efficient for large systems. Anyway you will most likely run multiple simulations of a given system, therefore you'll be able to make an educated guess after only a few !

• This method is somewhat unstable though. It's really derived for explicit methods. In DifferentialEquations.jl, we still use this anyways, but with a failure case that if it asks for something too small then you just use 1e-6. And FWIW, Fortran radau's default is to always start 1e-6. – Chris Rackauckas Nov 18 '20 at 23:23
• It's almost always a good idea to refer to Hairer & Wanner :-) – Wolfgang Bangerth Nov 19 '20 at 6:12
• @Laurent90 as far as I understand, regardless of whether we are using explicit methods, to derive the initial timestep, the Euler step is used as a rule of thumb, at the risk of being somewhat unstable as it was pointed out above, is that correct? – kostas1335 Nov 19 '20 at 7:45
• Yes this is the spirit. However it may indeed be too pessimist and let you choose an unnecessarily low value of the initial time step. I didn't write it here, but Hairer and Wanner also describe on the same page some "default" behaviour as @ChrisRackauckas mentionned. Lastly, you can also take a the initial time step as a fraction of some characteristic time of your system. – Laurent90 Nov 19 '20 at 8:35