Efficient composition conservative schemes for the Zakharov equations

1.   Introduction

The Zakharov system, introduced by Vladimir Zakharov in 1972 (see ^[16]), is a collection of coupled nonlinear partial differential equations. Its purpose is to simulate the interaction between high-frequency Langmuir waves and low-frequency ion-acoustic waves within a plasma. Comprising a complex equation for the electric field envelope and a real equation for the ion density fluctuation, this system effectively captures crucial plasma dynamics. Such dynamics include modulational instability and wave collapse phenomena. In the realm of plasma physics, the Zakharov system holds great significance. It details the energy transfer process between different wave modes, which ultimately results in the emergence of localized structures such as solitons and vortices. Moreover, its applications extend beyond plasma. In nonlinear optics, the Zakharov system functions as a fundamental model for comprehending the propagation of intense laser beams in nonlinear media. Additionally, it has relevance in fluid dynamics and Bose-Einstein condensates. In these fields, similar mathematical structures regulate wave interactions (references ^[6,21,27] attest to this). In this paper, we focus on the following Zakharov system.

where the complex-valued function $E(\mathbf{x}, t)$ represents the slowly varying envelope of a rapidly oscillating electric field. The real-valued function $N(\mathbf{x}, t)$ describes the low-frequency variations in ion density. The given initial functions $E_{0}(\mathbf{x})$, $N_{0}(\mathbf{x})$, and $N_{1}(\mathbf{x})$ satisfy the periodic boundary conditions on $\partial\Omega$. For further physical background, please refer to ^[12,17]. Additionally, the Zakharov system (1.1) can conserve the mass ^[29]

and the Hamiltonian energy ^[30]

where $\nabla v = N_t.$ In recent years, a wealth of numerical studies has been conducted on both classical and generalized Zakharov systems. These works encompass various methods such as the time-splitting technique ^[1,14], the scalar auxiliary variable approach ^[23], finite difference schemes ^[3,26], and discontinuous Galerkin methods ^[25], among others. Recently, there has been a significant interest in developing accurate and efficient numerical methods for various Zakharov systems. For example, Bao and Sun ^[1] developed a time-splitting spectral method for the ZS, which effectively captures the intricate dynamics of plasma interactions. Jin and Sulem ^[14] proposed numerical schemes that preserve the multi-scale nature of Zakharov systems, ensuring stability and accuracy over long-time simulations. Additionally, Xia and Shu ^[25] introduced discontinuous Galerkin methods for the quantum Zakharov systems, offering high-order accuracy and flexibility in handling complex geometries. These advancements contribute to a deeper understanding and more precise simulation of the nonlinear phenomena described by the Zakharov system. However, these mentioned schemes cannot conserve the mass (1.2) and the energy (1.3) of the system (1.1).

Structure-preserving algorithms are numerical methods specifically designed to maintain the intrinsic properties and structures of the original differential equations, such as symmetries, invariants, and conservation laws ^[5,7,11]. By preserving these features, the algorithms ensure that the numerical solutions exhibit behavior consistent with the underlying physical or mathematical models ^[13,24]. Advantages of structure-preserving algorithms include enhanced stability, improved long-term accuracy, and the ability to capture essential qualitative dynamics of the system ^[8,15]. They are particularly beneficial in simulations where conserving quantities like energy or mass is crucial for obtaining reliable and physically meaningful results ^[28]. Recently, there has been considerable progress in developing structure-preserving numerical schemes for Zakharov systems. These schemes are essential for accurately simulating the long-term dynamics of plasma waves, as they conserve the total energy of the system under discretization. By employing techniques such as discrete variational principles and implicit time-stepping methods, researchers have designed algorithms that maintain the inherent energy conservation property of the Zakharov system ^[4,20]. This not only enhances computational stability but also ensures that numerical solutions faithfully replicate the physical behavior of nonlinear wave interactions in plasma physics.

The averaged vector field (AVF) method can exactly preserve the energy of Hamiltonian systems ^[19,22] and has been applied to construct conservative schemes for energy-conserving partial differential equations. However, these schemes are fully implicit, requiring the development of nonlinear iteration procedures. Moreover, since all unknown variables are averaged simultaneously in the numerical schemes, the number of iterations for the nonlinear solver can be large. To reduce iterative costs at each time step while maintaining the desired energy-preserving property, Cai et al. ^[2] developed partitioned averaged vector field (PAVF) methods. These methods conserve energy and have been used to construct conservative schemes for Hamiltonian ODEs. Compared with the AVF method, the schemes resulting from the PAVF method are much simpler in their explicit expressions. They effectively reduce the original fully implicit schemes to semi-implicit or linearly implicit ones, significantly improving computational efficiency. However, as far as we known, the PAVF method has not yet been applied to construct conservative schemes for the Zakharov equations.

