Efficient Energy-Preserving Exponential Integrators for Multi-component Hamiltonian Systems

Abstract

In this paper, we develop a framework to construct energy-preserving methods for multi-component Hamiltonian systems, combining the exponential integrator and the partitioned averaged vector field method. This leads to numerical schemes with advantages of original energy conservation, long-time stability and excellent behavior for highly oscillatory or stiff problems. Compared to the existing energy-preserving exponential integrators (EP-EI) in practical implementation, our proposed methods are much efficient which can at least be computed by subsystem instead of handling a nonlinear coupling system at a time. Moreover, for most cases, such as the Klein-Gordon-Schrödinger equations and the Klein-Gordon-Zakharov equations considered in this paper, the computational cost can be further reduced. Specifically, one part of the derived schemes is totally explicit and the other is linearly implicit. In addition, we present rigorous proof of conserving the original energy of Hamiltonian systems, in which an alternative technique is utilized so that no additional assumptions are required, in contrast to the proof strategies used for the existing EP-EI. Numerical experiments are provided to demonstrate the significant advantages in accuracy, computational efficiency and the ability to capture highly oscillatory solutions.

Appendices

Appendix A Derivation of the EPAVF Scheme for 2D KGS Equations

Let the computation domain $\Omega =[a,b]\times [c,d]$ and $N_x$, $N_y$ be given even integers. The spatial steps are then defined as $h_x=(b-a)/N_x$, $h_y=(d-c)/N_y$ and the mesh grid is denoted by $\Omega _h=\{(x_i,y_j)\vert x_i=a+ih_x, y_j=c+jh_y,i=0,\cdots ,N_x, j=0,\cdots ,N_y\}$. Let $V_h=\{V\vert V=(V_{ij}),(x_i,y_j)\in \Omega _h\}$ be the space of grid function on $\Omega _h$. Applying the Fourier pseudospectral method to the 2D KGS equations and after some arrangements, we have

