Exponential and integrating factor integrators — Examples

For coordinate $i$ with ${\lambda }_i\in \{1,1000\}$, put $a={\lambda }_i+1>0$. The exact equation is $y_i^{\mathrm{\prime }}=-a{y}_i+1$, $y_i(0)=0$. The product rule [F6] and chain rule [F7] show $(e^{at}{y}_i(t){)}^{\mathrm{\prime }}=e^{at}$. Integrating by [F4] gives $e^{ah}{y}_i(h)={\int }_0^h{e}^{as}ds=(e^{ah}-1)\mathrm{/}a$, since [F3], [F7] give antiderivative $e^{as}\mathrm{/}a$. Dividing by the positive exponential using [F5] yields $y_i(h)=(1-e^{-ah})\mathrm{/}a$.

By [F8], powers of $A=\mathrm{diag}(-{\lambda }_i)$ are diagonal with entries $(-{\lambda }_i{)}^j$, so the series [F2] evaluates the matrix phi functions coordinatewise. The exact exponential-Euler map has $e^{-{\lambda }_{i}h}{y}_i(0)=0$ and reaction value $g_i(0)=1$, so its coordinate is $h{\phi }_1(-{\lambda }_{i}h)$. The ETD2 predictor has that same coordinate; hence $g_i(a)=1-h{\phi }_1(-{\lambda }_{i}h)$ and the ETD2 coordinate is $h[({\phi }_1(-{\lambda }_{i}h)-{\phi }_2(-{\lambda }_{i}h))+{\phi }_2(-{\lambda }_{i}h)(1-h{\phi }_1(-{\lambda }_{i}h))]=h({\phi }_1(-{\lambda }_{i}h)-h{\phi }_1(-{\lambda }_{i}h){\phi }_2(-{\lambda }_{i}h))$.

In particular the stiff coordinate uses ${\lambda }_{2}h=1$: its exact, exponential-Euler, and ETD2 values are respectively $(1-e^{-1.001})\mathrm{/}1001$, ${10}^{-3}{\phi }_1(-1)$, and ${10}^{-3}({\phi }_1(-1)-{10}^{-3}{\phi }_1(-1){\phi }_2(-1))$. All exponential and phi actions are at explicit matrix/scalar arguments; any numerical implementation must return the certificates specified by Certified-action implementation of exponential Euler and ETD2.

Write $A=-I+N$ with $N=\left(\begin{array}{cc}0& 2e\\ 0& 0\end{array}\right)$; direct multiplication gives $N^2=0$ [F4]. The formula holds for $n=1$ since $A=-I+N$. If $A^n=(-1{)}^{n}I+n(-1{)}^{n-1}N$, multiplying by $A=-I+N$ and using $N^2=0$ gives $A^{n+1}=(-1{)}^{n+1}I+(n+1)(-1{)}^{n}N$, so induction proves it for every $n\ge 1$. Substituting in the absolutely convergent matrix-exponential series [F2] gives $e^{tA}=e^{-t}I+{\sum }_{n\ge 1}{t}^n(-1{)}^{n-1}N\mathrm{/}(n-1)!=e^{-t}(I+tN)$ by the scalar series [F3].

By [F7] and step 1.1, the solution is $y(t)=e^{-t}(2et,1{)}^{\mathsf{T}}$. Direct multiplication using [F4] gives $A{A}^{\ast }=\left(\begin{array}{cc}1+4e^2& -2e\\ -2e& 1\end{array}\right)$ and $A^{\ast }A=\left(\begin{array}{cc}1& -2e\\ -2e& 1+4e^2\end{array}\right)$, which differ, so $A$ is nonnormal. By [F5] and [F6], its initial norm is $1$, whereas $\mathrm{\parallel }y(1){\mathrm{\parallel }}_2=e^{-1}\sqrt{4e^2+1}>2$. For $t>0$, the exponential series gives $e^t\ge t^2\mathrm{/}2$, and [F6] gives $e^{-t}=1\mathrm{/}{e}^t\le 2\mathrm{/}{t}^2$; hence $e^{-t}\to 0$ and $t{e}^{-t}\le 2\mathrm{/}t\to 0$, so $y(t)\to 0$.