In this paper, we develop a structure-preserving scheme by using the PAVF method in time and the Fourier pseudo-spectral method in space for the system. The constructed numerical scheme offers several significant advantages. It strictly conserves both energy and mass, preserving the inherent physical properties of the Zakharov system and ensuring reliable long-term simulations. By simplifying the fully implicit approach to a semi-implicit or linearly implicit scheme, it substantially reduces computational costs and enhances efficiency. The scheme maintains good stability even with larger time steps due to the implicit handling of stiff terms. Additionally, it is easier to implement since it avoids solving complex nonlinear systems. Incorporating spectral methods provides high accuracy, enabling precise capture of solution details and dynamic behaviors. These features make the scheme highly applicable to various nonlinear partial differential equations with similar structures.

The rest of the paper is organized as follows. In Section 2, we reformulate the original Zakharov equation as a Hamiltonian system and provid a brief introduction to the PAVF method. Section 3 presents the development of the semi-discrete structure-preserving numerical schemes by using the Fourier pseudo-spectral method. In Section 4, we derive the fully discrete scheme by using the PAVF method to approximate the semi-discrete system in time and prove the conservative of the schemes. Section 5 provides numerical experiments to demonstrate the efficiency and accuracy of our methods, comparing them with existing approaches. Finally, Section 6 concludes the paper with a summary of the main results and suggestions for future research.

2.   The Hamilton form of the Zakharov system and PAVF method

Constructing numerical schemes using the PAVF method is based on the Hamiltonian form of the system. First, we present the Hamiltonian form of the equations.

2.1.   The Hamilton form of the Zakharov system

To this end, we decompose the complex-valued function $E$ into its real and imaginary parts, namely

where $q$ and $p$ are real-valued functions. Based on this, we rewrite the original system of Eq (1.1) as the following real-valued first-order system:

Denote $z = (q, p, N, v)^T$, based on the variational principle, and we obtain

where $\dfrac{\delta \mathcal{H}}{\delta u}$ denotes the variational derivative of $\mathcal{H}$ with respect to $u$, $u = p, q, N, v$. Then, system (2.1) can be written into the following energy-conserving system:

where

is the skew-adjoint operator, and the Hamiltonian energy functional is

and the mass functional is

2.2.   The PAVF method

We consider the Hamiltonian system:

where $w \in \mathbb{R}^k$, $S_k$ is a skew-symmetric $k \times k$ matrix with even $k$, and $H(w)$ is a sufficiently differentiable Hamiltonian function.

The second-order AVF method for system (2.6) is defined as

where $\tau$ is the time step.

Setting $k = 2d$ and partitioning $w$ as $w = (y, z)^{T}$ with $y, z \in \mathbb{R}^{d}$, the system (2.6) becomes

and the Hamiltonian $H(y, z)$ remains conserved along any continuous flow. The PAVF method for system (2.8) is given by

Denote the PAVF method (2.9) by $\Psi_{\tau}$. Its adjoint method $\Psi_{\tau}^{*}$ is defined as

While noting that both scheme $\Psi_{\tau}$ and its adjoint scheme $\Psi_{\tau}^{*}$ are semi-implicit first-order one-step methods with high efficiency, they only possess first-order accuracy and cannot meet the high-accuracy requirements of numerical simulations. Therefore, based on composition techniques, we combine these two methods to obtain a numerical scheme with second-order accuracy. By combining $\Psi_{\tau}$ and its adjoint $\Psi_{\tau}^{*}$, we define the PAVF composition (PAVF-C) method:

Remark 2.1. The PAVF method (2.9) is a first-order one-step method. But the PAVFC method $\Upsilon_{\tau}$ in (2.11) has second-order methods that exactly conserve the energy of the Hamiltonian system (2.6).

3.   Spatial discretization

In this section, we apply the Fourier pseudo-spectral method in space for the Zakharov system, resulting in a semi-discrete conservative scheme. Note that other space discretization methods can also be used, such as finite difference, finite element, or finite volume methods.

3.1.   Fourier spectral discretization in space

