Comparison of numerical methods for the nonlinear Klein-Gordon equation in the nonrelativistic limit regime

TL;DR

This paper systematically compares various numerical methods for solving the nonlinear Klein-Gordon equation in the nonrelativistic limit, focusing on efficiency, accuracy, and scalability as the parameter approaches zero.

Contribution

It provides a comprehensive comparison of multiple numerical methods for the NKGE, highlighting their performance and $ ext{ε}$-resolution capabilities in the nonrelativistic limit regime.

Findings

Multiscale time integrator shows superior $ ext{ε}$-resolution.

Different methods exhibit varying efficiency and accuracy in the nonrelativistic limit.

Convergence rates of the NKGE to limiting models are analyzed.

Abstract

Different efficient and accurate numerical methods have recently been proposed and analyzed for the nonlinear Klein-Gordon equation (NKGE) with a dimensionless parameter $\varepsilon\in (0,1]$, which is inversely proportional to the speed of light. In the nonrelativestic limit regime, i.e. $0<\varepsilon\ll1$, the solution of the NKGE propagates waves with wavelength at $O(1)$ and $O(\varepsilon^2)$ in space and time, respectively, which brings significantly numerical burdens in designing numerical methods. We compare systematically spatial/temporal efficiency and accuracy as well as $\varepsilon$-resolution (or $\varepsilon$-scalability) of different numerical methods including finite difference time domain methods, time-splitting method, exponential wave integrator, limit integrator, multiscale time integrator, two-scale formulation method and iterative exponential integrator. Finally, we adopt the multiscale time integrator to study the convergence rates from the NKGE to its limiting models when $\varepsilon\to0^+$.

Figures

Equations

