An efficient Fourier spectral method and error analysis for the fourth order problem with periodic boundary conditions and variable coefficients

1.   Introduction

Though periodic boundary conditions do not belong to three traditional boundary conditions (Dirichlet, Neumann, and Robin boundary conditions), which are commonly used in mathematical physics, they still can be found in the research of scientific and engineering problems, such as the interaction between solutions of the nonlinear Schrödinger equation ^[1] or KdV equation ^[2], isotropic uniform turbulence problem ^[3], and so on. In addition, under the polar, cylindrical, and spherical coordinates, we note that they are also periodic ^[4,5] in the $\theta$ direction. Thus, it's obvious that physical models with periodic boundary conditions also have significant application. As a model, we first consider the following two dimensional fourth-order problem with periodic boundary conditions and variable coefficients:

where $\alpha$ is a nonnegative bounded periodic function, $\beta$ is a positive bounded function, $\mathbf{x} = ({x_1}, {x_2})$, $L_{x_1} = x_{1 R}-x_{1 L}$, $L_{x_2} = x_{2 R}-x_{2 L}$, $\Omega = (x_{1 L}, x_{1 R})\times(x_{2 L}, x_{2 R}).$

The fourth-order problems can be found in the applications to thin beams and plates ^[6,7]. Besides, many complex nonlinear problems also need to solve a fourth order problem repeatedly ^{[8,9,10,11,12,13]}. In the past decades, there have been many existing results for the theoretical analysis and numerical research of the fourth-order problems, mainly including various finite element methods ^{[14,15,16,17]}, spectral methods and some high-order numerical methods ^{[18,19,20,21,22,23,24,25,26,27]}. However, to the best of our knowledge, there are few report on the fourth-order problems with periodic boundary conditions and variable coefficients ^[28]. As aforementioned, periodic boundary conditions have significant applications in some science and engineering ^[29,30]. Thus, it is meaningful to construct an efficient and high-order numerical scheme for the fourth-order problems with periodic boundary conditions and variable coefficients.

The aim of this paper is to propose an efficient algorithm based on the Fourier spectral-Galerkin approximation for the fourth-order elliptic equation with periodic boundary conditions and variable coefficients. First, by using the Lax-Milgram theorem, we prove the existence and uniqueness of weak solution and its approximate solution. Then we define a high-dimensional $L^2$ projection operator and prove its approximation properties. Combined with Céa lemma, we further prove the error estimate of the approximate solution. In addition, from the Fourier basis function expansion and the properties of the tensor, we establish the equivalent matrix form based on tensor product for the discrete scheme. Finally, some numerical experiments are carried out to demonstrate the efficiency of the algorithm and correctness of the theoretical analysis.

The rest of this paper is organized as follows. In next section, we derive the weak form and associated discrete scheme. We give the error estimation of approximate solutions in section 3. In section 4, we present an efficient implementation of the algorithm. In section 5, we extend the algorithm to a three-dimensional case. In section 6, we carry out some numerical experiments. Finally, we make some concluding remarks in section 7.

2.   Weak form and discrete scheme

We shall derive the weak form and discrete scheme associated with problem (1.1)–(1.3). Denote by $H^m(\Omega)$ the usual m-order Sobolev space, $\|\cdot\|_m$ and $|\cdot|_m$ denote the norm and semi-norm in $H^m(\Omega)$, respectively. In particular, we have

with the following inner product and norm

where $\bar{\varphi}$ is the complex conjugate of $\varphi$. Define

with the following inner product, norm and semi-norm:

where $D^{\mathit{\boldsymbol{\alpha}}} = \frac{\partial^{|\mathit{\boldsymbol{\alpha}}|}}{\partial {x_1}^{{\alpha}_{1}}\partial {x_2}^{{\alpha}_{2}}}, \; \boldsymbol{\alpha} = (\alpha_1, \alpha_2), \; |\mathit{\boldsymbol{\alpha}}| = \alpha_1+\alpha_2$. We further denote by $H_{p}^{m}(\Omega)$ the subspace of $H^m(\Omega)$, which consists of functions with derivatives of order up to $m-1$ being $2\pi$-periodic.