$$\begin{aligned} {\left\{ \begin{array}{ll} {Q}_t + \beta (\mathbb {D}_2^x {P}+{P}\mathbb {D}_2^y) + {P} \odot {U} = 0, \\ {P}_t - \beta (\mathbb {D}_2^x {Q}+{Q}\mathbb {D}_2^y)- {Q} \odot {U} = 0, \\ {U}_t = {V},\\ \varepsilon ^{2}{V}_t - (\mathbb {D}_2^x {U} +{U}\mathbb {D}_2^y ) + \dfrac{1}{\varepsilon ^{2}} {U} - {Q}^{2} - {P}^{2} = 0, \end{array}\right. } \end{aligned}$$

(A1)

where $\mathbb {D}_2^\alpha , \alpha =x,y$ are the second-order spectral differential matrices in x, y directions, respectively. For convenience of the derivation, we reshape U by columns into a vector u of dimension $N_x\times N_y$, etc. Then the above system (A1) can be reformed as

$$\begin{aligned} {\left\{ \begin{array}{ll} \mathbf{q }_t + \beta \mathbb {D} \mathbf{p } + \mathbf{p } \odot \mathbf{u } = 0, \\ \mathbf{p }_t - \beta \mathbb {D} \mathbf{q }- \mathbf{q } \odot \mathbf{u } = 0, \\ \mathbf{u }_t = \mathbf{V },\\ \varepsilon ^{2}{v}_t - \mathbb {D} \mathbf{U } + \dfrac{1}{\varepsilon ^{2}} \mathbf{u } - \mathbf{q }^{2} - \mathbf{p }^{2} = 0, \end{array}\right. } \end{aligned}$$

(A2)

where $\mathbb {D} = I_{N_y}\otimes \mathbb {D}_2^x + \mathbb {D}_2^y \otimes I_{N_x}$. The only differences of the semi-discretization (21) and (A2) appear in the dimension of variables and the differential matrix $\mathbb {D}$. Therefore, the corresponding EPAVF scheme can be similarly obtained as (28), where the matrix exponentials are calculated by

$$\begin{aligned} \begin{aligned} \exp (V_1) = \left( \begin{array}{cc} \cos (\tau \beta \mathbb {D}) &{} -\sin (\tau \beta \mathbb {D})\\ \sin (\tau \beta \mathbb {D}) &{} \cos (\tau \beta \mathbb {D}) \end{array}\right) , \quad \exp (V_2)= \left( \begin{array}{cc} \cos (\tau \widetilde{\mathbb {D}})&{} \frac{\sin ( \tau \widetilde{\mathbb {D}}) }{ \widetilde{\mathbb {D}}} \\ -\widetilde{\mathbb {D}}\sin ( \tau \widetilde{\mathbb {D}}) &{} \cos (\tau \widetilde{\mathbb {D}}) \end{array} \right) , \\ \end{aligned} \end{aligned}$$

(A3)

and

$$\begin{aligned} \begin{aligned}&\varphi (V_1) = \left( \begin{array}{cc} \frac{\sin (\tau \beta \mathbb {D})}{\tau \beta \mathbb {D}} &{} \frac{I_{N_x\times N_y} - \cos (\tau \beta \mathbb {D})}{\tau \beta \mathbb {D}} \\ -\frac{I_{N_x\times N_y} - \cos (\tau \beta \mathbb {D})}{\tau \beta \mathbb {D}} &{} \frac{\sin (\tau \beta \mathbb {D})}{\tau \beta \mathbb {D}} \end{array}\right) ,\\&\varphi (V_2) = \left( \begin{array}{cc} \frac{\sin ( \tau \widetilde{\mathbb {D}} ) }{\tau \widetilde{\mathbb {D}}} &{} \frac{ I_{N_x\times N_y} - \cos (\tau \widetilde{\mathbb {D}}) }{\tau \widetilde{\mathbb {D}}^{2}} \\ \frac{\cos ( \tau \widetilde{\mathbb {D}} ) - I_{N_x\times N_y} }{\tau } &{} \frac{\sin ( \tau \widetilde{\mathbb {D}} ) }{\tau \widetilde{\mathbb {D}}} \end{array}\right) . \end{aligned} \end{aligned}$$

(A4)

Notice that the differential matrix $\mathbb {D}$ can be decomposed as

$$\begin{aligned} \mathbb {D}= (F_{N_y}^{-1} \otimes F_{N_x}^{-1})(I_{N_y}\otimes \Lambda _2^x+\Lambda _2^y \otimes I_{N_x} )(F_{N_y} \otimes F_{N_x}), \end{aligned}$$

so that the calculation of trigonometric functions of matrix $\mathbb {D}$ can also be accelerated by FFT. For example, we have

$$\begin{aligned} \cos {(\tau \beta \mathbb {D})} = (F_{N_y}^{-1} \otimes F_{N_x}^{-1})\cos ( \tau \beta ( I_{N_y}\otimes \Lambda _2^x + \Lambda _2^y \otimes I_{N_x}) ) (F_{N_y} \otimes F_{N_x}). \end{aligned}$$

Moreover, due to the fact that $(B^\top \otimes A)\mathrm{vec}(X)=\mathrm{vec}(AXB)$ for any suitable matrices A, B, X, where the function ‘$\mathrm{vec}(\cdot )$’ represents the vectorization of a matrix, we have

$$\begin{aligned} \cos {(\tau \beta \mathbb {D})} {\mathbf {u}} =\mathrm{vec}\Big ( F_{N_x}^{-1}( \cos {(\tau \beta \Lambda )}\odot \widetilde{{U}} ) F_{N_y}^{-\top } \Big ), \end{aligned}$$

where $\Lambda _{ij} = (\Lambda _2^x)_{ii} + (\Lambda _2^y)_{jj}$, $\widetilde{{U}} = F_{N_x} {U} F_{N_y}^\top$ and $(\cos {(\tau \beta \Lambda )})_{ij} = \cos {(\tau \beta \Lambda _{ij})}$. As a consequence, in practical computation we implement the EPAVF scheme (28) for 2D KGS equations just in the following matrix form

$$\begin{aligned} {\left\{ \begin{array}{ll} {\widetilde{Q}}^{n+1} = \exp _{11}^1\odot {\widetilde{Q}}^{n} + \exp _{12}^1\odot \widetilde{{P}}^{n} + \frac{\tau }{2} \varphi _{12}^1 \odot {\widetilde{G}}_1 - \frac{\tau }{2} \varphi _{11}^1 \odot {\widetilde{G}}_2, \\ {\widetilde{P}}^{n+1} = \exp _{21}^1\odot {\widetilde{Q}}^{n} + \exp _{22}^1\odot {\widetilde{P}}^{n} + \frac{\tau }{2} \varphi ^1_{22} \odot {\widetilde{G}}_1 - \frac{\tau }{2} \varphi ^1_{21} \odot {\widetilde{G}}_2, \\ {\widetilde{U}}^{n+1} =\exp _{11}^2\odot {\widetilde{U}}^{n} + \exp _{12}^2\odot {\widetilde{V}}^{n} + \frac{\tau }{\varepsilon ^{2}}\varphi _{12}^2 \odot {\widetilde{G}}_3, \\ {\widetilde{V}}^{n+1} = \exp _{21}^2\odot {\widetilde{U}}^{n} + \exp _{22}^2\odot {\widetilde{V}}^{n} + \frac{\tau }{\varepsilon ^{2}}\varphi _{22}^2 \odot {\widetilde{G}}_3, \end{array}\right. } \end{aligned}$$

(A5)

where ${G}_1=U^{n} \odot Q^{n+1/2}, {G}_2=U^{n} \odot P^{n+1/2}, {G}_3=\left( Q^{n+1}\right) ^2 + \left( P^{n+1}\right) ^2$. The components of matrix exponentials, i.e., $\exp _{ij}^k$, $\varphi _{ij}^k$ are obtained by replacing $\mathbb {D}$ with $\Lambda$ and the corresponding trigonometric functions are computed element-by-element.

Although we have only derived the practical EPAVF scheme (A5) for the 2D KGS equations here, its adjoint scheme as well the EPAVF schemes for 2D KGZ equations can be obtained similarly. Furthermore, the procedures of derivation can be directly generalized to 3D cases with few changes.

Appendix B Convergence Analysis of EPAVF for 1D KGS Equations

Due to the limitation of space, we only consider the convergence analysis of EPAVF for 1D KGS equations. In addition, since we aim to investigate the $\varepsilon$-scalability from the convergence analysis, the spatial discretizations are conducted by Fourier series such that the corresponding errors can be negligible.

We first express the unknowns into the Fourier series, such as

$$\begin{aligned} q(x, t) = \sum \limits _{k \in \mathbb {Z}} {\widehat{q}}_k(t) e^{i \mu _k (x-a)}, \quad p(x, t)u(x, t) = \sum \limits _{k \in \mathbb {Z}} \widehat{(pu)}_k(t) e^{i \mu _k (x-a)}, \end{aligned}$$

where $\mu _k = \frac{2\pi k}{b-a}$ and ${\widehat{q}}_k (t)$, $\widehat{(pu)}_k(t)$ are corresponing Fourier coefficients, respectively. Further noting that ${\widehat{\psi }}={\widehat{q}}+i{\widehat{p}}$, we then derive the semi-discretization of the KGS equations (20) in a compact form as follows:

$$\begin{aligned} \left\{ \begin{aligned}&i{\widehat{\psi }}_k^\prime (t) = -\beta \mu _k^2 {\widehat{\psi }}_k(t) + \widehat{(p\psi )}_k(t), \\&{\widehat{u}}_k^\prime (t) = {\widehat{v}}_k(t), \\&\varepsilon ^2 {\widehat{v}}_k^\prime (t) = -\mu _k^2 {\widehat{u}}_k(t) - \dfrac{1}{\varepsilon ^2}{\widehat{u}}_k(t) + \widehat{(q^2)}_k(t) + \widehat{(p^2)}_k(t). \end{aligned}\right. \end{aligned}$$

(B1)

Applying the EPAVF method to (B1) yields

$$\begin{aligned} \begin{aligned}&{\widehat{\psi }}^{n+1}_k = e^{-i\tau \beta \mu _k^2} {\widehat{\psi }}^n_k + \alpha _k(\tau ) \widehat{(u^n \psi ^{n+\frac{1}{2}})}_k, \\&{\widehat{u}}^{n+1}_k = \cos {(\tau \omega _k)} {\widehat{u}}^n_k + \frac{\sin {(\tau \omega _k)}}{\omega _k} {\widehat{v}}^n_k + \beta _k(\tau ) \widehat{(\vert \psi ^{n+1}\vert ^2)}_k, \\&{\widehat{v}}^{n+1}_k = -\omega _k \sin {(\tau \omega _k)} {\widehat{u}}^n_k + \cos {(\tau \omega _k)}{\widehat{v}}^n_k + \gamma _k(\tau ) \widehat{(\vert \psi ^{n+1}\vert ^2)}_k, \\ \end{aligned} \end{aligned}$$

(B2)

where $\omega _k = \frac{1}{\varepsilon ^2} \sqrt{1 + \varepsilon ^2\mu _l^2}$ and the notations $\alpha _k(\tau )$, $\beta _k(\tau )$, $\gamma _k(\tau )$ are analogous to the elements in (4.6) and (4.7). In order to obtain their bound, we adopt the integral forms as

$$\begin{aligned} \begin{aligned}&\alpha _k(\tau ) = \int _0^\tau e^{-i(s- \tau )\beta \mu _k^2} ds, \quad \beta _k(\tau ) = \int _0^\tau \dfrac{\sin {(\omega _k (\tau - s))}}{\varepsilon ^2 \omega _k}ds, \\&\gamma _k(\tau ) = \int _0^\tau \dfrac{\cos {(\omega _k (\tau - s))}}{\varepsilon ^2} ds, \end{aligned} \end{aligned}$$

such that one can quickly obtain

$$\begin{aligned} \vert \alpha _k(\tau )\vert \le \tau , \ \vert \beta _k(\tau )\vert \le \dfrac{\tau }{\varepsilon ^2 \omega _k}, \ \vert \gamma _k(\tau )\vert \le \dfrac{\tau }{\varepsilon ^2}. \end{aligned}$$

(B3)

We denote $H_p^m(\Omega ) := \big \{ u \in H^m(\Omega ) : u^{(k)}(a) = u^{(k)}(b), \ 0 \le k \le m-1 \big \}.$ The semi-norm and the norm of $H_p^m(\Omega )$ are characterized in the frequency space [40] by

$$\begin{aligned} \vert u\vert _m = \big (\sum \limits _{k \in \mathbb {Z}} \mu _k^{2m} \vert {\widehat{u}}_k \vert ^2\big )^\frac{1}{2}, \ \Vert u\Vert _m = \big (\sum \limits _{k \in \mathbb {Z}} (1 + \mu _k^2)^m \vert {\widehat{u}}_k \vert ^2 \big )^\frac{1}{2}. \end{aligned}$$

We will use C to represent a generic positive constant independent of $\varepsilon$ and the discretization parameters. Motivated by the asymptotic results in [3] as well as the numerical results, we make the following assumptions on the analytic solutions of (19), that is, $u, u_t \in L^\infty ([0, T]; H_p^1(\Omega ))$, $\psi \in L^\infty ([0, T]; H_p^3(\Omega ))$, $\psi _t \in L^\infty ([0, T]; H_p^1(\Omega ))$ and they are bounded in the sense

$$\begin{aligned} \begin{aligned}&\Vert \psi \Vert _{L^\infty ([0, T]; H^3)} + \Vert \partial _t \psi \Vert _{L^\infty ([0, T]; H^1)} + \Vert u\Vert _{L^\infty ([0, T]; H^1)} + \varepsilon ^2 \Vert \partial _t u\Vert _{L^\infty ([0, T]; H^1)} \le M, \\ \end{aligned}\nonumber \\ \end{aligned}$$

(B4)

where M is a positive number independent of $\varepsilon$ and mesh size. Notice that from the Schrödinger part in (19), we can derive

$$\begin{aligned} i \partial _t \rho + \partial _x [ {\overline{\psi }}\partial _x \psi - \psi \partial _x {\overline{\psi }} ] = 0, \end{aligned}$$

where $\rho := \vert \psi \vert ^2$. Consequently, from the assumption (B4), we have $\partial _t \rho \in L^\infty ([0, T]; H^1_p(\Omega ))$ and

$$\begin{aligned} \Vert \rho \Vert _{L^\infty ([0, T]; H^1)} + \Vert \partial _t\rho \Vert _{L^\infty ([0, T]; H^1)} \le M. \end{aligned}$$

(B5)

We first give an estimate of the local error. According to the variation-of-constant formula, the exact solution of (B1), from $t_n$ to $t_{n+1}$, can be expressed as

$$\begin{aligned} \begin{aligned}&{\widehat{\Psi }}^{n+1}_k = e^{-i\tau \beta \mu _k^2}{\widehat{\Psi }}^n_k + \int _0^\tau e^{-i(s-\tau )\beta \mu _k^2} \widehat{(U \Psi )}_k(t_n + s) ds, \\&{\widehat{U}}^{n+1}_k = \cos {(\tau \omega _k)} {\widehat{U}}^n_k + \dfrac{\sin {(\tau \omega _k)}}{\omega _k}{\widehat{V}}^n_k + \int _0^\tau \dfrac{\sin {(\omega _k(\tau - s))}}{\varepsilon ^2 \omega _k} \widehat{(\vert \Psi \vert ^2)}_k(t_n + s) ds, \\&{\widehat{V}}^{n+1}_k = -\omega _k \sin {(\tau \omega _k)}{\widehat{U}}^n_k + \cos {(\tau \omega _k)} {\widehat{V}}^n_k + \int _0^\tau \dfrac{\cos {(\omega _k (\tau - s))}}{\varepsilon ^2} \widehat{(\vert \Psi \vert ^2)}_k (t_n + s) ds, \\ \end{aligned} \end{aligned}$$

(B6)

where the capital letters represent the analytic solutions. We then define the local error functions as

$$\begin{aligned} \xi _\psi ^n = \sum \limits _{k \in \mathbb {Z}} \widehat{(\xi _\psi ^n)}_k e^{i\mu _kx}, \ \xi _u^n = \sum \limits _{k \in \mathbb {Z}} \widehat{(\xi _u^n)}_k e^{i\mu _kx}, \ \xi _v^n= \sum \limits _{k \in \mathbb {Z}} \widehat{(\xi _v^n)}_k e^{i\mu _kx}, \end{aligned}$$

where

$$\begin{aligned} \begin{aligned}&\widehat{(\xi _\psi ^n)}_k := {\widehat{\Psi }}^{n+1}_k - e^{-i\tau \beta \mu _k^2} {\widehat{\Psi }}^n_k - \alpha _k(\tau ) \widehat{(U^n\Psi ^{n+\frac{1}{2}})}_k, \\&\widehat{(\xi _u^n)}_k := {\widehat{U}}^{n+1}_k - \cos {(\tau \omega _k)} {\widehat{U}}^n_k - \dfrac{\sin {(\tau \omega _k)}}{\tau \omega _k} {\widehat{V}}^n_k - \beta _k(\tau ) \widehat{(\vert \Psi ^{n+1} \vert ^2)}_k, \\&\widehat{(\xi _v^n)}_k := {\widehat{V}}^{n+1}_k + \omega _k \sin {(\tau \omega _k)} {\widehat{U}}^n_k - \cos {(\tau \omega _k)} {\widehat{V}}^n_k - \gamma _k(\tau ) \widehat{(\vert \Psi ^{n+1} \vert ^2)}_k. \end{aligned} \end{aligned}$$

(B7)

Lemma 2

(Estimates for local errors). Under the assumption (B4) and (B5), we have the following estimates for the local errors:

$$\begin{aligned} \varepsilon ^2\Vert \xi _\psi ^{n}\Vert _1 + \Vert \xi _u^{n} \Vert _{1} + \varepsilon \Vert \partial _x \xi _u^{n} \Vert _1 + \varepsilon ^2 \Vert \xi _v^{n}\Vert _1 \le C \tau ^2. \end{aligned}$$

(B8)

Proof

It follows from (B6) and (B7) that

$$\begin{aligned} \begin{aligned}&\widehat{(\xi _\psi ^{n})}_k = \int _0^\tau e^{-i(s - \tau )\beta \mu _k} \widehat{(\eta _1)}_k ds, \ \widehat{(\xi _u^{n})}_k = \int _0^\tau \dfrac{\sin {(\omega _k (\tau - s))}}{\varepsilon ^2 \omega _k} \widehat{(\eta _2)}_k ds, \\&\widehat{(\xi _v^{n})}_k = \int _0^\tau \dfrac{\cos {(\omega _k(\tau - s))}}{\varepsilon ^2} \widehat{(\eta _2)}_k ds, \end{aligned} \end{aligned}$$

(B9)

where

$$\begin{aligned} \eta _1 = (U\Psi )(t_n+s) - U^n \Psi ^{n+\frac{1}{2}}, \ \eta _2 = \vert \Psi (t_n+s) \vert ^2 - \vert \Psi ^{n+1} \vert ^2. \end{aligned}$$

By the Taylor’s expansion with integral remainder, we obtain

$$\begin{aligned} \begin{aligned} {\eta }_1&= s \int _0^1 \partial _t(U\Psi )(t_n + \theta s) d\theta - \dfrac{\tau }{2} U^n \int _0^1 \partial _t \Psi (t_n + \theta \tau ) d\theta , \\ {\eta }_2&= (s - \tau ) \int _0^1 \partial _t (\vert \Psi \vert ^2) (t_n+\tau + (s-\tau )\theta ) d\theta . \end{aligned} \end{aligned}$$

Squaring on both sides of the first equation in (B9) and multiplying the resulting equation by $1 + \mu _k^2$ and then summing up with respect to k, we get

$$\begin{aligned} \begin{aligned} \Vert \xi _{\psi }^{n}\Vert _{1}^2&\le 2\tau ^4 \Vert \partial _t(u \psi )\Vert ^2_{L^\infty ([0, T]; H^1)} + \dfrac{\tau ^4}{4} \Vert u\Vert ^2_{L^\infty ([0, T]; H^1)}\\&\quad + \dfrac{\tau ^4}{4} \Vert \partial _t\psi \Vert ^2_{L^\infty ([0, T]; H^1)} \le C \dfrac{\tau ^4}{\varepsilon ^4}, \end{aligned} \end{aligned}$$

(B10)

where the assumption (B4) is used. Similarly, from the rest two equations in (B9) and the assumption (B5), we have

$$\begin{aligned} \Vert \xi _u^{n}\Vert _{1}^2 \le C \tau ^4, \ \Vert \partial _x \xi _u^{n} \Vert _{1}^2 \le C \dfrac{\tau ^4}{\varepsilon ^2}, \ \Vert \xi _v^{n}\Vert _{1}^2 \le C \dfrac{\tau ^4}{\varepsilon ^4}. \end{aligned}$$

(B11)

Combining (B10) and (B11), we get the estimate (B8). $\square$

Let $\psi ^n, u^n, v^n$ be the numerical solutions of the EPAVF method (B2). Since the spatial discretization is carried out by Fourier serials without any truncation, the spatial errors can be negligible and thus we can only consider the error estimate induced by the time discretization.

Theorem 4

(Convergence) Under the assumptions (B4)-(B5) and for sufficiently small time step $\tau$ such that $C \dfrac{\tau }{\varepsilon ^2} \le 1$, we have

$$\begin{aligned} \Vert \psi ^n - \Psi ^n \Vert _1 \le C \dfrac{\tau }{\varepsilon ^2}, \ \Vert u^n - U^n \Vert _1 \le C \dfrac{\tau }{\varepsilon ^2}, \ \Vert v^n - V^n \Vert _1 \le C \dfrac{\tau }{\varepsilon ^4}. \end{aligned}$$

(B12)

Consequently, the numerical solution $\psi ^n$, $u^n$ and $v^n$ are uniformly bounded in the sense

$$\begin{aligned} \begin{aligned} \Vert \psi ^n\Vert _1 \le M + 1, \ \Vert u^n\Vert _1 \le M + 1, \ \Vert v^n\Vert _1 \le \dfrac{M + 1}{\varepsilon ^2}. \end{aligned} \end{aligned}$$

(B13)

Proof

The following proof will be proceeded by the mathematical induction. Define the errors and notations as

$$\begin{aligned} \begin{aligned}&\varepsilon _\psi ^n := \psi ^n - \Psi ^n, \ \varepsilon _u^n := u^n - U^n, \ \varepsilon _v^n := v^n - V^n, \\&\delta _1^{n+1} := u^n \psi ^{n+\frac{1}{2}} - U^n \Psi ^{n+\frac{1}{2}}, \ \delta _2^{n+1} = \vert \psi ^{n+1}\vert ^2 - \vert \Psi ^{n+1} \vert ^2, \end{aligned} \end{aligned}$$

which satisfies the error equations

$$\begin{aligned} \begin{aligned} \widehat{(\varepsilon _\psi ^{n+1})}_k&= e^{-i\tau \beta \mu _k^2} \widehat{(\varepsilon _\psi ^n)}_k + \alpha _k(\tau ) \widehat{(\delta _1^n)}_k + \widehat{(\xi _\psi ^{n})}_k, \\ \widehat{(\varepsilon _u^{n+1})}_k&= \cos {(\tau \omega _k)} \widehat{(\varepsilon _u^n)}_k + \dfrac{\sin {(\tau \omega _k)}}{\omega _k} \widehat{(\varepsilon _v^n)}_k + \beta _k (\tau ) \widehat{(\delta _2^n)}_k + \widehat{(\xi _u^{n})}_k, \\ \widehat{(\varepsilon _v^{n+1})}_k&= -\omega _k \sin {(\tau \omega _k)} \widehat{(\varepsilon _u^n)}_k + \cos {(\tau \omega _k)} \widehat{(\varepsilon _v^n)}_k + \gamma _k (\tau )\widehat{(\delta _2^n)}_k + \widehat{(\xi _v^{n})}_k. \end{aligned} \end{aligned}$$

(B14)

It is clear that the results (B12) and (B13) are valid for $n = 0$. Suppose these results hold for $n = 1$, $\cdots$, l. Then we will prove for $n = l+1$, we also have the estimates for the errors and the corresponding boundedness.

We first prove $\psi ^{l+1}$ is bounded under the $H^1$-norm. Squaring on both sides of the first equation of (B2) and by (B13), we have

$$\begin{aligned} \begin{aligned} \vert {\widehat{\psi }}^{l+1}_k\vert ^2&\le \vert {\widehat{\psi }}^l_k\vert ^2 + 2\tau \vert {\widehat{\psi }}^l_k \vert \vert \widehat{(u^l \psi ^{l+\frac{1}{2}})}_k \vert + \tau ^2 \vert \widehat{(u^l \psi ^{l+\frac{1}{2}})}_k \vert ^2 \\&\le (1 + \tau ) \vert {\widehat{\psi }}^l_k\vert ^2 + \tau (1 + \tau ) \vert \widehat{(u^l \psi ^{l+\frac{1}{2}})}_k \vert ^2. \end{aligned} \end{aligned}$$

(B15)

Let $\tau \le 1$ and sum up the above equation for $k \in \mathbb {Z}$. We obtain

$$\begin{aligned} \begin{aligned} \Vert \psi ^{l+1}\Vert ^2&\le (1 + \tau ) \Vert \psi ^l\Vert ^2 + 2\tau \Vert u^l \psi ^{l+\frac{1}{2}}\Vert ^2 \\&\le (1 + \tau ) \Vert \psi ^l\Vert ^2 + 2\tau \Vert u^l\Vert _\infty ^2 \Vert \psi ^{l+\frac{1}{2}} \Vert ^2 \\&\le \Vert \psi ^l\Vert ^2 + C \tau (\Vert \psi ^{l+1}\Vert _1^2 + \Vert \psi ^l\Vert _1^2), \end{aligned} \end{aligned}$$

(B16)

where the one-dimensional Sobolev inequality $\Vert u\Vert _\infty ^2 \le C \Vert u\Vert _1^2$ is utilized. Similarly, according to the definition of semi-norm and from (B15), we can derive

$$\begin{aligned} \begin{aligned} \vert \psi ^{l+1}\vert _1^2&\le (1 + \tau )\vert \psi ^l\vert _1^2 + 4 \tau \left( \Vert \partial _x u^l \psi ^{l+\frac{1}{2}} \Vert ^2 + \Vert u^l \partial _x \psi ^{l+\frac{1}{2}} \Vert ^2 \right) \\&\le (1 + \tau )\vert \psi ^l\vert _1^2 + 4 \tau \left( \vert u^l\vert _1^2 \Vert \psi ^{l+\frac{1}{2}} \Vert _\infty ^2 + \Vert u^l\Vert _\infty ^2 \vert \psi ^{l+\frac{1}{2}}\vert _1^2 \right) \\&\le \vert \psi ^l\vert _1^2 + C\tau (\Vert \psi ^{l+1}\Vert _1^2 + \Vert \psi ^l\Vert _1^2). \end{aligned} \end{aligned}$$

(B17)

Adding (B16) and (B17) together, we have

$$\begin{aligned} \Vert \psi ^{l+1}\Vert _1^2 \le \Vert \psi ^l\Vert _1^2 + C \tau (\Vert \psi ^{l+1}\Vert _1^2 + \Vert \psi ^l\Vert _1^2), \end{aligned}$$

which implies

$$\begin{aligned} \Vert \psi ^{l+1}\Vert _1^2 \le C \Vert \psi ^l\Vert _1^2, \end{aligned}$$

providing $\tau$ is sufficiently small. Then based on the induction, we obtain

$$\begin{aligned} \Vert \psi ^{l+1}\Vert _1^2 \le C, \quad \Vert \psi ^{l+1}\Vert _\infty ^2 \le C. \end{aligned}$$

Next, we will give the $H^1$-estimates for the nonlinear terms $\delta _1^l$ and $\delta _2^l$. In view of their definitions, we can easily get

$$\begin{aligned} \begin{aligned} \Vert \delta _1^l \Vert ^2&\le 2 (\Vert u^l\Vert _\infty ^2\Vert \varepsilon _\psi ^{l+\frac{1}{2}}\Vert ^2 + \Vert \Psi ^{l+\frac{1}{2}}\Vert _\infty ^2 \Vert \varepsilon _u^l\Vert ^2 ) \le C( \Vert \varepsilon _\psi ^l\Vert ^2 + \Vert \varepsilon _\psi ^{l+1}\Vert ^2 + \Vert \varepsilon _u^l\Vert ^2), \end{aligned} \end{aligned}$$

and

$$\begin{aligned} \begin{aligned} \Vert \partial _x \delta _1^l\Vert ^2&\le 4(\vert u^l\vert _1^2 \Vert \varepsilon _\psi ^{l+\frac{1}{2}}\Vert _\infty ^2 + \vert \varepsilon _u^l\vert _1^2 \Vert \Psi ^{l+\frac{1}{2}} \Vert _\infty ^2 + \Vert u^l\Vert _{\infty }^2 \vert \varepsilon _\psi ^{l+\frac{1}{2}}\vert _1^2 + \Vert \varepsilon _u^l\Vert ^2 \Vert \Psi ^{l+\frac{1}{2}}\Vert _\infty ^2) \\&\le C ( \Vert \varepsilon _\psi ^l\Vert _1^2 + \Vert \varepsilon _\psi ^{l+1}\Vert _1^2 + \Vert \varepsilon _u^l\Vert _1^2 ). \end{aligned} \end{aligned}$$

A combination of the above two equations yields

$$\begin{aligned} \Vert \delta _1^l\Vert _1^2 \le C ( \Vert \varepsilon _\psi ^l\Vert _1^2 + \Vert \varepsilon _\psi ^{l+1}\Vert _1^2 + \Vert \varepsilon _u^l\Vert _1^2 ). \end{aligned}$$

(B18)

Similarly, we can derive

$$\begin{aligned} \begin{aligned} \Vert \delta _2^l \Vert ^2 \le C \Vert \varepsilon _\psi ^{l+1}\Vert ^2, \Vert \partial _x \delta _2^l \Vert ^2 \le C \Vert \varepsilon _\psi ^{l+1}\Vert _1^2, \end{aligned} \end{aligned}$$

and subsequently,

$$\begin{aligned} \Vert \delta _2^l\Vert _1^2 \le C \Vert \varepsilon _\psi ^{l+1}\Vert _1^2. \end{aligned}$$

(B19)

Finally, we shall prove the convergence of the EPAVF method (B2). Squaring on both sides of the first error equation in (B14) and by the boundedness results (B13), we have

$$\begin{aligned} \begin{aligned} \vert \widehat{(\varepsilon _\psi ^{l+1})_k}\vert ^2&\le \vert \widehat{(\varepsilon _\psi ^l)}_k\vert ^2 + 2 \vert \widehat{(\varepsilon _\psi ^l)}_k\vert \vert \alpha _k(\tau ) \widehat{(\delta _1^l)}_k + \widehat{(\xi _\psi ^{l})}_k\vert + \vert \alpha _k(\tau ) \widehat{(\delta _1^l)}_k + \widehat{(\xi _\psi ^{l})}_k\vert ^2 \\&\le (1 + \tau )\vert \widehat{(\varepsilon _\psi ^l)}_k\vert ^2 + \dfrac{1 + \tau }{\tau } \vert \alpha _k(\tau ) \widehat{(\delta _1^l)}_k + \widehat{(\xi _\psi ^{l})}_k\vert ^2 \\&\le (1 + \tau )\vert \widehat{(\varepsilon _\psi ^l)}_k\vert ^2 + 2\tau (1 + \tau ) \vert \widehat{(\delta _1^l)}_k \vert ^2 + 2\dfrac{1 + \tau }{\tau } \vert \widehat{(\xi _\psi ^{l})}_k\vert ^2. \end{aligned} \end{aligned}$$

Multiplying the above inequality by $1 + \mu _k^2$ and summing up for $k\in \mathbb {Z}$, then utilizing the results in (B18) and local error estimates in Lemma 2, we get

$$\begin{aligned} \Vert \varepsilon _\psi ^{l+1}\Vert _1^2 \le \Vert \varepsilon _\psi ^l\Vert _1^2 + C \tau (\Vert \varepsilon _\psi ^{l+1}\Vert _1^2 + \Vert \varepsilon _\psi ^l\Vert _1^2 + \Vert \varepsilon _u^l\Vert _1^2) + C \dfrac{\tau ^3}{\varepsilon ^4}. \end{aligned}$$

(B20)

Similar processes can be employed for the rest error equations and give

$$\begin{aligned} \begin{aligned} \vert \widehat{(\varepsilon _u^{l+1})}_k\vert ^2&\le (1 + \tau ) \left| \cos {(\tau \omega _k)} \widehat{(\varepsilon _u^l)}_k + \dfrac{\sin {(\tau \omega _k)}}{\omega _k} \widehat{(\varepsilon _v^l)}_k \right| ^2 \\&\quad + \dfrac{2\tau (1 + \tau )}{\varepsilon ^4 \omega _k^2} \vert \widehat{(\delta _2^l)}_k\vert ^2 + 2\dfrac{1+\tau }{\tau } \vert \widehat{(\xi _u^{l})}_k\vert ^2, \end{aligned} \end{aligned}$$

(B21)

and

$$\begin{aligned} \begin{aligned} \vert \widehat{(\varepsilon _v^{l+1})}_k\vert ^2&\le (1 + \tau ) \left| -\omega _k\sin {(\tau \omega _k)} \widehat{(\varepsilon _u^l)}_k + \cos {(\tau \omega _k)}\widehat{(\varepsilon _v^l)}_k \right| ^2 \\&\quad + \dfrac{2\tau (1 + \tau )}{\varepsilon ^4} \vert \widehat{(\delta _2^l)}_k\vert ^2 + 2\dfrac{1+\tau }{\tau } \vert \widehat{(\xi _v^{l})}_k\vert ^2. \end{aligned} \end{aligned}$$

(B22)

Multiplying (B21) and (B22) by $(1 + \mu _k^2)(1 + \varepsilon ^2 \mu _k^2)$ and $(1 + \mu _k^2)\varepsilon ^4$, respectively and adding them together, then summing up for $k\in \mathbb {Z}$, we obtain

$$\begin{aligned} \begin{aligned}&\Vert \varepsilon _u^{l+1}\Vert _1^2 + \varepsilon ^2 \Vert \partial _x \varepsilon _u^{l+1}\Vert _1^2 + \varepsilon ^4 \Vert \varepsilon _v^{l+1}\Vert _1^2 \\&\quad \le (1 + \tau )(\Vert \varepsilon _u^l\Vert _1^2 + \varepsilon ^2 \Vert \partial _x \varepsilon _u^l\Vert _1^2 + \varepsilon ^4 \Vert \varepsilon _v^l\Vert _1^2) \\&\qquad +4\tau (1+\tau ) \Vert \delta _2^l\Vert _1^2 + \dfrac{2(1+\tau )}{\tau }( \Vert \xi _u^{l}\Vert _1^2 + \varepsilon ^2 \Vert \partial _x \xi _u^{l}\Vert _1^2 + \varepsilon ^4 \Vert \xi _v^{l}\Vert _1^2 ) \\&\quad \le (1+\tau ) ( \Vert \varepsilon _u^l\Vert _1^2 + \varepsilon ^2 \Vert \partial _x \varepsilon _u^l\Vert _1^2 + \varepsilon ^4 \Vert \varepsilon _v^l\Vert _1^2 ) + C \tau \Vert \varepsilon _\psi ^{l+1}\Vert _1^2 + C \tau ^3. \end{aligned} \end{aligned}$$

(B23)

Adding (B20) and (B23) together and denoting the error energy functional

$$\begin{aligned} {\mathcal {E}}(\varepsilon _\psi , \varepsilon _u, \varepsilon _v) := \Vert \varepsilon _\psi \Vert _1^2 + \Vert \varepsilon _u\Vert _1^2 + \varepsilon ^2 \Vert \partial _x \varepsilon _u\Vert _1^2 + \varepsilon ^4 \Vert \varepsilon _v\Vert _1, \end{aligned}$$

we derive

$$\begin{aligned} \begin{aligned}&{\mathcal {E}} (\varepsilon _\psi ^{l+1}, \varepsilon _u^{l+1}, \varepsilon _v^{l+1}) - {\mathcal {E}}(\varepsilon _\psi ^l, \varepsilon _u^l, \varepsilon _v^l) \le C \tau ( {\mathcal {E}} (\varepsilon _\psi ^{l+1}, \varepsilon _u^{l+1}, \varepsilon _v^{l+1}) + {\mathcal {E}} (\varepsilon _\psi ^{l}, \varepsilon _u^{l}, \varepsilon _v^{l}) ) + C \dfrac{\tau ^3}{\varepsilon ^4}. \end{aligned} \end{aligned}$$

Since the time step is assumed to be sufficient small, an application of the discrete Gronwall inequality yields

$$\begin{aligned} {\mathcal {E}}(\varepsilon _\psi ^l, \varepsilon _u^l, \varepsilon _v^l) \le C \dfrac{\tau ^2}{\varepsilon ^4}. \end{aligned}$$

Consequently, we have

$$\begin{aligned} \Vert \varepsilon _\psi ^{l+1}\Vert _1 \le C \dfrac{\tau }{\varepsilon ^2}, \ \Vert \varepsilon _u^{l+1}\Vert _1 \le C \dfrac{\tau }{\varepsilon ^2}, \ \Vert \varepsilon _v^{l+1}\Vert _1 \le C \dfrac{\tau }{\varepsilon ^4}, \end{aligned}$$

as well as the boundedness (B13). Therefore, we complete the proof for the case when $n = l+1$ and thus finish the whole convergence analysis. $\square$

Remark 9

We emphasize that the above error estimate for the EPAVF method can not only be generalized to the 2D case by replacing the error bound under the $H^2$-norm, but also provide the main procedures to analyze the convergence when the spatial derivatives are discretized by the Fourier pseudospectral method as we have employed in the paper. For the detailed process, one can refer to the related parts in [3].

Remark 10

The presented convergence analysis is only for the EPAVF method. We note that similar results can also be obtained for its adjoint with slight modifications of the proof. However, for the composition method, i.e., EPAVF-C, additional techniques have to be utilized and for the limitation of space, we leave it to our future works.