Numerical method

Two-scale formulation

First, rewrite equation (1) using the variable change $w(t)=\exp(-(t-t_{start})A/\varepsilon) u(t)$ to obtain

\[\frac{d w(t)}{dt} = F\Big(\frac{t-t_{start}}{\varepsilon}, w(t) \Big), \;\;\; w(t_{start})=u_{in}, \;\; \varepsilon\in ]0, 1],\]

where the function $F$ is expressed from the data of the original problem (1)

\[F\Big( \frac{s}{\varepsilon}, w \Big) = \exp(-sA/\varepsilon) \; f( \exp(sA/\varepsilon), \; w).\]

We then introduce the function $U(t, \tau), \tau\in [0, 2 \pi]$ such that $U(t, \tau=(t-t_{start})/\varepsilon) = w(t)$. The two-scale function is then the solution of the following equation

\[\frac{\partial U}{\partial t} + \frac{1}{\varepsilon} \frac{\partial U}{\partial \tau} = F( \tau, U), \;\;\; U(t=t_{start}, \tau)=\Phi(\tau), \;\; \varepsilon\in ]0, 1], \;\;\;\;\;\;\;\;\;\; (2)\]

where $\Phi$ is a function checking $\Phi(\tau=0)=u_{0}$ chosen so that the $U$ solution of (2) is smooth (see Philippe Chartier, Nicolas Crouseilles, Mohammed Lemou, Florian Méhats (2015) and Philippe Chartier, Mohammed Lemou, Florian Méhats, Xiaofei Zhao (2020).

The numerical method is based on a discretization of equation (2). In the direction $\tau$, a spectral method is used, while for the time $t$, an exponential Adams-Bashforth method allows to build a high order method (see Philippe Chartier, Mohammed Lemou, Florian Méhats, Xiaofei Zhao (2020)). The initialization is based on a "butterfly" technique (going back and forth around the initial time).

Initialization

Let r be the order of the method $AB_r$.
Let $\Delta t$ the time step, for $i \in \{r, -(r-1), \ldots, r-1, r\}$, we note $u_i = u(t_{start}+i \Delta t)$.
Let $r'$ be the orders of the intermediate AB methods we will use.
If $u_{k}$ is known with a precision of ${\mathcal O}(\Delta t^{r'+1})$, and for $r' \geq 2, u_{k-1}, \ldots, u_{k-r'+1}$ are known with a precision of ${\mathcal O}(\Delta t^{r'})$ then we can calculate $u_{k+1}$ with a precision of ${\mathcal O}(\Delta t^{r'+1})$ with the method $AB_{r'}$.
Similarly, if $u_{k}$ is known with a precision of ${\mathcal O}(\Delta t^{r'+1})$, and for $r' \geq 2, u_{k+1}, \ldots, u_{k+r'-1}$ are known with a precision of ${\mathcal O}(\Delta t^{r'})$ then we can calculate $u_{k-1}$ with a precision of ${\mathcal O}(\Delta t^{r'+1})$ with the method $AB_{r'}$.

Algorithm

With the method $AB_1$, from $u_0$ we calculate $u_{-1}$ with a precision of ${\mathcal O}(\Delta t^2)$
With the method $AB_2$, starting from $u_{0}$ and $u_{-1}$, we calculate $u_{1}$ with a precision of ${\mathcal O}(\Delta t^3)$
For $r' = $3 to $r' = r$.
- For $k=1$ to $k=r'-1$
  - With the method $AB_{r'-1}$, from $u_{1-k}, u_{2-k}, \ldots,u_{r'-1-k}$, we calculate $u_{-k}$ with a precision of ${\mathcal O}(\Delta t^{r'})$
- For $k=1$ to $k=r'-1$
  - With the method $AB_{r'}$, from $u_{k-1}, u_{k-2}, \ldots,u_{k-r'}$, we calculate $u_{k}$ with a precision of ${\mathcal O}(\Delta t^{r'+1})$

At the end of this algorithm, the values $u_0, u_1, \ldots u_{r-1}$ are known with a precision of ${\mathcal O}(\Delta t^{r+1})$, we can launch the algorithm $AB_r$. Note that this technique requires that the function $f$ has to be smooth enough (namely $f\in {\mathcal C}^r([t_{start} - r \Delta t, t_{end}])$).

The Adams-Bashforth Method

For an Adams-Bashforth method of order $r$ in time and spectral $\tau$, we first introduce a mesh in the $\tau$ direction.
$\tau_{\ell} = \ell \Delta \tau, \ell = 0, \ldots, N_{\tau}-1$. Where $N_{\tau}$ is the number of points of discretization. If we apply the Fourier transform to the two-scale equation, we obtain

\[\frac{\partial \hat{U}_\ell}{\partial t} + \frac{i\ell}{\varepsilon}\hat{U}_\ell = \hat{F}_\ell(t), \;\; \ell=-N_\tau/2, \dots, N_\tau/2-1,\]

with

\[U(t, \tau_k) = \sum_{\ell=-N_{\tau}/2}^{N_{\tau}/2-1} \hat{U}_{\ell}(t) e^{i\ell k 2\pi/N_{\tau}} \;\;\; \text{ and } F(\tau_k, U(t, \tau_k)) = \sum_{\ell=-N_{\tau}/2}^{N_{\tau}/2-1} \hat{F}_{\ell}(t) e^{i\ell k 2\pi/N_\tau}.\]

and the inverse (discrete) Fourier transform formulae

\[\hat{U}_{\ell}(t) = \frac{1}{N_{\tau}}\sum_{k=0}^{N_{\tau}-1} U(t, \tau_k) e^{-i\ell k 2\pi/N_{\tau}}\;\;\; \text{ and } \hat{F}_{\ell}(t) = \frac{1}{N_{\tau}}\sum_{k=0}^{N_{\tau}-1} F(\tau_k, U(t, \tau_k))e^{-i\ell k 2\pi/N_{\tau}}.\]

If we wish to calculate $\hat{F}_{\ell}$ from $\hat{U}_{\ell}$ we have the following formula

\[F(\tau_k,U(t, \tau_k)) = e^{-\tau_k A}f(e^{\tau_k A}U(t, \tau_k)) \;\;\]

from which the Fourier transform is calculated in $\tau$ from the discrete Fourier transform formulas above.

Now, given a time step $\Delta t>0$ and a discretization in time $t_n=n\Delta t$, we can write the following Duhamel formula ($n\geq 0$)

\[\hat{U}_{\ell}(t_{n+1}) = e^{-i\ell\Delta t/\varepsilon}\hat{U}_{\ell}(t_{n}) + \int_0^{\Delta t} e^{-i\ell(\Delta t -s)/\varepsilon} \hat{F}_\ell(t_n+s)ds.\]

Thus, to obtain a $(r+1)$ scheme, we can approaches the function $\hat{F}_\ell(t_n+s)$ by the Lagrange polynomial of order $r$ interpolator at points $t_{n-j}, j=0, \dots, r$. This polynomial is written

\[\hat{F}_\ell(t_n+s) \approx \sum_{k=0}^r \Big(\Pi_{j=0, j\neq k}^r \frac{s+j \Delta t}{(j-k)\Delta t} \Big) \hat{F}_\ell(t_n-t_j), \;\; n\geq 0.\]

Thus, from this approximation, we integrate exactly, which requires the following formulas

\[p^{[r]}_{\ell, j} = \int_0^{\Delta t}e^{-i\ell(\Delta t -s)/\varepsilon}\Big( \Pi_{j=0, j\neq k}^r \frac{s+j \Delta t}{(j-k)\Delta t}\Big) ds,\]

for each $j$ and $\ell$ such that $0\leq j\leq r, \; \ell=-N_\tau/2, \dots, N_\tau/2-1$. These coefficients $p^{[r]}_{\ell, j}$ can be pre-calculated and stored once and for all. Thus, the schema is finally written

\[\hat{U}_{\ell}^{n+1}= e^{-i\ell\Delta t/\varepsilon}\hat{U}_{\ell}^n + \sum_{j=0}^r p^{[r]}_{\ell, j} \hat{F}_\ell^{n-j},\]

with $\hat{U}_{\ell}^n \approx \hat{U}_{\ell}(t_n)$ and $\hat{F}_\ell^{n-j}\approx \hat{F}_\ell(t_{n-j})$.

We can verify that the truncation error in this schema is ${\mathcal O}(\Delta t^{r+1})$, once the initial values $\hat{U}_\ell^1, \dots, \hat{U}_\ell^r$ have been calculated.

Non-homogeneous case $f(t, u)$

Here we consider the case where $f$ depends on the variable $t$.

\[\frac{d u(t)}{dt} = \frac{1}{\varepsilon} A u(t) + f(t, u(t)), \;\;\; u(t=t_{start})=u_{in}, \;\; \varepsilon\in ]0, 1] \;\;\;\; (3)\]

The non-homogeneous case (3) falls under (1), by entering the variable $\theta : t \in [t_{start}, t_{end}] \mapsto \theta(t) = t\in \mathbb{R}$ which allows us to reformulate the non-homogeneous case (3) into a homogeneous problem of the form (1). Indeed, we rephrase (3) as follows

\[\begin{aligned} \frac{d u(t) }{dt} & = \frac{1}{\varepsilon} Au(t) + f(\theta(t), u(t)), \\ \frac{d \theta(t) }{dt} & = 1 \end{aligned}\;\;\;\;(4)\]

with the initial condition $u(t_{start})=u_{in}, \theta(t_{start})=t_{start}$. Thus, the problem (4) is rewritten into an equation
satisfied by $y: t\in [t_{start}, t_{end}] \mapsto y(t) =(u(t), \theta(t))\in \mathbb{R}^{n+1}$

\[\frac{d y}{dt} = \frac{1}{\varepsilon} \tilde{A} y + g(y), \;\; y(t_{start})=(u_{in}, t_{start}),\]

with $\tilde{A}\in{\mathcal M}_{n+1, n+1}(\mathbb{R})$