Then a weak form of problem (1.1)–(1.3) is: Find $\psi\in H_{p}^{2}(\Omega)$, such that

where

Define an approximation space of $H_{p}^{2}(\Omega)$ as follows:

Then the corresponding discrete scheme of the weak form (2.1) is$:$ Find $\psi_{ M}\in X_M(\Omega)$, such that

3.   Error estimation of the approximation solution

In this section, we shall first prove the existence and uniqueness of weak solution and its approximate solution, and then further prove the error estimate between them.

3.1.   Existence and uniqueness

For the sake of brevity, we denote by $a\lesssim b$ that $a\leq cb$, where $c$ is a positive constant. Without loss of generality, we shall confine our discussion to the following assumptions:

Lemma 1. For any $\psi, \varphi\in H^{2}_{p}(\Omega)$, the following equalities hold:

Proof. Using the integration by parts, we have

From (1.2) and (1.3), we derive that

Then we have

It follows from the above equality that

We can obtain the following equalities in the same way

Lemma 2. $For\; any\; \psi\in H^{2}_{p}(\Omega)$, the following inequalities hold:

Proof. We derive from integration by parts that

Similarly, we can obtain

This finishes our proof.

Lemma 3. Let $\alpha(\mathbf{x})$, $\beta(\mathbf{x})\in L^\infty(\Omega)$ satisfy the assumptions (3.1) and (3.2). Then $a(\psi, \varphi)$ is a continuous and coercive bilinear form in $H_{p}^2(\Omega)\times H_{p}^2(\Omega)$, i.e.,

where $Q^* = \max\{1, \alpha^*, \beta^*\}$, $Q_* = \frac{1}{2}\min\{1, \beta_*\}$.

Proof. Employing Lemma 1, we obtain that

Then, using Cauchy-Schwarz inequality, we can derive that

On the other hand, using Lemma 2, we can derive that

This finishes our proof.

Lemma 4. If $f(\mathbf{x})\in L^2(\Omega)$, then $F(\varphi)$ is a bounded linear functions on $H^{2}_{p}(\Omega)$, i.e.,

Proof. In light of definition of $F(\varphi)$ and Cauchy-Schwarz inequality, we have

The proof is completed.

From Lemma 3, Lemma 4 and Lax-Milgram theorem, we have following theorem:

Theorem 1. If $f(\mathbf{x})\in L^2(\Omega)$, then problems (2.1) and (2.2) have unique solutions $\psi(\mathbf{x})$ and $\psi_{ M}(\mathbf{x})$, respectively.

3.2.   Error estimation

Theorem 2. Let $\psi(\mathbf{x})$ and $\psi_{ M}(\mathbf{x})$ be the solutions of the variational form (2.1) and discrete scheme (2.2), respectively. Then it holds that

Proof. We obtain from (2.1) and (2.2) that

Then we have

Form Lemma 3 and (3.4), we arrive at

which is equivalent to the following form

From (3.5) and the arbitrariness of $\varphi_{ M}$, the desired result follows.

Let $\Pi_{ M}:L^2(\Omega)\rightarrow {X_M}(\Omega)$ be a $L^2$-orthogonal projection:

Theorem 3. For any $\psi(\mathbf{x})\in H^m_p(\Omega)$ and $0\leq\mu\leq m$, there is a constant $C$ such that the following inequality holds:

Proof. We first derive that

For any $|\mathit{\boldsymbol{\alpha}}|:0\leqslant|\mathit{\boldsymbol{\alpha}}|\leqslant\mu\leq m$, taking $\alpha_1\leq m_1, \; \alpha_2\leq m-m_1$, we have

by making a summation for $|\mathit{\boldsymbol{\alpha}}|$ from $0$ to $\mu$, we can obtain the expected results.

