C.1 Brief remarks on the mid-point rule numerical technique

Let us consider the generic ordinary differential equation:

$$\frac{d}{dt}x(t)=f(x(t),t) \quad,$$ (C.1)

with $f(x,t):\mathbb{R}^2\rightarrow\mathbb{R}$. Let us consider the time interval $[t,t+\Delta t]$ where $\Delta t$ is the time step. The latter equation can be written at the time instant $t+\Delta t/2$:

$$\frac{dx}{dt}\Big\vert _{t+\Delta t/2}=f(x(t+\Delta t/2),t+\Delta t/2) \quad,$$ (C.2)

Let us develop the function $x(t)$ in Taylor series with respect to the initial point $t$:

$$x(t+\Delta t)=x(t+\Delta t/2)+\frac{dx}{dt}\Big\vert _{t+\Delta t... ...}\Big\vert _{t+\Delta t/2} \frac{\Delta t^2}{4}+\mathcal{O}(\Delta t^3) \quad,$$ (C.3)

$$x(t)=x(t+\Delta t/2)-\frac{dx}{dt}\Big\vert _{t+\Delta t/2} \frac... ...}\Big\vert _{t+\Delta t/2} \frac{\Delta t^2}{4}+\mathcal{O}(\Delta t^3) \quad,$$ (C.4)

where the symbol $\mathcal{O}(\Delta t^3)$ indicates the terms of the third order and greater in $\Delta t$. By subtracting the latter expressions, one obtains:

$$\frac{x(t+\Delta t)-x(t)}{\Delta t}=\frac{dx}{dt}\Big\vert _{t+\Delta t/2} +\mathcal{O}(\Delta t^2) \quad,$$ (C.5)

meaning that the substitution in Eq. (C.2) of the derivative in the mid-point of the interval $[t,t+\Delta t]$ with the left-hand side of Eq. (C.5), leads to a truncation error of the second order with respect to the time step. By summing Eqs. (C.3)-(C.4) and using simple algebra the following mid-point formula can be derived:

$$x(t+\Delta t/2)=\frac{x(t+\Delta t)+x(t)}{2}+\mathcal{O}(\Delta t^2) \quad.$$ (C.6)

In addition, by using Eq. (C.6), it can be shown with the very same line of reasoning that:

$$f\left(x(t+\Delta t/2),t+\frac{\Delta t}{2}\right)=f\left(\frac{... ...Delta t)+x(t)}{2},t+\frac{\Delta t}{2}\right) + \mathcal{O}(\Delta t^2) \quad.$$ (C.7)

Therefore, the numerical scheme obtained from Eq. (C.2) with the second-order approximations (C.5) and (C.7)

$$\frac{dx}{dt}\Big\vert _{t+\Delta t/2} =\frac{x(t+\Delta t)-x(t)}{\Delta t} +\mathcal{O}(\Delta t^2)$$ (C.8)

$$\frac{x(t+\Delta t)+x(t)}{2} =x(t+\Delta t/2) +\mathcal{O}(\Delta t^2)$$ (C.9)

with the position $t_n=t_0+n\Delta t$ and $x_n=x(t_n)$, can be written in the following way:

$$\frac{x_{n+1}-x_n}{\Delta t}=f\left(\frac{x_{n+1}+x_n}{2},t_n+\frac{\Delta t}{2}\right) \quad.$$ (C.10)

This scheme is commonly refereed to as mid-point rule numerical technique and is second-order accurate with respect to the time step $\Delta t$. Now we want to investigate the stability property of the mid-point rule scheme (C.10). We refer, for sake of simplicity, to the scalar initial value problem:

$$\frac{dx}{dt}=\lambda x \quad,\quad x(t_0)=x_0 \quad,\quad\lambda\in\mathbb{C}: [\lambda]<0 \quad.$$ (C.11)

The latter equation can be discretized according to the mid-point rule:

$$x_{n+1}-x_n=\frac{\lambda\Delta t }{2} (x_{n+1}+x_n) \quad.$$ (C.12)

With some straightforward algebra, one can obtain the following time-stepping algorithm:

$$x_{n+1}=\frac{1-\frac{\lambda\Delta t }{2}}{1+\frac{\lambda\Delta t }{2}} x_n \quad.$$ (C.13)

Now, if we study the evolution of two solutions of Eq. (C.11), one starting from the initial condition $x_0$ and the other starting from $y_0=x_0+e_0$, the evolution of the perturbation $e_n=y_n-x_n$ can be found with the same time-stepping as Eq. (C.13):

$$e_{n+1}=\frac{1-\frac{\lambda\Delta t }{2}}{1+\frac{\lambda\Delta t }{2}} e_n \quad,$$ (C.14)

which can be rewritten with respect to the initial perturbation $e_0$:

$$e_{n}=\left[\frac{1-\frac{\lambda\Delta t }{2}}{1+\frac{\lambda\Delta t }{2}}\right]^n e_0 \quad.$$ (C.15)

In particular, the modulus of the perturbation evolves according to the following equation:

$$\vert e_{n}\vert=\left\vert\frac{1-\frac{\lambda\Delta t }{2}}{1+\frac{\lambda\Delta t }{2}}\right\vert^n \vert e_0\vert \quad.$$ (C.16)

It turns out that, in order that the error vanishes for $n\rightarrow\infty$, the following constraint is required:

$$\left\vert\frac{1-\frac{z }{2}}{1+\frac{z}{2}}\right\vert< 1 \quad \quad.$$ (C.17)

where the complex variable $z=\lambda\Delta t$ has been defined. It can be shown that the complex function

$$g(z)=\frac{1-\frac{z }{2}}{1+\frac{z}{2}}$$ (C.18)

fulfills the constraint (C.17) $\forall z\in\mathbb{C}:$Re$[z]<0$. Therefore, the mid-point rule numerical scheme is stable for any $\lambda\in\mathbb{C}:$Re$[\lambda]<0$ and for any time step $\Delta t$, namely is unconditionally stable. In particular, this property is referred in literature [95] to as A-stability of the mid-point rule numerical method.