This subsection applies the Fourier pseudo-spectral method to the Zakharov system. Note that other space discretization methods can also be used, such as finite difference, finite element, or finite volume methods. Without losing generality, we consider the Zakharov system in the 2D case with the periodic boundary condition and $\Omega = (-L, L)^2$. For a positive even integer $M$, we partition the domain uniformly with mesh sizes $h = h_x = h_y = 2L/M$. Let the spatial grid points

with $x_i = -{L}+ih, y_j = -{L}+jh$. We denote

For ${\mathit{\boldsymbol{u}}}, {\mathit{\boldsymbol{v}}}\in \mathcal{U}_h$, we define the corresponding discrete inner product and norms as follows:

Let $X_{i}(x)$ be the interpolation basis functions given by

where $\mu = \pi/L$ and

Then, we can define the two-dimensional interpolation space

and the interpolation operator $I_{M}:C(\overline{\Omega})\rightarrow S_{M}$

The corresponding second-order spectral differential matrices for the $x$- and $y$-directions are uniformly calculated by

Owing to the circulant property of these differential matrices, one can utilize the fast Fourier transform (FFT) to accelerate the computation of matrix-vector multiplication with $D_2$. In fact, we have the decomposition that

where ${F}_{M}$ denotes the matrix of the discrete Fourier transform, ${F}_{M}^H$ is the conjugate transpose matrix of ${F}_{M}$ and ${F}_{M}^H = {F}_{M}^{-1}$ ^[10]. The diagonal matrix $\Lambda$ corresponds to the eigenvalues of $D_1$ and $D_2$ with the elements given by

Moreover, the approximation of the two-dimensional Laplace operator by the pseudo-spectral method yields the second-order differential matrix $\mathcal{D}: = I_M\otimes{D}+{D}\otimes{I}_M$, which by (3.1) also admits a diagonal decomposition

Therefore, the practical computation associated with the differential matrix will be efficiently carried out.

3.2.   Construction of the semi-discrete conservative schemes

Applying the pseudo-spectral method to discretize the spatial derivative of system (2.1), the semi-discrete scheme is as follows:

Similarly to the continuous case, the above system can be reformulated in Hamiltonian form:

where the discrete Hamiltonian is defined as

and $S$ is a skew-symmetric matrix given by

Then, the conservation laws of the semi-discrete scheme are presented in the following theorem.

Theorem 3.1. The scheme (3.4) satisfies the semi-discrete mass and energy conservation laws

and

where

and $H({\mathit{\boldsymbol{z}}})$ is given by (3.6).

Proof. Based on the definition of the mass, we have

In addition to mass conservation, the Hamiltonian $H({\mathit{\boldsymbol{z}}})$ of the system is conserved. To verify this, we compute the time derivative of $H({\mathit{\boldsymbol{z}}})$ can obtain

where $S$ is a skew-symmetric matrix. The final equality holds because for any skew-symmetric matrix $S$ and vector $\nabla H({\mathit{\boldsymbol{z}}})$

This property ensures that the Hamiltonian $H({\mathit{\boldsymbol{z}}})$ remains constant over time, signifying energy conservation within the system. □

4.   Fully discrete conservative schemes

Building upon the analysis of the semi-discrete scheme (3.4), this subsection aims to develop a class of fully discrete conservative schemes for both temporal and spatial discretization.

4.1.   Conservative PAVF scheme

First, we apply the PAVF method for system (3.4) can derive

which can be further integrated as

where $\mathit{\boldsymbol{z}}^{n+\frac{1}{2}} = \frac{{\mathit{\boldsymbol{z}}}^{n}+{\mathit{\boldsymbol{z}}}^{n+1}}{2}$, ${\mathit{\boldsymbol{z}}} = {\mathit{\boldsymbol{p}}}, {\mathit{\boldsymbol{q}}}$ etc. The proposed PAVF schemes (4.5)–(4.8) are denoted as $\Psi^{}_{\tau}$. Similarly, we can obtain the adjoint PAVF method $\Psi^{*}_{\tau}$, namely

Theorem 4.1. The proposed schemes (4.5)–(4.8) and their adjoint schemes (4.9)–(4.12) possess the following discrete total mass and energy conservation laws:

where the mass is defined as

and the energy is given by

Proof. For the scheme $\Psi^{}_{\tau}$, taking the inner product of (4.7) and (4.8) with $\mathit{\boldsymbol{p}}^{n+1} + \mathit{\boldsymbol{p}}^n$ and $\mathit{\boldsymbol{q}}^{n+1} + \mathit{\boldsymbol{q}}^n$, respectively, we obtain