By [F1], every exact-action step is multiplication by $e^{hA}$. Repeated use of the semigroup identity in [F2] gives $(e^{hA}{)}^n=e^{nhA}$, so the method returns $y(t_n)$ at every grid point. The transient growth is therefore present in the exact linear flow and is not removed by exponentiating this nonnormal matrix.

By [F3], the exact solution is $y(t)={\int }_0^{t}sds=t^2\mathrm{/}2$, so $y(h)=h^2\mathrm{/}2$.

Since $A=0$, $e^{hA}=1$ and $g(0,y(0))=0$. The exponential-Euler step is therefore $0+h{\phi }_1(0)0=0$, and its predictor for ETD2 is also $a=0$.

At the ETD2 endpoint, $g(h,a)=h$. Its update is $0+h({\phi }_1(0)-{\phi }_2(0))0+h{\phi }_2(0)h=h^2\mathrm{/}2$, which equals the exact value from step 1.1. This calculation evaluates ${\phi }_k(0)$ from the series and never divides by $A$ or its eigenvalues.

The matrix series [F5] for $A_{\mathrm{i}\mathrm{m}}=-H$ is diagonal term by term and its diagonal series are the real exponential series from [F7], so $e^{\tau A_{\mathrm{i}\mathrm{m}}}=\mathrm{diag}(e^{-\tau },e^{-2\tau })$. By [F6], $u(\tau )=(e^{-\tau },e^{-2\tau }{)}^{\mathsf{T}}$. For $\tau >0$, [F7] gives $0<e^{-\tau }<1$ and $0<e^{-2\tau }<1$, hence $\mathrm{\parallel }u(\tau ){\mathrm{\parallel }}_2^2=(e^{-\tau }{)}^2+(e^{-2\tau }{)}^2<2=\mathrm{\parallel }u(0){\mathrm{\parallel }}_2^2$; this imaginary-time propagator is not unitary.

The matrix series [F5] for $A_{\mathrm{r}\mathrm{t}}=-iH$ is diagonal term by term and agrees with [F4], so $e^{t{A}_{\mathrm{r}\mathrm{t}}}=\mathrm{diag}(e^{-it},e^{-2it})$. Since $A_{\mathrm{r}\mathrm{t}}^{\ast }=-A_{\mathrm{r}\mathrm{t}}$, [F1] gives a unitary propagator; [F6] yields $y(t)=(e^{-it},e^{-2it}{)}^{\mathsf{T}}$ and $\mathrm{\parallel }y(t){\mathrm{\parallel }}_2=\sqrt{2}$.

By [F2], both exact-action methods reduce to multiplication by $e^{h{A}_{\mathrm{r}\mathrm{t}}}$. The semigroup identity [F5] gives $(e^{h{A}_{\mathrm{r}\mathrm{t}}}{)}^n=e^{nh{A}_{\mathrm{r}\mathrm{t}}}$, so every grid value is the exact $y(t_n)$ from step 1.2 and has norm $\sqrt{2}$.

By [F8], $\mathbb{R}\times {\mathbb{C}}^2$ is identified with ${\mathbb{R}}^5$ using the realification and product Euclidean norm. For a finite certified real-time run, take $D=\mathbb{R}\times {\mathbb{C}}^2$ and $g=0$. Choose $R>\sqrt{2}$ larger than the absolute values of every real and imaginary coordinate of its finitely many numerical states and ETD2 predictor stages. The box $\mathcal{T}=[0,T]\times [-R,R{]}^4$ contains those states and stages, and the exact path from step 1.2 lies in it because its norm is $\sqrt{2}$. By [F9], $\mathcal{T}$ is compact; since $D$ is the whole time-state space, it has positive margin. Here $M=1$, $\omega =0$, $K_{+}=K_{-}=1$, $L_g=Q_1=Q_2=0$, and therefore $C_p=L_p=0$ in [F3]. With exact starting data, [F3] gives ${\max }_n\mathrm{\parallel }{\tilde{y}}_n-y(t_n){\mathrm{\parallel }}_2\le {\sum }_{j=0}^{N-1}{\delta }_j$.