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

Question

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:

\begin{equation} r = \frac{||e_l||}{||y_{n+1}||\epsilon_R+\epsilon_A} \tag{4} \\ \end{equation} 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?

Answer 1

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 !