and

Adding them together, we deduce

which implies that

Then, we take the discrete inner products of (4.5)–(4.8) with $\delta_t \mathit{\boldsymbol{v}}^n$, $\delta_t \mathit{\boldsymbol{N}}^n$, $-\delta_t \mathit{\boldsymbol{q}}^n$, and $\delta_t \mathit{\boldsymbol{p}}^n$, respectively, and sum them together to obtain the energy conservation law

Similarly, we can prove the adjoint schemes (4.9)–(4.12) also can conserve the mass and energy in fully discrete sense. This concludes the proof. □

4.2.   Composition conservative scheme

It is noteworthy that the PAVF scheme only achieves first-order accuracy in time, and this subsection will construct second-order numerical schemes that achieve energy and mass conservation based on the composition technique.

PAVFC scheme. For given $\mathit{\boldsymbol{z}}^n = \left({\mathit{\boldsymbol{p}}}^n, {\mathit{\boldsymbol{q}}}^n, {\mathit{\boldsymbol{N}}}^n, {\mathit{\boldsymbol{v}}}^n \right)^T$, we update $\mathit{\boldsymbol{z}}^{n+1} = \left({\mathit{\boldsymbol{p}}}^{n+1}, {\mathit{\boldsymbol{q}}}^{n+1}, {\mathit{\boldsymbol{N}}}^{n+1}, {\mathit{\boldsymbol{v}}}^{n+1} \right)^T$ via the following composition steps:

with

and $\Psi_{\frac{\tau}{2}}$, $\Psi^{*}_{\frac{\tau}{2}}$ represent replacing $\tau$ in the PAVF scheme and its adjoint scheme with $\frac{\tau}{2}$, which can be further expressed as

Based on the proof of Theorem 4.1, we can immediately obtain the following theorem.

Theorem 4.2. The proposed PAVFC scheme can conserve the following energy and mass conservation law, namely

To facilitate a comparison between the standard AVF method and the PAVFC methods, we discretize the semi-discrete system in time using the original second-order AVF method, thereby obtaining the following AVF scheme for the equation.

AVF scheme.

Remark 4.1. Theorem 4.2 demonstrates that the PAVFC scheme simultaneously preserves both the energy and mass conservation of the system, whereas the traditional AVF scheme only maintains energy conservation and fails to preserve mass. We will validate these properties through numerical examples.

5.   Numerical results

This section aims to evaluate the validity, conservation properties, and efficiency of the newly proposed schemes for solving system (1.1) under periodic boundary conditions. Denote ${\mathit{\boldsymbol{E}}}^n$ and ${\mathit{\boldsymbol{N}}}^n$ be the numerical solutions obtained by the proposed scheme with mesh size $h$ and time step $\tau$ for approximating the exact solutions $E(\cdot, t_n)$ and $N(\cdot, t_n)$. To quantify the numerical methods, we define the error functions as

Meanwhile, we also define the relative residuals on the mass and energy as

Example 5.1. Consider that the one-dimensional Zakharov system admits a solitary-wave solution given by ^[9]

The initial conditions are defined as

The computations are performed within the interval $\Omega = [-128,128]$. This domain is chosen to ensure adequate coverage of the solution's behavior and to capture the dynamics accurately.

We first propose the numerical errors and convergence orders for four schemes in Table 1. This table presents the numerical errors and convergence orders for four different numerical schemes: the AVF, PAVF, PAVFC, and CN schemes when solving a system with a grid size of $N = 1024$ and a final time $T = 1$. Here, the CN scheme represents using the Crank-Nicolson method to construct the conservative schemes for the Zakharov system ^[26]. The errors and convergence rates are shown for varying time step sizes $\tau$, ranging from 0.01 to 0.01/8. The AVF, PAVFC, and CN schemes demonstrate second-order accuracy, as indicated by the consistent convergence order of approximately 2.00 across all time steps. On the other hand, the PAVF scheme shows first-order accuracy, with the convergence order remaining close to 1.00 for each time step. These results confirm that the AVF, PAVFC, and CN schemes provide higher accuracy than the PAVF scheme. Figure 1 compares the relative mass errors (RM) and relative energy errors (RH) for four different numerical schemes: AVF, PAVF, PAVFC, and CN, over time. The graphs are plotted with a spatial step size $h = 0.5$ and time step $\tau = 0.001$. In the left panel, the relative mass errors (RM) are shown. The PAVFC and AVF schemes demonstrate stable and small errors, with the AVF scheme maintaining higher stability in mass conservation. The PAVF scheme exhibits larger mass errors, indicating less accurate mass conservation compared to the others. The CN scheme also shows relatively stable mass error behavior but it is slightly larger than the AVF and PAVFC schemes. The right panel shows the relative energy errors (RH). The AVF scheme exhibits larger energy errors, which increase over time. Both the PAVFC and CN schemes show superior energy conservation, maintaining smaller and more stable errors. The PAVF scheme, while improving in mass conservation, shows higher energy errors than PAVFC and CN. The PAVFC scheme, in particular, stands out with high accuracy while preserving both energy and mass conservation, making it an effective method for numerical simulations.