Theorem 4. Let $\psi_{ M}(\mathbf{x})$ be the approximation solution of $\psi(\mathbf{x})$. If $\psi(\mathbf{x}) \in H^{m}_{p}(\Omega)$, the following inequality holds

Proof. According to Theorem 2, we have

We derive form Theorem 3 that

The proof is completed.

4.   Efficient implementation of the algorithm

In this section, we will describe the implementation process of the algorithm in detail, and give a brief pseudo code. To solve (2.2) by Fourier spectral method, we shall look for

Let

We denote by $\overline{\Psi}$ a column vectors with $(2M+1)^2$ elements, which consist of $2M+1$ columns of $\Psi$. Let $\overline{\varphi}_{ M}(\mathbf{x}) = e^{-ik{x_1}}e^{-il{x_2}}, (|k|, |l| = 0, 1, \cdots, M)$, then we have

where

$S(k, :)$ indicates the k-th row of the matrix $S$, $M(k, :)$, $O(k, :)$ and $G(k, :)$ are similar to $S(k, :)$.

$\otimes$ represents the tensor product of matrix, i.e. $M\otimes S = (m_{ kt}S)^{2M+1}_{|k|, |t| = 0}$.

where

Then the equivalent matrix form based on tensor product for the discrete scheme (2.2) is as follows:

where

Note that when $\alpha, \beta$ are constants, we know from the orthogonal property of Fourier basis functions that the stiffness matrix and mass matrix in (4.2) are all sparse, so we can solve (4.2) efficiently. However, for general variable coefficients $\alpha, \beta$, the stiffness matrix and mass matrix are usually full. In that case, we can use the preconditioned iteration method or Schur-complement approach, i.e., block Gaussian elimination to solve (4.2).

5.   Extension to three-dimensional case

In this section, we shall extend our algorithm to three-dimensional case. As a model, we consider the following three dimensional fourth-order problem with periodic boundary conditions:

where $\alpha$ and $\beta$ are constant coefficients, $L_{x_1} = x_{1 R}-x_{1 L}$, $L_{x_2} = x_{2 R}-x_{2 L}$, $L_{x_3} = x_{3 R}-x_{3 L}$, $\mathbf{x} = ({x_1}, {x_2}, {x_3})$, $\Omega = (x_{1 L}, x_{1 R})\times(x_{2 L}, x_{2 R})\times(x_{3 L}, x_{3 R}).$

Similar to two-dimensional case, we can derive the weak form and discrete scheme for the three-dimensional case. Define a Sobolev space:

Then a weak form of (5.1)–(5.4) is to find $\psi\in H_{p}^{2}(\Omega)$, such that

Define an approximation space:

Then the corresponding discrete scheme for the weak form (5.5) is to find $\psi_{ M}\in X_M(\Omega)$, such that

We shall derive the equivalent matrix form based on tensor product for the discrete scheme (5.6). Let

We denote by $\tilde{\Psi}^j$ a column vectors with $(2M+1)^2$ elements consisting of $2M+1$ columns of $\Psi^j$. Let $\Psi = (\tilde{\Psi}^{ -M}, \tilde{\Psi}^{ -M+1}, \cdots, \tilde{\Psi}^{ M})$, and denote by $\tilde{\Psi}$ a column vectors with $(2M+1)^3$ elements consisting of $2M+1$ columns of $\Psi$. Taking $\overline{\varphi}_{ M}({\mathbf{x}}) = e^{-ik{x_1}}e^{-il{x_2}}e^{-ip{x_3}}, (|k|, |l|, |p| = 0, 1, \cdots, M)$, then we have

Then the equivalent matrix form based on tensor product for the discrete scheme (5.6) is as follows:

where

6.   Numerical experiments

In this section, we shall perform some numerical experiments to confirm the correctness of theoretical analysis and the effectiveness of our algorithm. The programs are compiled and operated in MATLAB 2018b.

6.1.   Two dimensional case