$$\{\begin{array}{rl}& {\epsilon }^2{\partial }_{tt}u(\mathbf{x},t)-\Delta u(\mathbf{x},t)+\frac{1}{{\epsilon }^2}u(\mathbf{x},t)+\lambda |u(\mathbf{x},t){|}^{2}u(\mathbf{x},t)=0,\mathbf{x}\in {\mathbb{R}}^d,t>0,\\ & u(\mathbf{x},0)={\varphi }_1(\mathbf{x}),{\partial }_{t}u(\mathbf{x},0)=\frac{1}{{\epsilon }^2}{\varphi }_2(\mathbf{x}),\mathbf{x}\in {\mathbb{R}}^d,\end{array}$$

$$E(t)$$

$$\equiv {\int }_{{\mathbb{R}}^d}\left[\frac{1}{{\epsilon }^2}|{\varphi }_2(\mathbf{x}){|}^2+|\nabla {\varphi }_1(\mathbf{x}){|}^2+\frac{1}{{\epsilon }^2}|{\varphi }_1(\mathbf{x}){|}^2+\frac{\lambda }{2}|{\varphi }_1(\mathbf{x}){|}^4\right]d\mathbf{x}=E(0),t\ge 0.$$

$${\varphi }_1(x)=\frac{3\sin (x)}{e^{x^2/2}+e^{-x^2/2}},{\varphi }_2(x)=\frac{2e^{-x^2}}{\sqrt{\pi }},x\in \mathbb{R}.$$

$$u(\mathbf{x},t)=e^{it/{\epsilon }^2}z(\mathbf{x},t)+e^{-it/{\epsilon }^2}\overline{z}(\mathbf{x},t)+O({\epsilon }^2),\mathbf{x}\in {\mathbb{R}}^d,t\ge 0,$$

$$\{\begin{array}{rl}& 2i{\partial }_{t}z(\mathbf{x},t)+{\epsilon }^2{\partial }_{tt}z(\mathbf{x},t)-\Delta z(\mathbf{x},t)+3\lambda |z(\mathbf{x},t){|}^{2}z(\mathbf{x},t)=0,\mathbf{x}\in {\mathbb{R}}^d,t>0,\\ & z(\mathbf{x},0)=\frac{1}{2}\left[{\varphi }_1(\mathbf{x})-i{\varphi }_2(\mathbf{x})\right]=:z_0(\mathbf{x}),\mathbf{x}\in {\mathbb{R}}^d,\\ & {\partial }_{t}z(\mathbf{x},0)=\frac{i}{2}\left[-\Delta z_0(\mathbf{x})+3\lambda |z_0(\mathbf{x}){|}^2{z}_0(\mathbf{x})\right].\end{array}$$

$$\{\begin{array}{rl}& 2i{\partial }_{t}z(\mathbf{x},t)-\Delta z(\mathbf{x},t)+3\lambda |z(\mathbf{x},t){|}^{2}z(\mathbf{x},t)=0,\mathbf{x}\in {\mathbb{R}}^d,t>0,\\ & z(\mathbf{x},0)=\frac{1}{2}\left[{\varphi }_1(\mathbf{x})-i{\varphi }_2(\mathbf{x})\right]:=z_0(\mathbf{x}),\mathbf{x}\in {\mathbb{R}}^d.\end{array}$$

$$\{\begin{array}{rl}& {\epsilon }^2{\partial }_{tt}u(x,t)-{\partial }_{xx}u(x,t)+\frac{1}{{\epsilon }^2}u(x,t)+\lambda (u(x,t){)}^3=0,x\in \Omega ,t>0,\\ & u(x,0)={\varphi }_1(x),{\partial }_{t}u(x,0)=\frac{1}{{\epsilon }^2}{\varphi }_2(x),x\in \overline{\Omega },\\ & u(a,t)=u(b,t),{\partial }_{x}u(a,t)={\partial }_{x}u(b,t).\end{array}$$

$$\Vert u^n{\Vert }_{l^2}^2=h\sum _{j=0}^{N-1}|u_j^n{|}^2,n\ge 0.$$

$${\delta }_t^2{u}_j^n=\frac{u_j^{n+1}-2u_j^n+u_j^{n-1}}{{\tau }^2},{\delta }_x^{+}{u}_j^n=\frac{u_{j+1}^n-u_j^n}{h},{\delta }_x^2{u}_j^n=\frac{u_{j+1}^n-2u_j^n+u_{j-1}^n}{h^2}.$$

$${\epsilon }^2{\delta }_t^2{u}_j^n+\left[-\frac{1}{2}{\delta }_x^2+\frac{1}{2{\epsilon }^2}+\frac{\lambda }{4}[(u_j^{n+1}{)}^2+(u_j^{n-1}{)}^2]\right](u_j^{n+1}+u_j^{n-1})=0,0\le j\le N-1,n\ge 0.$$

$${\epsilon }^2{\delta }_t^2{u}_j^n-\frac{1}{2}{\delta }_x^2(u_j^{n+1}+u_j^{n-1})+\frac{1}{2{\epsilon }^2}(u_j^{n+1}+u_j^{n-1})+\lambda (u_j^n{)}^3=0,0\le j\le N-1,n\ge 0;$$

$${\epsilon }^2{\delta }_t^2{u}_j^n-{\delta }_x^2{u}_j^n+\frac{1}{{\epsilon }^2}{u}_j^n+\lambda (u_j^n{)}^3=0,0\le j\le N-1,n\ge 0.$$

$$\begin{array}{rl}& u_0^n=u_N^n,u_{-1}^n=u_{N-1}^n,n\ge 0;u_j^0={\varphi }_1(x_j),j=0,\dots ,N,\\ & u_j^1={\varphi }_1(x_j)+\sin \left(\frac{\tau }{{\epsilon }^2}\right){\varphi }_2(x_j)+\frac{\tau }{2}\sin \left(\frac{\tau }{{\epsilon }^2}\right)\left[{\delta }_x^2{\varphi }_1(x_j)-\frac{1}{\tau }\sin \left(\frac{\tau }{{\epsilon }^2}\right){\varphi }_1(x_j)-\lambda ({\varphi }_1(x_j){)}^3\right].\end{array}$$

$$u_j^1={\varphi }_1(x_j)+\frac{\tau }{{\epsilon }^2}{\varphi }_2(x_j)+\frac{{\tau }^2}{2{\epsilon }^2}\left[{\delta }_x^2{\varphi }_1(x_j)-\frac{1}{{\epsilon }^2}{\varphi }_1(x_j)-\lambda ({\varphi }_1(x_j){)}^3\right],j=1,2,\dots ,N-1,$$

$$E^n$$

$$\equiv E^0,n\ge 0,$$

$$\Vert e^n{\Vert }_{l^2}+\Vert {\delta }_x^{+}{e}^n{\Vert }_{l^2}\lesssim h^2+\frac{{\tau }^2}{{\epsilon }^6},0\le n\le \frac{T}{\tau },$$

$$u_j^{n+1}=\sum _{l=-N/2}^{N/2-1}{\widetilde{(u^{n+1})}}_l{e}^{i{\mu }_l(x_j-a)}=\sum _{l=-N/2}^{N/2-1}{\widetilde{(u^{n+1})}}_l{e}^{2ijl\pi /N},j=0,1,\dots ,N,n\ge 0,$$

$$\widetilde{(u^{n+1})}_{l}=\left\{\begin{array}[]{ll}\left[\cos(\omega_{l}^{0}\tau)+\displaystyle\frac{\alpha^{0}\left(1-\cos(\omega_{l}^{0}\tau)\right)}{(\varepsilon\omega_{l}^{0})^{2}}\right]\widetilde{(u^{0})}_{l}+\displaystyle\frac{\sin(\omega_{l}^{0}\tau)}{\omega_{l}^{0}}\widetilde{(\dot{u}^{0})}_{l}+\displaystyle\frac{\cos(\omega_{l}^{0}\tau)-1}{(\varepsilon\omega_{l}^{0})^{2}}\widetilde{(f^{0})}_{l},&n=0,\\ \\ -\widetilde{(u^{n-1})}_{l}+2\left[\cos(\omega_{l}^{n}\tau)+\displaystyle\frac{\alpha^{n}\left(1-\cos(\omega_{l}^{n}\tau)\right)}{(\varepsilon\omega_{l}^{n})^{2}}\right]\widetilde{(u^{n})}_{l}+\displaystyle\frac{2\left(\cos(\omega_{l}^{n}\tau)-1\right)}{(\varepsilon\omega_{l}^{n})^{2}}\widetilde{(f^{n})}_{l},&n\geq 1,\\ \end{array}\right.$$

$${\tilde{v}}_l=\frac{1}{N}\sum _{j=0}^{N-1}{v}_j{e}^{-i{\mu }_l(x_j-a)}=\frac{1}{N}\sum _{j=0}^{N-1}{v}_j{e}^{-2ijl\pi /N},l=-\frac{N}{2},-\frac{N}{2}+1,\dots ,\frac{N}{2}-1.$$

$${\dot{u}}_j^{n+1}=\sum _{l=-N/2}^{N/2-1}{\widetilde{({\dot{u}}^{n+1})}}_l{e}^{i{\mu }_l(x_j-a)}=\sum _{l=-N/2}^{N/2-1}{\widetilde{({\dot{u}}^{n+1})}}_l{e}^{2ijl\pi /N},j=0,1,\dots ,N,n\ge 0,$$

$$\widetilde{(\dot{u}^{n+1})}_{l}=\left\{\begin{array}[]{ll}-\omega_{l}\sin(\omega_{l}\tau)\widetilde{(u^{0})}_{l}+\cos(\omega_{l}\tau)\widetilde{(\dot{u}^{0})}_{l}-\displaystyle\frac{\sin(\omega_{l}\tau)}{\varepsilon^{2}\omega_{l}}\widetilde{(f^{0})}_{l},&n=0,\\ \\ \widetilde{(\dot{u}^{n-1})}_{l}-2\omega_{l}\sin(\omega_{l}\tau)\widetilde{(u^{n})}_{l}-2\displaystyle\frac{\sin(\omega_{l}\tau)}{\varepsilon^{2}\omega_{l}}\widetilde{(f^{n})}_{l},&n\geq 1.\\ \end{array}\right.$$

$$\Vert u(\cdot ,t_n)-I_N{u}^n{\Vert }_{L^2}\lesssim h^{m_0}+\frac{{\tau }^2}{{\epsilon }^4},\Vert {\partial }_x[u(\cdot ,t_n)-I_N{u}^n]{\Vert }_{L^2}\lesssim h^{m_0-1}+\frac{{\tau }^2}{{\epsilon }^4},0\le n\le \frac{T}{\tau },$$

$$\{\begin{array}{rl}& {\partial }_{t}u=0,\\ & {\partial }_{t}v+\frac{\lambda }{{\epsilon }^2}{u}^3=0,\end{array}\text{mbox}and\{\begin{array}{rl}& {\partial }_{t}u-v=0,\\ & {\partial }_{t}v-\frac{1}{{\epsilon }^2}{\partial }_{xx}u+\frac{1}{{\epsilon }^4}u=0.\end{array}$$

$$\begin{array}{rl}& v_j^{(1)}=v_j^n-\frac{\lambda \tau }{2{\epsilon }^2}(u_j^n{)}^3,\\ & v_j^{(2)}=\sum _{l=-N/2}^{N/2-1}\left[-{\omega }_l\sin \left({\omega }_l\tau \right){\widetilde{(u^n)}}_l+\cos \left({\omega }_l\tau \right){\widetilde{(v^{(1)})}}_l\right]e^{i{\mu }_l(x_j-a)},\\ & u_j^{n+1}=\sum _{l=-N/2}^{N/2-1}\left[\cos \left({\omega }_l\tau \right){\widetilde{(u^n)}}_l+\frac{\sin \left({\omega }_l\tau \right)}{{\omega }_l}{\widetilde{(v^{(1)})}}_l\right]e^{i{\mu }_l(x_j-a)},\\ & v_j^{n+1}=v_j^{(2)}-\frac{\lambda \tau }{2{\epsilon }^2}(u_j^{n+1}{)}^3,\end{array}j=0,1,\dots ,N,n\ge 0,$$

$$\Vert u(\cdot ,t_n)-I_N{u}^n{\Vert }_{L^2}\lesssim h^{m_0}+\frac{{\tau }^2}{{\epsilon }^2},\Vert {\partial }_x[u(\cdot ,t_n)-I_N{u}^n]{\Vert }_{L^2}\lesssim h^{m_0-1}+\frac{{\tau }^2}{{\epsilon }^2},0\le n\le \frac{T}{\tau }.$$

$$u_j^{n+1}=e^{i{t}_{n+1}/{\epsilon }^2}{z}_j^{n+1}+e^{-i{t}_{n+1}/{\epsilon }^2}\overline{z_j^{n+1}},j=0,1,\dots ,N,n\ge 0,$$

$$\begin{array}{rl}& z_j^{(1)}=\sum _{l=-N/2}^{N/2-1}{e}^{i{\mu }_l^2\tau /4}{\widetilde{(z^n)}}_l{e}^{i{\mu }_l(x_j-a)},\\ & z_j^{(2)}=e^{3i\lambda \tau |z_j^{(1)}{|}^2/2}{z}_j^{(1)},\\ & z_j^{n+1}=\sum _{l=-N/2}^{N/2-1}{e}^{i{\mu }_l^2\tau /4}{\widetilde{(z^{(2)})}}_l{e}^{i{\mu }_l(x_j-a)},\end{array}j=0,1,\dots ,N,n\ge 0.$$

$$u(x,t)=e^{it/{\epsilon }^2}z(x,t)+e^{-it/{\epsilon }^2}\overline{z}(x,t)+{\epsilon }^{2}w(x,t)+O({\epsilon }^4),x\in \Omega ,t\ge 0,$$

$$w(x,t)=$$

Full text

††journal: J. Comput. Phys.

1 Introduction

Consider the dimensionless nonlinear Klein-Gordon equation (NKGE) in $d$-dimensions ($d=1,2,3$) with cubic nonlinearity [10, 11, 30, 35, 49]:

[TABLE]

where $t$ is time, $\mathbf{x}\in\mathbb{R}^{d}$ is the spatial coordinate, $u:=u(\mathbf{x},t)$ is a complex-valued scalar field, $0<\varepsilon\leq 1$ is a dimensionless parameter inversely proportional to the speed of light, $\lambda\in\mathbb{R}$ is a given dimensionless parameter (positive and negative for defocusing and focusing self-interaction, respectively), and $\phi_{1}$ and $\phi_{2}$ are given complex-valued $\varepsilon$-independent initial data.

When $\lambda=0$, the above Klein-Gordon equation is known as the relativistic version of the Schrödinger equation for correctly describing the spinless relativistic composite particles, like the pion and the Higgs boson [35]. When $\lambda\neq 0$, the NKGE was widely adapted in plasma physics for modeling interaction between Langmuir and ion sound waves [21, 38] and in cosmology as a phonological model for dark-matter and/or black-hole evaporation [53, 75]. The NKGE (1.1) is time symmetric and conserves the energy

[TABLE]

For the derivation and nondimensionlization of (1.1), we refer to [11, 35, 61, 62] and references therein. For the well-posedness of the Cauchy problem (1.1) with a fixed $\varepsilon\in(0,1]$, e.g. $\varepsilon=1$, we refer to [35, 61, 62] and references therein. We remark here that when the initial data $\phi_{1}$ and $\phi_{2}$ are real-valued, then the solution $u$ of (1.1) is also real-valued; and this case has been widely studied analytically and numerically in the literature [1, 3, 22, 30, 36, 37, 41, 45, 46, 49, 54, 55, 56, 58, 59, 73, 72]. For simplicity of notations and without loss of generality, from now on, we assume that $\phi_{1}$ and $\phi_{2}$ are real-valued, and thus the solution $u$ of (1.1) is real-valued too. Our methods and results can be straightforwardly extended to the case when $\phi_{1}$ and $\phi_{2}$ are complex-valued and/or general nonlinearity in (1.1) [10, 11], and the conclusion will be remained the same.

When $\varepsilon\to 0^{+}$ in (1.1), due to that the energy $E(t)=O(\varepsilon^{-2})$ in (1) becomes unbounded, the analysis of the nonrelativistic limit of the solution $u$ becomes challenging and quite complicated. Fortunately, convergence from the NKGE (1.1) to a nonlinear Schrödinger equation has been extensively studied in the mathematics literature [10, 61, 62, 65]. Based on their results, the solution $u$ of the NKGE (1.1) propagates waves with wavelength at $O(\varepsilon^{2})$ and $O(1)$ in time and space, respectively, in the nonrelativistic limit regime, i.e. $0<\varepsilon\ll 1$. To illustrate this, Figure 1 shows the solution of (1.1) with $d=1$ and $\lambda=1$ and the initial data

[TABLE]

In fact, when $0<\varepsilon\ll 1$, formally by taking the ansatz [61, 62]

[TABLE]

where $z:=z(\mathbf{x},t)$ is a complex-valued function and $\overline{z}$ denotes the complex conjugate of $z$, the NKGE (1.1) can be formally reduced to – semi-limiting model – the nonlinear Schrödinger equation with wave operator (NLSW) under well-prepared initial data [4, 6]

[TABLE]

In addition, by dropping the small term $\varepsilon^{2}\partial_{tt}z$ in (1.5), one gets – limiting model – the nonlinear Schrödinger equation (NLSE) [4, 6, 61, 62]

[TABLE]

When $\varepsilon=1$ in (1.1), i.e. $O(1)$-wave speed regime, several numerical methods have been proposed and analyzed for the Cauchy problem (1.1) in the literature, see [22, 41, 60, 55, 59, 73] and references therein. Specifically, the finite difference time domain (FDTD) methods [41, 55, 72] have been demonstrated excellent performance in terms of efficiency and accuracy for (1.1) when $\varepsilon=1$. However, when $0<\varepsilon\ll 1$, i.e. in the nonrelativistic limit regime, it becomes much more challenging in designing and analyzing efficient and accurate numerical methods for (1.1) due to the highly oscillatory nature of the solution in time (cf. Fig. 1). To address this issue, Bao and Dong [11] established rigorous error bounds of the FDTD methods for (1.1), which depends explicitly on the mesh size $h$ and time step $\tau$ as well as the small parameter $\varepsilon\in(0,1]$. Based on their results [11], in order to obtain the ‘correct’ numerical solution of (1.1), the $\varepsilon$-resolution (or $\varepsilon$-scalability or meshing strategy) of the FDTD methods is $\tau=O(\varepsilon^{3})$ and $h=O(1)$, which is under-resolution in time with respect to $\varepsilon\in(0,1]$ regarding to the Shannon’s information theory [57, 68, 69] – to resolve a wave one needs a few points per wavelength – since the wavelength in time is at $O(\varepsilon^{2})$. To overcome the temporal under-resolution of the FDTD methods, they [11] proposed to adapt the exponential wave integrator (EWI) [52] for discretizing temporal derivatives in (1.1) and showed rigorously the $\varepsilon$-resolution of EWI is $\tau=O(\varepsilon^{2})$ and $h=O(1)$, which is optimal-resolution in time with respect to $\varepsilon\in(0,1]$. Later, the time-splitting (TS) method [63] was also applied to discretize the NKGE (1.1), and the method was shown as equivalent to one type EWI and thus it retains the same $\varepsilon$-resolution as that of the EWI but with an improved error bound regarding to the small parameter $\varepsilon\in(0,1]$ [40]. In fact, FDTD, EWI and TS methods perform very well when $\tau\to 0$ under $\varepsilon=\varepsilon_{0}$ fixed and they lose accuracy when $\varepsilon\to 0$ under $\tau=\tau_{0}$ fixed. At the meantime, Faou and Schratz [42] presented a class of limit integrators (LI) for (1.1) via solving numerically the limiting model NLSE (1.6) and obtained their error bounds. On the contrary, the LI methods perform very well when $\varepsilon\to 0$ under $\tau=\tau_{0}$ fixed and they lose accuracy when $\tau\to 0$ under $\varepsilon=\varepsilon_{0}$ fixed.

It is a natural question to ask on whether one can design a numerical method for the NKGE (1.1) such that it is uniformly accurate for $\varepsilon\in(0,1]$, i.e. super-resolution in time, especially in the nonrelativistic limit regime, since we have the solution structure (1.4) of the NKGE (1.1) via the limiting model NLSW (1.5) or NLSE (1.6). Recently, different uniformly accurate (UA) numerical methods have been designed and analyzed for the NKGE (1.1) including a multiscale time integrator (MTI) via a multiscale decomposition of the solution [10] and a two-scale formulation (TSF) method [26] as well as two uniformly and optimally accurate (UOA) methods [19, 20]. The main aim of this paper is to carry out a systematical comparison of different numerical methods which have been proposed for the NKGE (1.1) in terms of temporal/spatial accuracy and efficiency as well as $\varepsilon$-resolution for $\varepsilon\in(0,1]$, especially in the nonrelativistic limit regime.

The rest of the paper is organized as follows. The FDTD, EWI and TS methods as well as the LI schemes for the NKGE (1.1) are briefly reviewed in Section 2; the uniformly accurate methods are briefly reviewed in Section 3; and the uniformly and optimally accurate methods are briefly reviewed in Section 4. In Section 5, we present detailed comparison of different numerical methods; and in Section 6, we report convergence rates of the NKGE (1.1) to its limiting models NLSW (1.5) and NLSE (1.6) and show wave interactions of NKGE in two dimensions (2D). Finally, some conclusions are drawn in Section 7. Throughout this paper, we adopt the notation $A\lesssim B$ to represent that there exists a generic constant $C>0$, which is independent of $\tau$ (or $n$), $h$ and $\varepsilon$, such that $|A|\leq CB$.

2 Non-uniformly accurate numerical methods

In this section, we briefly review the FDTD, EWI and TS methods [11] as well as the LI methods [42] which have been proposed in the literature for discretizing the NKGE (2.1) (or (1.1)).

For simplicity of notation and without loss of generality, we only present the numerical methods in one dimension (1D). Generalization to high dimensions is straightforward by tensor product. As adapted in the literature [10, 20, 26], the NKGE (1.1) with $d=1$ is usually truncated onto a bounded interval $\Omega=(a,b)$ ($|a|$ and $b$ are usually taken large enough such that the truncation error is negligible) with periodic boundary condition

[TABLE]

Choose $\tau>0$ be the time step and $h=(b-a)/N$ be the mesh size with $N$ an even positive integer, denote the grid points as $x_{j}=a+jh$ for $j=0,1,\ldots,N$ and time steps as $t_{n}=n\tau$ for $n\geq 0$. Let $u_{j}^{n}$ be the numerical approximation of $u(x_{j},t_{n})$ for $0\leq j\leq N$ and $n\geq 0$ and $u^{n}=(u_{0}^{n},u_{1}^{n},\ldots,u_{N}^{n})^{T}$ be the solution vector at $t=t_{n}$, and define

[TABLE]

2.1 Finite difference time domain (FDTD) methods

Introduce the finite difference operators as

[TABLE]

As used in [11], the Crank-Nicolson finite difference (CNFD) method for discretizing (2.1) reads

[TABLE]

Similarly, the semi-implicit finite difference (SIFD) method is [11]

[TABLE]

and the leap-frog finite difference (LFFD) method is [11]

[TABLE]

The initial and boundary conditions in (2.1) are discretized as [11, 13]

[TABLE]

We remark here that we adapt (2.6) to compute the approximation at $t=t_{1}$ instead of the classical method

[TABLE]

i.e. replacing $\tau/\varepsilon^{2}$ by $\sin(\tau/\varepsilon^{2})$, such that the numerical solution $u^{1}$ is uniformly bounded for $\varepsilon\in(0,1]$.

As observed and stated in [11], the above CNFD, SIFD and LFFD methods are time symmetric and their memory cost is $O(N)$. The LFFD method is explicit and its computational cost per step is $O(N)$. It is conditionally stable and there is a severe stability condition which depends on both $h$ and $\varepsilon$, especially when $0<\varepsilon\ll 1$ [11]. In fact, it is the most efficient and accurate method among all FDTD methods for the NKGE when $\varepsilon=1$. The SIFD method is implicit, but it can be solved efficiently via the fast Fourier transform (FFT) and thus its computational cost per step is $O(N\,\ln N)$. It is conditionally stable and the stability condition depends on $\varepsilon$ and is independent of $h$ [11]. The CNFD method is implicit and at every time step a fully nonlinear coupled system needs to be solved. One main advantage is that it conserves the energy (1) in the discrete level as [11]

[TABLE]

which immediately implies that it is unconditionally stable when $\lambda\geq 0$. In addition, under proper regularity of the solution $u$ of the NKGE (2.1) and stability conditions for the SIFD and LFFD methods, the following rigorous error bound was established for the three FDTD methods [11]

[TABLE]

where $T>0$ is a fixed time and the error function $e^{n}$ is defined as $e_{j}^{n}=u(x_{j},t_{n})-u_{j}^{n}$ for $0\leq j\leq N$ and $n\geq 0$. This error bound suggests that the FDTD methods are second order in both space and time discretization for any fixed $\varepsilon=\varepsilon_{0}$ and the $\varepsilon$-resolution of the FDTD methods is $h=O(1)$ and $\tau=O(\varepsilon^{3})$ in the nonrelativistic limit regime, i.e. $0<\varepsilon\ll 1$, which immediately show that the temporal resolution with respect to $\varepsilon\in(0,1]$ of the FDTD methods is under-resolution in time since the wavelength in time is at $O(\varepsilon^{2})$.

2.2 Exponential wave integrator (EWI)

As it has been proposed in [11], the NKGE (2.1) is discretized in space by the Fourier (pseudo)spectral method and followed by adapting an exponential wave integrator (EWI) in time which has been widely used for discretizing second order oscillatory differential equations in the literature [31, 39, 44, 47, 48, 50, 51, 52].

Let $u_{j}^{n}$ and $\dot{u}_{j}^{n}$ be the approximations of $u(x_{j},t_{n})$ and $\partial_{t}u(x_{j},t_{n})$, respectively, for $0\leq j\leq N$ and $n\geq 0$ and take $u_{j}^{0}=\phi_{1}(x_{j}),\dot{u}_{j}^{0}=\phi_{2}(x_{j})/\varepsilon^{2}$ for $j=0,1,\ldots,N$. When the EWI is taken as the Gautschi’s quadrature [11, 44, 47, 48, 51], a Gautschi-type exponential wave integrator Fourier pseudospectral (EWI-FP) method [11] reads as:

[TABLE]

where

[TABLE]

with $f(u)=\lambda u^{3}$, $f^{n}:=f(u^{n})=(f(u_{0}^{n}),f(u_{1}^{n}),\ldots,f(u_{N}^{n}))^{T}$, $\mu_{l}=2l\pi/(b-a)$, $\omega_{l}^{n}=\frac{1}{\varepsilon^{2}}\sqrt{1+\varepsilon^{2}(\mu_{l}^{2}+\alpha^{n})}$ ($l=-N/2,\ldots,N/2-1$) with $\alpha^{n}=\max\left\{\alpha^{n-1},\;\max_{0\leq j\leq N}\{\lambda(u_{j}^{n})^{2}\}\right\}$ for $n\geq 0$ and $\alpha^{-1}=0$ being the stabilization constants, and $\widetilde{v}_{l}$ ($-N/2\leq l\leq N/2-1$) being the discrete Fourier transform coefficients of the vector $v=(v_{0},v_{1},\ldots,v_{N})^{T}$ with $v_{0}=v_{N}$ defined as

[TABLE]

Of course, in practice if the approximation of the first order derivative in time is needed, then they can be obtained as [11]

[TABLE]

where

[TABLE]

As it can be seen, EWI-FP is explicit and time symmetric. The memory cost is $O(N)$ and computational cost per step is $O(N\,\ln N)$. In addition, the EWI-FP is unconditionally stable due to the stabilization constant $\alpha^{n}$ [11]. Under proper regularity of the solution $u$ of the NKGE (2.1) and the assumption $\tau\lesssim\varepsilon^{2}$, the following rigorous error bound was established for the EWI-FP method [11]

[TABLE]

where $m_{0}\geq 2$ depends on the regularity of the solution $u$ of (2.1) and $I_{N}$ is the standard interpolation operator [71]. The error bounds suggest that EWI-FP is spectral order in space if the solution is smooth and is second order in time for any fixed $\varepsilon=\varepsilon_{0}$ and the $\varepsilon$-resolution is $h=O(1)$ and $\tau=O(\varepsilon^{2})$ in the nonrelativistic limit regime, i.e. $0<\varepsilon\ll 1$, which immediately show that EWI-FP is optimal-resolution in time with respect to $\varepsilon\in(0,1]$ since the wavelength in time is at $O(\varepsilon^{2})$. Recently, the EWI-FP method has been extended to arbitrary even order in time [60, 74].

2.3 Time-splitting (TS) method

The time-splitting method [63] has been widely used to solve different (partial) differential equations and it has shown great advantages in many cases, such as for the (nonlinear) Schrödinger equation [5, 63]. As proposed in [40], in order to adapt the TS method for solving the NKGE, the NKGE (2.1) is first re-formulated into a first order system by introducing $v:=v(x,t)=\partial_{t}u(x,t)$. Then the first order system is split into

[TABLE]

Let $u_{j}^{n}$ and $v_{j}^{n}$ be the approximations of $u(x_{j},t_{n})$ and $v(x_{j},t_{n})$, respectively, for $0\leq j\leq N$ and $n\geq 0$ and take $u_{j}^{0}=\phi_{1}(x_{j})$, $v_{j}^{0}=\phi_{2}(x_{j})/\varepsilon^{2}$ for $j=0,1,\ldots,N$. Then a second-order time splitting Fourier pseudospectral method (TS-FP) [40] reads as:

[TABLE]

where $\omega_{l}=\frac{1}{\varepsilon^{2}}\sqrt{1+\varepsilon^{2}\mu_{l}^{2}}$ for $l=-N/2,\ldots,N/2-1$.

Again, the TS-FP (2.14) is explicit and time symmetric. Its memory cost is $O(N)$ and computational cost per time step is $O(N\,\ln N)$. We remark here that the TS-FP (2.14) is mathematically equivalent to an EWI via trapezoidal quadrature (or known as Deuflhard-type exponential integrator [39]) for solving the NKGE (2.1) (or (1.1)) [40, 50]. Under the condition $\tau\lesssim\varepsilon^{2}$, the following error bound was observed for the TS-FP in [40]:

[TABLE]

The above error bound could be rigorously obtained by the super-convergence analysis in [8, 9, 29]. It can be seen that the error bound (2.15) is an improved error bound compared to the error bound (2.13) regarding to the small parameter $\varepsilon$ when $0<\tau\lesssim\varepsilon^{2}$ and $0<\varepsilon\ll 1$ (cf. Tabs. 9&10). Of course, due to the convergence restriction $\tau\lesssim\varepsilon^{2}$, the $\varepsilon$-resolution of TS-FP is still $h=O(1)$ and $\tau=O(\varepsilon^{2})$ in the nonrelativistic limit regime. Thus the TS-FP is also optimal-resolution in time with respect to $\varepsilon\in(0,1]$.

2.4 Limit integrators (LIs)

As presented in [42], a class of limit integrators (LIs) has been designed for the NKGE (2.1) (or (1.1)) with different order of accuracy in terms of $\varepsilon$ when $0<\varepsilon\ll 1$. In the LIs, the limiting equation of the NKGE (1.1), e.g. (1.6), is solved numerically and the numerical solution of the NKGE (1.1) is constructed via the ansatz (1.4).

In practice, the NLSE (1.6) in 1D is truncated on a bounded computational domain $\Omega=(a,b)$ with periodic boundary condition and then it is discretized by the second order time-splitting Fourier pseudospectral (TSFP) method [2, 5, 14, 15]. Let $u_{j}^{n}$ and $z_{j}^{n}$ be the approximations of $u(x_{j},t_{n})$ and $z(x_{j},t_{n})$, respectively, for $0\leq j\leq N$ and $n\geq 0$ and take $u_{j}^{0}=\phi_{1}(x_{j})$, $z_{j}^{0}=z_{0}(x_{j})$ for $j=0,1,\ldots,N$. Then a first order (with respect to the small parameter $\varepsilon$) limit integrator Fourier pseudospectral (LI-FP1) method was proposed in [42] as:

[TABLE]

where $z^{n+1}$ is a numerical approximation of (1.6) by a TSFP method [2, 5, 14, 15] and is given as

[TABLE]

Again, as presented in [42], when $0<\varepsilon\ll 1$, formally by taking the following ansatz (found by the modulated Fourier expansion [31, 32, 43, 50]) which is more accurate than (1.4) for approximating the solution of NKGE (2.1),

[TABLE]

one can obtain $z:=z(x,t)$ still satisfies the NLSE (1.6) and $w:=w(x,t)$ is given by [42]

[TABLE]

with $v:=v(x,t)$ satisfying [42]

[TABLE]

and the initial condition [42]

[TABLE]

The NLSE (1.6) can be solved by the second order TSFP method [2, 5, 14, 15] (cf. (2.17) for the case of 1D) as before. In order to solve (2.20) numerically [42], it is split into a kinetic part

[TABLE]

and a potential part

[TABLE]

and then the flow is composed by a second order splitting scheme as $\Phi(\tau)\approx\Phi_{k}\left(\frac{\tau}{2}\right)\Phi_{p}\left(\tau\right)\Phi_{k}\left(\frac{\tau}{2}\right).$ The kinetic part can be integrated as usual, while the potential part is integrated in its vector form by an exponential trapezoidal rule in [42].

Here we present the method in 1D and truncate (2.22) on the bounded domain $\Omega=(a,b)$ with periodic boundary condition. Let $u_{j}^{n}$, $z_{j}^{n}$, $w_{j}^{n}$ and $v_{j}^{n}$ be the approximations of $u(x_{j},t_{n})$, $z(x_{j},t_{n})$, $w(x_{j},t_{n})$ and $v(x_{j},t_{n})$, respectively, for $0\leq j\leq N$ and $n\geq 0$ and take $u_{j}^{0}=\phi_{1}(x_{j})$, $z_{j}^{0}=z_{0}(x_{j})$, $v_{j}^{0}=v_{0}(x_{j})\approx-\frac{\lambda}{2}(z_{j}^{0})^{3}+\frac{\lambda}{4}(\overline{z_{j}^{0}})^{3}+\frac{3\lambda}{2}|z_{j}^{0}|^{2}\overline{z_{j}^{0}}+\frac{1}{2}\left(\partial_{xx}^{\mathcal{F}}z^{0}_{j}-\partial_{xx}^{\mathcal{F}}\overline{z^{0}_{j}}\right)$ for $j=0,1,\ldots,N$, where $\partial_{xx}^{\mathcal{F}}$ is the standard Fourier pseudospectral approximation of the operator $\partial_{xx}$ on the bounded domain $\Omega=(a,b)$, e.g. $\partial_{xx}^{\mathcal{F}}z_{j}^{0}=\partial_{xx}(I_{N}z^{0})(x_{j})$ [71]. Then a second order (with respect to the small parameter $\varepsilon\in(0,1]$) limit integrator Fourier pseudospectral (LI-FP2) method is given [42] as

[TABLE]

where $z_{j}^{n+1}$ ($j=0,1,\ldots,N$) are given in (2.17) and $w^{n+1}$ is an approximation of (2.4) as

[TABLE]

Here $v^{n+1}$ is a numerical solution of (2.20) by a TSFP method [42] and is given as:

[TABLE]

where

[TABLE]

with

[TABLE]

Here $Re(f)$ and $Im(f)$ denote the real and imaginary parts of $f$, respectively, and $\partial_{xxxx}^{\mathcal{F}}$ is the Fourier pseudospectral approximation of $\partial_{xxxx}$ on the bounded domain $\Omega=(a,b)$ [71].

As stated and proved in [42], both LI-FP1 and LI-FP2 are explicit, unconditionally stable and time symmetric, and their memory cost is $O(N)$ and computational cost per step is $O(N\,\ln N)$. In addition, under proper regularity of the solution $u$ of the NKGE (2.1) and $z$ of the NLSE (1.6), the following rigorous error bound was established for the LI-FP1 method [42]

[TABLE]

where $m_{1}\geq 1$ depends on the regularity of the solution $u$ of (2.1). Similarly, the following rigorous error bound was established for the LI-FP2 method [42]

[TABLE]

These error bounds suggest that both LI-FP1 and LI-FP2 methods are spectral order in space if the solution is smooth and when $\varepsilon\to 0^{+}$, and LI-FP1 and LI-FP2 are second order in time when $0<\varepsilon\lesssim\tau$ and $0<\varepsilon\lesssim\tau^{1/2}$, respectively. The $\varepsilon$-resolution of the two methods is $h=O(1)$ and $\tau=O(1)$ in the nonrelativistic limit regime, i.e. $0<\varepsilon\ll 1$, which immediately show that both LI-FP1 and LI-FP2 are super-resolution in time with respect to $0<\varepsilon\ll 1$ since the wavelength in time is at $O(\varepsilon^{2})$. On the contrary, when $\varepsilon=\varepsilon_{0}$ is fixed, e.g. $\varepsilon=1$, there is no convergence of LI-FP1 and LI-FP2 for the NKGE (2.1) (or (1.1)).

3 Uniformly accurate (UA) methods

In this section, we review the uniformly accurate MTI [13, 10] and TSF method [26] which have been proposed in the literature for discretizing the NKGE (2.1) (or (1.1)).

3.1 A multiscale time integrator (MTI)

As proposed in [10], the MTI was designed via a multiscale decomposition of the solution of the NKGE (1.1) and adapting the EWI-FP method for discretizing the decomposed sub-problems.

For any fixed $n\geq 0$, by assuming that the initial data at $t=t_{n}$ is given as

[TABLE]

and decomposing the solution $u(\mathbf{x},t)=u(\mathbf{x},t_{n}+s)$ of the NKGE (1.1) on the time interval $t\in[t_{n},t_{n+1}]$ as [10, 61, 62]

[TABLE]

then a multiscale decomposition by the $\varepsilon$-frequency (MDF) of the NKGE (1.1) can be given as [10, 12]

[TABLE]

with the well-prepared initial data for $z^{n}$ and small initial data for $r^{n}$ as [10, 12]

[TABLE]

where

[TABLE]

Then the problem (3.3) with (3.4) is truncated on a bounded domain with periodic boundary condition and then discretized by the EWI-FP method [10] with details omitted here for brevity. After solving numerically the decomposed problem (3.3) with (3.4), the solution of the NKGE (1.1) at $t=t_{n+1}$ is reconstructed by the ansatz (3.2) by setting $s=\tau$ [10].

For the convenience of the reader and simplicity of notation, here we only present the method in 1D on $\Omega=(a,b)$ with the periodic boundary condition. Let $u_{j}^{n}$ and $\dot{u}_{j}^{n}$ be the approximations of $u(x_{j},t_{n})$ and $\partial_{t}u(x_{j},t_{n})$, respectively; and let $z_{j}^{n+1},$ $\dot{z}_{j}^{n+1},$ $r_{j}^{n+1}$ and $\dot{r}_{j}^{n+1}$ be the approximations of $z^{n}(x_{j},\tau),$ $\partial_{s}z^{n}(x_{j},\tau),$ $r^{n}(x_{j},\tau)$ and $\partial_{s}r^{n}(x_{j},\tau)$, respectively, for $j=0,1,\ldots,N$ and $n\geq 0$. Choosing $u_{j}^{0}=\phi_{1}(x_{j})$ and $\dot{u}_{j}^{0}=\phi_{2}(x_{j})/\varepsilon^{2}$, then a multiscale time integrator Fourier pseudospectral (MTI-FP) method for discretizing the NKGE (2.1) is given as [10]

[TABLE]

where $z^{n+1},$ $\dot{z}^{n+1}$, $r^{n+1}$ and $\dot{r}^{n+1}$ are numerical approximations of (3.3) with (3.4) as

[TABLE]

with

[TABLE]

and

[TABLE]

Here we adopt the following functions [10]

[TABLE]

The MTI-FP is explicit, unconditionally stable, and its memory cost is $O(N)$ and computational cost per step is $O(N\,\ln N)$. As established in [10], under proper regularity of the solution $u$ of the NKGE (2.1), the following two error bounds were established by using two different techniques for the MTI-FP method [10]

[TABLE]

which imply a uniform error bound for $\varepsilon\in(0,1]$ [10]

[TABLE]

where $m_{2}\geq 1$ depends on the regularity of the solution $u$ of (2.1).

These error bounds suggest that the MTI-FP method is uniformly spectral order in space if the solution is smooth and is uniformly first order in time for $0<\varepsilon\leq 1$. The $\varepsilon$-resolution is $h=O(1)$ and $\tau=O(1)$ in the nonrelativistic limit regime, which immediately show that the MTI-FP is super-resolution in time with respect to $\varepsilon\in(0,1]$ since the time step can be chosen as independently of $\varepsilon$ although the solution is highly oscillatory in time at wavelength $O(\varepsilon^{2})$.

3.2 Two-scale formulation (TSF) method

As presented in [26], the TSF method was constructed by separating the slow time scale $t$ and the fast time scale $\xi=t/\varepsilon^{2}$ and re-formulating the NKGE (1.1) into a two-scale formulation. This approach offers a general strategy to design uniformly accurate schemes for highly oscillatory differential equations and PDEs which contain general nonlinearity or strong couplings [27, 33, 34].

Introducing [26]

[TABLE]

then the NKGE (1.1) can be re-written as a first order PDE

[TABLE]

Let

[TABLE]

so as to filter out the main oscillation in the above PDE. Then one gets [26]

[TABLE]

Introduce $U:=U(\mathbf{x},t,\xi)$ with $\xi$ interpreted as another ‘space’ variable on torus $\mathbb{T}=\mathbb{R}/(2\pi\mathbb{Z})$ such that

[TABLE]

with $t$ the slow time variable and $\xi=t/\varepsilon^{2}$ the fast time variable [26]. Noticing (3.14), one needs to request $U$ satisfies the following PDE [26]

[TABLE]

where $U_{0}(\mathbf{x},\xi)$ to be determined later satisfying $U_{0}(\mathbf{x},0)=v_{0}(\mathbf{x})$ and

[TABLE]

with $\phi:=\phi(\mathbf{x})$.

The initial data $U_{0}(\mathbf{x},\xi)$ in (3.16) is only prescribed at one point, i.e. $\xi=0$, so there is some freedom to choose the initial data in order to bound the time derivatives of $U$. By using the Chapman-Enskog expansion, the initial data $U_{0}(\mathbf{x},\xi)$ was obtained at different order of accuracy in term of $\varepsilon$ [26]. For example, the initial data at first order of accuracy was given as [26]

[TABLE]

such that $\partial_{t}^{2}U(x,t,\xi)=O(1)$ for fixed $t\geq 0$ as $\varepsilon\to 0$ [26], where

[TABLE]

with the operators ${\mathcal{L}}$ and $\Pi$ defined as

[TABLE]

for some periodic function $\varphi:=\varphi(\xi)$ on $\mathbb{T}$.

Let $U^{n}(\mathbf{x},\xi)$ be the numerical approximation of $U(\mathbf{x},t_{n},\xi)$ for $n\geq 0$. Then (3.16) with (3.18) can be discretized in time as

[TABLE]

with $U^{n}(\mathbf{x},\xi)=U_{0}^{\rm 1st}(\mathbf{x},\xi)$. In the $\xi$-direction, thanks to the periodicity, one can further discretize (3.19) by the Fourier pseudospectral method as: Let $h_{\xi}=2\pi/N_{\xi}$ with $N_{\xi}$ an even positive integer, $\xi_{m}=mh_{\xi}$ and $U^{n}_{m}(\mathbf{x})$ be the approximation of $U^{n}(\mathbf{x},\xi_{m})$ for $m=0,1,\ldots,N_{\xi}$, denote $U^{n}(\mathbf{x})=(U_{0}^{n}(\mathbf{x}),U^{n}_{1}(\mathbf{x}),\ldots,$ $U^{n}_{N_{\xi}}(\mathbf{x}))^{T}$, take $U^{0}_{m}(\mathbf{x})=U_{0}^{\rm 1st}(\mathbf{x},\xi_{m})$ for $m=0,1,\ldots,N_{\xi}$, then one can get

[TABLE]

with $\widetilde{U}^{n}_{l}(\mathbf{x})=\sum_{m=0}^{N_{\xi}-1}U^{n}_{m}(\mathbf{x})\,e^{-il\xi_{m}}$ and

[TABLE]

Then (3.20) with (3.21) will be first truncated (in $\mathbf{x}$) on a bounded computational domain with periodic boundary condition and then discretized by the standard Fourier pseudospectral method with details omitted here for brevity [26]. Finally, noticing (3.11), (3.13) and (3.15), one can reconstruct the approximation of the solution $u$ of the NKGE (1.1) (or (2.1)). For the simplicity of notations, here we only present a first order two-scale formulation Fourier pseudospectral (TSF-FP1) method in 1D as [26]:

[TABLE]

where $u_{j}^{0}=\phi_{1}(x_{j})$ for $j=0,1,\ldots,N$.

Similarly, by taking the initial data as

[TABLE]

such that $\partial_{t}^{3}U(x,t,\xi)=O(1)$ for fixed $t\geq 0$ as $\varepsilon\to 0$ [26], where

[TABLE]

Then (3.16) with (3.23) can be discretized by a second order scheme in time as

[TABLE]

with $U^{n}(\mathbf{x},\xi)=U_{0}^{\rm 2nd}(\mathbf{x},\xi)$. Similarly, (3.24) can be discretized in $\xi$-direction via the Fourier pseudospectral method, truncated in $\mathbf{x}$-direction onto a bounded computational domain with periodic boundary condition and then discretized via the Fourier pseudospecral method with details omitted here for brevity [26]. Finally one can obtain a second order two-scale formulation Fourier pseudospectral (TSF-FP2) for the NKGE (1.1) (or (2.1)) via the reconstruction (3.22).

As shown in [26], both TSF-FP1 and TSF-FP2 are explicit, unconditionally stable, and its memory cost is $O(N_{\xi}N)$ and computational cost per step is $O(N_{\xi}N\,\ln(N_{\xi}N))$. Under proper regularity of the solution $U$ of the PDE (3.16), the following error bound was established for TSF-FP1 [26]

[TABLE]

and respectively, for TSF-FP2 [26]

[TABLE]

where $m_{0}$ and $m_{1}$ are two positive integers which depend on the regularity of solution $U$ of (3.16) in $\mathbf{x}$- and $\xi$-direction, respectively.

These error bounds suggest that both TSF-FP1 and TSF-FP2 methods are uniformly spectral order in space and in $\xi$-direction if the solution is smooth, and the TSF-FP1 and TSF-FP2 are uniformly first and second order, respectively, in time. Again, the $\varepsilon$-resolution is $h=O(1)$ and $\tau=O(1)$ in the nonrelativistic limit regime, which immediately show that both TSF-FP1 and TSF-FP2 are super-resolution in time with respect to $\varepsilon\in(0,1]$ since the time step can be chosen independently on $\varepsilon$ although the solution is highly oscillatory in time at wavelength $O(\varepsilon^{2})$. We remark that the finite difference integrator (3.19) or (3.24) from [26] is not the unique choice for discritizating the two-scale system (3.16). The formulation (3.16) and the well-prepared initial data (3.18) or (3.23) are essential for the TSF approach, and other numerical discretizations such as EWI-FP can also be applied to solve (3.16), see [27, 34, 76].

4 Uniformly and optimally accurate (UOA) methods

In this section, we review briefly two UA methods with optimal convergence rate in time and/or computational costs, i.e. uniformly second-order in time without solving a problem in one more spatial dimension. One is the iterative exponential-type integrator in [20], and the other is a MTI based on higher order multiscale expansion by frequency in [19].

4.1 An iterative exponential integrator (IEI)

Without using higher order approximations or extra dimensions, a second order UOA method for the NKGE (1.1) (or (2.1)) was very recently proposed in [20] by using an iterative exponential integrator. By reformulating the NKGE (1.1) into the first order PDE (3.12) and then introducing

[TABLE]

one finds

[TABLE]

which based on the Duhammel’s formula gives

[TABLE]

Here $D_{\varepsilon}$ is defined the same as in (3.17). As used in [20], an exponential integrator is proposed by plugging (4.2) iteratively into the cubic terms (see this technique also in [66]). To describe the scheme, the following functions and operators are introduced.

[TABLE]

For simplicity of notations, here we only present the method in 1D. In this case, $\Delta=\partial_{xx}$. Let $v_{*}^{n}$ be an approximation of $v_{*}(x,t_{n})$ for $n\geq 0$ via time integration, and take $u^{0}(x)=\phi_{1}(x)$ and $v_{*}^{0}=\phi_{1}(x)-i(1-\varepsilon^{2}\partial_{xx})^{-1/2}\phi_{2}(x)$. Then an iterative exponential integrator (IEI) scheme for approximating (4.2) reads:

[TABLE]

where

[TABLE]

with

[TABLE]

Combining (3.11) and (4.1), a semi-discretized approximation of the NKGE (2.1) is given as

[TABLE]

In practice, performing all the differential operations in the above IEI method by the Fourier pseudospectral approximation with details omitted here for brevity [66], we obtain the iterative exponential integrator Fourier pseudospectral (IEI-FP) scheme.

The IEI-FP is explicit, unconditionally stable, and its memory cost is $O(N)$ and computational cost per step is $O(N\,\ln N)$. As established in [66], under proper regularity of the solution $u$ of the NKGE (2.1) [66] (which is weaker than that is needed for MTI-FP or TSF-FP), the following error bound was established for IEI-FP [66]

[TABLE]

4.2 A higher order MTI

Another uniform second-order in time UOA method for solving the NKGE (1.1) was very recently presented in [19] via a higher order multiscale expansion of the solution of the NKGE. As obtained in [19], the solution $u(\mathbf{x},t)$ of (1.1) is expanded as

[TABLE]

where $c.c.$ represents the complex conjugate of the whole expression before it within the bracket, and $v:=v(\mathbf{x},t)$ solves the following nonlinear Schrödinger equation with wave operator (NLSW) under well-prepared initial data [6, 10, 19]

[TABLE]

where

[TABLE]

and $R:=R(\mathbf{x},t)$ solves the following NKGE with small initial data [19]

[TABLE]

with

[TABLE]

It was shown that $R(\mathbf{x},t)=O(\varepsilon^{2})$ [19] and thus the multiscale expansion (4.5) without the last term is a higher order multiscale expansion of the solution of the NKGE (1.1) at $O(\varepsilon^{4})$ [19]. Then the decomposed problems (4.6) and (4.7) are solved numerically by the EWI-FP method [19].

Again, for the convenience of the reader and simplicity of notations, here we only present the method in 1D on $\Omega=(a,b)$ with the periodic boundary condition. Let $u_{j}^{n}$ and $\dot{u}_{j}^{n}$ be the approximations of $u(x_{j},t_{n})$ and $\partial_{t}u(x_{j},t_{n})$, respectively; and let $v_{j}^{n},$ $\dot{v}_{j}^{n},$ $R_{j}^{n}$ and $\dot{R}_{j}^{n}$ be the approximations of $v(x_{j},t_{n}),$ $\partial_{t}v(x_{j},t_{n}),$ $R(x_{j},t_{n})$ and $\partial_{t}R(x_{j},t_{n})$ respectively, for $j=0,1,\ldots,N$ and $n\geq 0$. Choosing $u_{j}^{0}=\phi_{1}(x_{j})$ and $\dot{u}_{j}^{0}=\phi_{2}(x_{j})/\varepsilon^{2}$, $v^{0}_{j}=v_{0}(x_{j})$, $\dot{v}^{0}_{j}=v_{1}(x_{j})$, $R^{0}_{j}=R_{0}(x_{j})$ and $\dot{R}^{0}_{j}=R_{1}(x_{j})$ for $j=0,1\ldots,N$, then a high-order multiscale time integrator Fourier pseudospectral (MTI-FP2) method for discretizing the NKGE (2.1) is given as [19]

[TABLE]

Here $v^{n+1}$ and $\dot{v}^{n+1}$ are approximations of the NLSW (4.6) by an EWI-FP [6, 10, 19] as

[TABLE]

where

[TABLE]

[TABLE]

with $g^{n}=(g_{0}^{n},g_{1}^{n},\ldots,g_{N}^{n})^{T}$ given as

[TABLE]

Similarly, here $R^{n+1}$ and $\dot{R}^{n+1}$ are approximations of the NKGE (4.7) by an EWI-FP [10, 11, 19] as

[TABLE]

where

[TABLE]

with

[TABLE]

and

[TABLE]

Here we adopt the following functions from [19]: for $n\geq 0,\,l=-N/2,\ldots,N/2-1$ and $k=3,5$,

[TABLE]

Again, the MTI-FP2 is explicit, unconditionally stable, and its memory cost is $O(N)$ and computational cost per step is $O(N\,\ln N)$. As established in [19], under proper regularity of the solution $u$ of the NKGE (2.1), the following error bound was established for MTI-FP2 [19]

[TABLE]

We remark here that two new UOA schemes named as micro-macro method and pull-back method were proposed very recently in [28] based on the averaging theory. The accuracy of the micro-macro method is very similar to TSF-FP2 and the efficiency is very similar to MTI-FP2 and IEI-FP. The pull-back is implicit but with superior long-time behaviour over other methods [28].

5 Numerical comparisons and results

In this section, we report the performance of different numerical methods reviewed in previous sections and carry out a systematical comparison.

In order to do so, we take $d=1$ and $\lambda=1$ in (1.1) and choose the initial data as

[TABLE]

The problem is solved numerically on a bounded computational domain $\Omega=(-16,16)$. The ‘exact’ solution of the NKGE (2.1) is obtained numerically by TSF-FP2 with a very small step size, i.e. $\tau=10^{-6},h=1/64,h_{\xi}=\pi/64$. Define error

[TABLE]

where $P_{N}$ is the standard projection operator [71]. We depict the errors at $t=1$. The temporal error and spatial error of each numerical method are studied and shown separately in the following.

5.1 Spatial errors

We first test and compare the spatial discretization error of different numerical methods. For spatial error analysis, the time step $\tau$ is chosen small enough such that the discretization error in time is negligible, e.g. $\tau=10^{-6}$.

The three finite difference methods share almost the same discretization error in space. Thus, we only give the spatial error of ECFD in Tab. 1 as a representative and omit the results of SIFD and LFFD. The errors of LI-FP1 and LI-FP2 are given in Tab. 2 and Tab. 3, respectively. These errors are the spatial errors mixed with the model reduction errors. The spatial error of MTI-FP is given in Tab. 4. The results of EWI-FP, TS-FP, TSF-FP2, IEI-FP and MTI-FP2 behave similarly as that of MTI-FP since they share the same Fourier discretization, so the corresponding results have been omitted for brevity as well. The error of TSF-FP1 in the extra space direction $\xi$ is given in Tab. 5, which represents the very similar corresponding results of TSF-FP2.

From Tabs. 1-5, we can draw the following observations:

(i) ECFD, SIFD and LFFD have second order accuracy in space error, while EWI-FP, TS-FP, MTI-FP, TSF-FP1 and TSF-FP2 have spectral accuracy in space. The errors are uniform in space in terms of $\varepsilon$ with spatial $\varepsilon$-scalability $h=O(1)$. Thus, when the initial data of NKGE is smooth enough, the Fourier pseudospectral discretization in space is obviously more efficient than finite difference discretization.

(ii) The spatial errors of the LI-FP1 and LI-FP2 are mixed with the residue of $O(\varepsilon^{2})$ and $O(\varepsilon^{4})$, respectively from the model reductions. Thus, for a fixed $0<\varepsilon\leq 1$, the spectral accuracy of space discretization is broken and the error is bounded from blow.

(iii) The errors of TSF-FP1 and TSF-FP2 in the extra space direction $\xi$ are of spectral accuracy. The error is uniformly bounded for all $\varepsilon\in(0,1]$ which allows the use of $h_{\xi}=O(1)$. The computational resource needed in the $\xi$-direction is not very heavy, i.e. $N_{\xi}=32$ is enough to get machine accuracy for all the $\varepsilon$ in this example. The smaller $\varepsilon$ is, the less grid points are needed in $\xi$ direction to reach the machine accuracy.

5.2 Temporal errors

For temporal error analysis, the mesh size $h$ (and so is $h_{\xi}$ for TSF-FP) is chosen small enough such that the discretization error in space is negligible. The detailed data of the used $h,h_{\xi}$ is given case by case below.

The results of ECFD, SIFD, LFFD, EWI-FP, TS-FP, LI-FP1, LI-FP2, MTI-FP, TSF-FP1, TSF-FP2, IEI-FP and MTI-FP2 are shown in Tabs. 6-17, respectively. For the LFFD method, in order to show the temporal discretization error but meanwhile to satisfy the stability condition, we choose [8, 9]

[TABLE]

in Tab. 8 for the temporal error. To illustrate the UA property of MTI-FP, TSF-FP1 and TSF-FP2, we also define the error

[TABLE]

The convergence rate of each method is shown along with the error. In the tables, we highlight the error of each classical method in the table when the time step is chosen according to its $\varepsilon$-scalability. The efforts made here are to illustrate the convergence rate of each method and the $\varepsilon$-scalability in the limit regime.

As a summary, a detailed table on the comparison of computational complexity of each method has been given in Tab. 18. The comparison of the temporal error of each method in the classical regime, i.e. $\varepsilon=O(1)$, has been given in Tab. 19 together with the computational time. The comparisons of the temporal error in the limit regime are given in Tab. 20. Tab. 21 shows the temporal error of different methods under the natural mesh strategy, i.e. $\tau=O(\varepsilon^{2})$ which fully resolves the temporal wavelength of the oscillation. All methods are programmed with Matlab and run on an Intel i3-3120M 2.5GHz CPU laptop.

5.3 Comparison of different methods

From Tabs. 6-21, we can draw the following conclusions:

(i) All FDTD methods have temporal $\varepsilon$-scalability $\tau=O(\varepsilon^{3})$ (cf. Tabs. 6-8). In the classical regime, LFFD is the most accurate and efficient method among the three. However, all FDTD methods become inefficient when $\varepsilon$ becomes small.

(ii) Both EWI-FP and TS-FP have temporal $\varepsilon$-scalability $\tau=O(\varepsilon^{2})$ (cf. Tabs. 9&10). While when $\tau\lesssim\varepsilon^{2}$, TS-FP has an improved error bound at $\tau^{2}/\varepsilon^{2}$ with respect to the small parameter $\varepsilon\in(0,1]$. Both methods perform very well in the classical regime (cf. Tab. 19), but they are unsatisfactory in the limit regime. They have very similar efficiency but TS-FP is more accurate when $\varepsilon$ is small.

(iii) The two LI methods are accurate in the limit regime, but both of them do not have convergence in the classical or the intermediate regime (cf. Tabs. 11&12). LI-FP2 is more accurate than LI-FP1 due to the correction, but LI-FP2 is much less efficient.

(iv) MTI-FP, TSF-FP1, TSF-FP2, IEI-FP and MTI-FP2 have temporal $\varepsilon$-scalability $\tau=O(1)$ and they offer uniformly correct results for all $\varepsilon\in(0,1]$ (cf. Tabs. 13-15). All the five UA methods have some temporal convergence order reductions in the resonance regime. Between the two first order UA schemes, MTI-FP is more accurate than TSF-FP1 (cf. Tabs. 13&14). TSF-FP2 is the most accurate method among the five (cf. Tabs. 19&20), while from the computational cost point of view, TSF-FP2 is more expensive than the UOA methods due to the extra dimension (cf. Tabs. 18&19), especially the memory cost in high dimensions (cf. Tab. 18), and IEI-FP is found to be most efficient (cf. Tabs. 19&20&21).

(v) Among all the methods, in the $\varepsilon=O(1)$ regime, the EWI-FP method and TS-FP are the most accurate and efficient methods (cf. Tab. 19). While in the intermediate regime and the limit regime, the two UOA methods, i.e. IEI-FP and MTI-FP2 are significantly more powerful than the others (cf. Tabs. 20&21).

Remark 5.1**.**

The second order uniform accuracy of TSF-FP2 has been shown in [26] under condition $\tau\leq C$ with $C>0$ independent of $\varepsilon$. Here in our test (cf. Tab. 15), we are interested in performance of the scheme with a wide range of time step. Hence the order reduction here does not conflict with the theoretical results.

6 Applications

In this section, we apply the UOA MTI-FP2 method to study numerically the convergence rates from the NKGE (1.1) to its limiting models (1.5) and (1.6), and to simulate wave interaction in two dimensions (2D).

6.1 Convergence rates of NKGE to its limiting models

We take $d=1$ and $\lambda=1$ in the NKGE (1.1). Let $u$ be the solution of NKGE (1.1), $z_{\rm sw}$ be the solution of the NLSW (1.5) and $z_{\rm s}$ be the solution of the NLSE (1.6). Take the initial data as

[TABLE]

or

[TABLE]

where $m=1,2$. The solutions are obtained numerically with very fine mesh on a bounded interval $\Omega=(-128,128)$ with periodic boundary conditions. Define

[TABLE]

and define the error functions as

[TABLE]

Fig. 3 shows the errors defined in (6.3) as functions of time with the smooth initial data (6.1). Fig. 4 and Fig. 5 show the results from the nonsmooth initial data (6.2) with $m=2$ and $m=1$, respectively.

From Figs. 3-5, we can draw the following conclusions:

(i) The solution of the NKGE (1.1) converges to that of the NLSW (1.5) quadratically in $\varepsilon$ (and uniformly in time) provided that the initial data in (1.1) is smooth or at least satisfies $\phi_{1}$ and $\phi_{2}\in H^{2}(\Omega)$, i.e.

[TABLE]

where the constant $C_{0}>0$ is independent of $\varepsilon$ and time $t\geq 0$.

(ii) The solution of the NKGE (1.1) converges to that of the NLSE (1.6) quadratically in $\varepsilon$ (in general, not uniformly in time) provided that the initial data in (1.1) is smooth or at least satisfies $\phi_{1}$ and $\phi_{2}\in H^{3}(\Omega)$, i.e.

[TABLE]

where $C_{1}$ and $C_{2}$ are two positive constants which are independent of $\varepsilon$ and $T$. On the contrary, if the regularity of the initial data is weaker, e.g. $\phi_{1}$ and/or $\phi_{2}\in H^{2}(\Omega)$, then the convergence rate collapses to linear rate, i.e.

[TABLE]

where $C_{3}$ and $C_{4}$ are two positive constants which are independent of $\varepsilon$ and $T$. Rigorous mathematical justification for these numerical observations is on-going.

(iii) Under the same $\varepsilon$ and at the same time $t$, the error $\eta_{\rm sw}(t)$ is much small than $\eta_{\rm w}(t)$. It indicates that the NLSW (1.5) would be a better choice to design the LI scheme than the limit model (1.6), especially considering the long time behavior of the approximation.

6.2 Wave interactions in 2D

We take $d=2$ and $\lambda=1$ in the NKGE (1.1) and choose the initial data as

[TABLE]

The problem is solved numerically on a bounded computational domain $\Omega=(-16,16)\times(-16,16)$ with the periodic boundary condition. Fig. 6 shows contour plots of the solutions of the NKGE (1.1) in 2D under different $\varepsilon$.

7 Conclusions

We systematically studied and compared different numerical methods to solve the nonlinear Klein-Gordon equation (NKGE) in the nonrelativistic limit regime, while the solution is highly oscillatory in time in the limit regime. The numerical methods considered here include the classical finite difference time domain methods, the exponential wave integrator (EWI) spectral method, the time-splitting (TS) spectral method, the limit integrators, and the recently proposed uniformly accurate (UA) methods namely the multiscale time integrator (MTI) spectral method, the two-scale formulation (TSF) method and the iterative exponential integrator (IEI). We emphasized the finite time error bound of each method and the resolution capacity in terms of the oscillation wavelength in the limit regime. Systematical comparisons between the methods in the accuracy, computational complexity and other mathematical properties were carried out. Numerical experiments were done to show and compare the performance of each method from the classical regime to the limit regime. Our results show the EWI and TS methods are most efficient in the classical regime, while the UA methods are more powerful in the intermediate and limit regimes. Among the UA methods, the uniformly and optimally accurate methods are the most efficient and accurate for $\varepsilon\in(0,1]$ . Finally, the UA numerical methods were applied to study numerically the convergence rates of the NKGE (1.1) to its limiting models and to simulate wave interaction in two dimensions.

Acknowledgements

This work was supported by the Ministry of Education of Singapore grant R-146-000-223-112 (W. Bao) and the French ANR project MOONRISE ANR-14-CE23-0007-01 (X. Zhao).