Example 5.2. Let us study two-dimensional system (1.1) with the following initial conditions:

When quantum effects are neglected ($\varepsilon = 0$) and $2 \nu - \frac{\nu^2}{\mu} > 4 \pi$, the solution to the classical ZS model will blow up at some point in time, as demonstrated in ^[14,18]. This issue has been used in ^[1,14] to showcase the numerical behavior of their methods. Additionally, considering quantum effects, we employ the PAVFC scheme to analyze the conservation properties and dynamics of system (1.1). For simplicity, we set $\varepsilon = 0$ and $\mu = \nu = 20$. The problem is confined to the domain $\Omega = [-8, 8] \times [-8, 8]$.

Figure 2 shows the relative errors of mass (RM) and energy (RH) over time for the PAVFC scheme. The computations were performed with grid spacings $h_x = h_y = 0.5$ and a time step size $\tau = 0.01$. The red line represents the relative mass error (RM), which stays relatively small throughout the simulation. Meanwhile, the blue line shows the relative energy error (RH), which also remains at a low level but fluctuates slightly more compared to the mass error. Overall, the figure demonstrates that the PAVFC scheme maintains both mass and energy conservation effectively, with minimal errors in both quantities over the time interval depicted.

Figure 3 displays the evolution of the energy density $|E(x, y, t)|^2$ and ion density fluctuation $N(x, y, t)$ for the Zakharov system at different times. The left column $(a, c, e)$ shows the energy density $|{\mathit{\boldsymbol{E}}}|^2$, while the right column $(b, d, f)$ presents the ion density fluctuation ${\mathit{\boldsymbol{N}}}$. At $t = 0.1$ (plots $a$ and $b$), the initial energy density and ion fluctuation exhibit a concentrated, localized structure at the center of the domain. Over time, both quantities become more focused and sharper. By $t = 1.0$ (plots $c$ and $d$), the energy density $|{\mathit{\boldsymbol{E}}}|^2$ has condensed further, showing a sharp peak, while the ion density fluctuation ${\mathit{\boldsymbol{N}}}$ also intensifies, reflecting a more focused region of fluctuation. At $t = 1.16$ (plots $e$ and $f$), the energy density $|{\mathit{\boldsymbol{E}}}|^2$ and ion fluctuation ${\mathit{\boldsymbol{N}}}$ have reached their most concentrated state, with significantly higher peak values compared to earlier times. The figures effectively illustrate the dynamics of energy and ion fluctuations in the Zakharov system as they evolve over time, with both quantities becoming increasingly localized and peaking as time progresses.

6.   Conclusions

In this paper, we propose a novel class of structure-preserving methods based on the partitioned averaged vector field approach and Fourier pseudo-spectral method to solve the Zakharov system. These schemes are designed to accurately preserve the system's key properties, including mass and Hamiltonian energy conservation. By leveraging the PAVF framework and combining it with the composition technique, we derive second-order accuracy and a stability conservative scheme. Numerical experiments are conducted to validate the efficiency and accuracy of the proposed schemes. Furthermore, our method demonstrates second-order convergence of the numerical solution to the Zakharov system. The results confirm that the PAVF schemes offers a robust and reliable tool for preserving the structure of the Zakharov system.

Author contributions

Yan Zhang: Writing-original draft, Writing- review & editing; Yuyang Gao: Methodology, Writing-original draft, Writing-review & editing; Yayun Fu: Methodology; Xiaopeng Yue: Methodology, Writing-review & editing. All authors have read and approved the final version of the manuscript for publication.

Use of Generative-AI tools declaration

The authors declare that they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The research is supported by the key scientific research projects of colleges and universities in Henan Province (24A170030).

Conflict of interest

The authors declare there is no conflict of interest.