Example 1. We take $\alpha = 1, \; \beta = 10$ and choose the exact solution $\psi = \sin 4{x_1}\sin 8{x_2}$. Then $f$ can be obtained by plugging $\psi$ into the Eq (1.1). We shall solve (1.1)–(1.3) by using the algorithm proposed in section 4. We list in Table 1 the errors between the exact solution and the approximate solution under $H^2$ norm, $H^2$ seminorm and $L^2$ norm respectively for different $M$. In addition, we also present their comparison figures and absolute error figures for different $M$ in Figures 1 and 2.

We observe from Table 1 that the approximate solution $\psi_{ M} (\mathbf{x})$ reaches about $10^{-12}$ accuracy when $M\geq8$. Besides, we also see from Figures 1 and 2 that the approximation solution converges to the exact solution.

Remark 1. Though the proof of well-posedness of weak solution requires $\alpha$ to be a nonnegative bounded periodic function and $\beta$ to be a positive bounded function, our algorithm is still valid for some large and negative $\alpha$ and $\beta$, and the corresponding numerical results are listed in Table 2.

Example 2. We take $\alpha = \sin({x_1}+{x_2})+2, \; \beta = e^{\cos({x_1}+{x_2})}$, and choose the exact solution $\psi = e^{\sin({x_1}+{x_2})}$. We list in Table 3 the errors between the exact solution and the approximate solution under $H^2$ norm, $H^2$ seminorm and $L^2$ norm respectively for different $M$. Similarly, we also present their comparison figures and absolute error figures for different $M$ in Figures 3 and 4. In order to further show the spectral accuracy of our algorithm, we present in Figure 5 error figures between exact solution and approximation solutions under $L^2$ and $H^2$ norms for different $M$.

We observe from Table 3 that the approximate solution $\psi_{ M} (\mathbf{x})$ reach about $10^{-11}$ accuracy when $M\geq12$. We see from Figures 3–5 that the approximation solution exponentially converges to the exact solution.

6.2.   Three dimensional case

Example 3. We take $\alpha = 1, \; \beta = 1$ and choose the exact solution $\psi = \cos 3x_{1}\cos 4x_{2}\cos 5x_{3}$. We list in Table 4 the errors between the exact solution and the approximate solution under $H^2$ norm, $H^2$ seminorm and $L^2$ norm respectively for different $M$.

We observe from Table 4 that the approximate solution reach about $10^{-13}$ accuracy when $M\geq8$. That is to say, even in the three-dimensional case, our algorithm still has spectral accuracy.

Example 4. We take $\alpha = 1, \; \beta = 1$ and choose the exact solution $\psi = e^{\cos x_{1}+\cos x_{2}+\cos x_{3}}$. We list in Table 5 the errors between the exact solution and the approximate solution under $H^2$ norm, $H^2$ seminorm and $L^2$ norm respectively for different $M$.

We observe from Table 5 that the approximate solution reach about $10^{-10}$ accuracy when $M\geq12$. Again, our algorithm has spectral accuracy.

7.   Conclusions

We have developed an efficient Fourier spectral-Galerkin method to solve the fourth-order elliptic equation with periodic boundary conditions and variable coefficients. Firstly, we prove the error estimations between the weak solutions and approximation solutions. Then we derive the equivalent matrix form based on tensor product for the discrete scheme. Numerical experiments validate the theoretical analysis and algorithm. Besides, the method proposed in this paper can be extended to some more complex linear and nonlinear equations, such as fourth-order parabolic equation ^[31], Cahn-Hiliard equation, Gross Pitaevskii equation, which is our future research goal.

Acknowledgments

The authors would like to thank the anonymous referees for their valuable comments and suggestions, which led to considerable improvement of the article. The research is supported by National Natural Science Foundation of China (Grant No. 12061023), Guizhou Normal University Academic Young Talent Foundation(Qian Teacher Young Talent [2021] A04), and Research Foundation for Scientific Scholars of Moutai Institute (Grant number [2022] 140).

Conflict of interest

The authors declare that they have no competing interests.