Convergence analysis of Fourier pseudo-spectral schemes for three-dimensional incompressible Navier-Stokes equations

1.   Introduction

In this paper we consider the three-dimensional incompressible Navier-Stokes equations (NSEs) over the domain $\Omega = (0, 1)^3$

with a periodic boundary condition (for any ${{\boldsymbol{\alpha }}} = (\alpha_1, \alpha_2, \alpha_3)$, $\alpha_j \ge 0$ being non-negative integers)

Here ${{\boldsymbol{u}}} = (u, v, w)^{T}$ is the velocity field, $p$ the pressure, and $\nu = 1/{Re}$ where $Re$ is the Reynolds number. See the related discussions in [45].

The NSEs can also be formulated in terms of the pressure Poisson equation (PPE), given by

Note that the incompressibility constraint (1.2) has been replaced by the pressure Poisson equation (1.5). Taking the divergence of (1.4) and using the pressure Poisson equation (1.5), we arrive at

Due to the periodic boundary condition and initial value: $( \nabla \, {\!\cdot\!} \, {{\boldsymbol{u}}}) \vert_{t = 0} = 0$, we observe that such a heat equation results in a trivial solution

A proof of the equivalence between (1.4)-(1.5) and the classical formulation (1.1)-(1.2) within the framework of strong solutions is straightforward. See the related works [33,34], etc.

Spectral and pseudo-spectral methods have been actively studied since the 1970s [26]. For nonlinear equations, the theoretical foundation was laid in the pioneering work [37] for steady-state spectral solutions. For time-dependent nonlinear PDEs, there have been many related works of one-dimensional conservation laws [42,43], semi-discrete viscous Burgers' equation and Navier-Stokes equations by E [21,22], etc. On the other hand, most existing works considered the spatial approximation with the time variable kept continuous. Among the existing works for the fully discrete pseudospectral method applied to nonlinear problems, a time step constraint of the form ${\Delta t} \le C h^{d/2}$ (with ${\Delta t}$ the time-step, $h$ the spatial grid size, and $d$ the dimension) has often been imposed to ensure the numerical stability, such as [4] on the one-dimensional viscous Burgers' equation.

In this work, we consider the fully discrete Fourier pseudo-spectral spatial approximation for the three-dimensional incompressible Navier-Stokes equations, combined with a variety of suitably chosen semi-implicit multistep methods in the temporal discretization. The theoretical analysis for the Fourier Galerkin spectral schemes has been reported in a few earlier works [29,30], while the analysis for the pseudo-spectral schemes turns out to be more challenging. To overcome the well-known difficulties associated with the aliasing errors appearing in the Fourier pseudo-spectral method, an aliasing error control technique, which was originally developed in a recent work [28], is applied for the incompressible Navier-Stokes equations. In addition, the Fourier pseudo-spectral spatial approximation is combined with stable time-discretization of up to fourth order which are specially tailored for the numerical stability. In more details, we adopt an explicit multi-step Adams-Bashforth approach for the fluid convection term and pressure gradient term, combined with an implicit Adams-Moulton method for the diffusion term, following similar numerical ideas in [13,28]. Moreover, the pressure variable is solved by a pressure Poisson equation, instead of being teated as a Lagrange multiplier. As a result, the computed velocity vector is proved to be divergence-free at a discrete level, so that it is $L^2$ orthogonal to the pressure gradient part, which would greatly facilitate the convergence estimate. This approach has the advantage of handling the nonlinear terms in an inexpensive way, while providing the stability associated with implicit methods. Because of the coefficient distribution of the Adams-Moulton long stencil for the diffusion part, we see that certain stability condition is satisfied, which in turn leads to an optimal rate convergence estimate. In more details, for each time-discretization, we prove that when the equation is solved by the pseudospectral method up to some fixed final time $T^*$, the method is consistent, stable and convergent (to design order) in the $H^2$ norm. In addition, the maximum norm bound of the numerical solution is automatically obtained, because of the $H^2$ error estimate and the corresponding Sobolev embedding in 3-D. This approach avoids the use of the inverse inequality in the stability analysis, so that any scaling law between the time step ${\Delta t}$ and the space grid size $h$ is not required. Instead, the numerical stability is always preserved provided that ${\Delta t}$ and $h$ are bounded by two corresponding given constants over a finite time.

The rest of the article is organized as follows. In Section 2 we review of the Fourier pseudo-spectral spatial approximation, and recall the aliasing error control technique. The first order numerical scheme is proposed in Section 3, and the detailed stability and convergence analyses are provided in Section 4. The higher order numerical schemes (in temporal accuracy) are proposed and analyzed in Section 5. The numerical results of accuracy check are presented in Section 6. Finally, the concluding remarks are given in Section 7.

2.   A review of 3-D Fourier pseudo-spectral differentiation

The Fourier series of a function $f(x, y, z) \in L^2 (\Omega)$ with $\Omega = (0, 1)^3$ is defined by

In turn, the truncated series is the projection onto the space ${\mathcal B}^N$ of trigonometric polynomials in $x$, $y$, and $z$ of degree up to $N$, given by

If $f(x, y, z)$ and all its derivatives up to $m$-th order are continuous and periodic with $\left| f^{(k)} \right| \le M$, then the truncated series converges

in which $\| \cdot \|$ denotes the $L^2$ norm.

The projection operator is one approach to a Fourier series approximation. However, sometimes we want an approximation which matches the function at a given set of points. For this purpose, an interpolation operator ${\mathcal I}_N$ is introduced. Given a uniform 3-D numerical grid with $(2N+1)$ points in each dimension, and where each point is denoted by $(x_i, y_j, z_k)$, the interpolation of the function is

where the $(2N+1)^3$ collocation coefficients $(\hat{f}_c^N)_{l, m, n}$ are given by the interpolation condition $f(x_i, y_j, z_k) = \left( {\mathcal I}_N f \right)(x_i, y_j, z_k)$. These collocation coefficients can be efficiently computed using the fast Fourier transform (FFT). Note that the interpolation coefficients depend on the number of points: increasing $N$ gives a completely different set of coefficients. Also, the collocation coefficients are not equal to the actual Fourier coefficients; the difference between them is known as the aliasing error. In general, ${\mathcal P}_N f(x, y, z) \ne {\mathcal I}_N f(x, y, z)$, and even ${\mathcal P}_N f(x_i, y_j, z_k) \ne {\mathcal I}_N f(x_i, y_j, z_k)$, except of course in the case that $f \in {\mathcal B}^N$.

The Fourier series and the formulas for its projection and interpolation allow us to easily take derivatives in the $x$, $y$, or $z$ direction by simply multiplying the appropriate Fourier coefficients $(\hat{f}_c^N)_{l, m, n}$ by $2 l \pi {\rm i}$, $2 m \pi {\rm i}$, or $2 n \pi {\rm i}$, respectively. Furthermore, we can take subsequent derivatives in the same way, so that differentiation in physical space is accomplished via multiplication in Fourier space. As long as $f$ and all is derivatives (up to $m$-th order) are continuous and periodic on $\Omega$, the convergence of the derivatives of the projection and interpolation is given by

Also see the related discussion of approximation theory [5] by Canuto and Quarteroni. Recall that the collocation coefficients $(\hat{f}_c^N)_{l, m, n}$ differ from the actual Fourier coefficients $\hat{f}_{l.m.n}$. Due to this difference, interpolation of the derivative is no longer equal to the derivative of the interpolation.

These properties of the Fourier projection and interpolation form the basis of the Fourier spectral and pseudospectral methods. The Fourier spectral method for (1.4)-(1.5) relies on the projection operator ${\mathcal P}_N$: the method is defined by the requirement that the projection of the residual onto the space ${\mathcal B}^N$ will be zero. This requirement produces a system of ordinary differential equations which are then approximated numerically. This approach is known as the Galerkin approach. An alternative to the Galerkin spectral approach is the pseudo-spectral (or collocation) approach. The Fourier pseudos-pectral method for (1.4)-(1.5) relies on the interpolation operator ${\mathcal I}_N$: the method is defined by the requirement that the interpolation of the residual onto the uniform grid will be zero. Once again, this requirement produces a system of ordinary differential equations which are then integrated numerically.

The major advantage of the collocation method is that it easier to implement, and very efficient due to the fast Fourier transform. The ease of implementation comes from the fact that the collocation approach is point-wise, which avoids many difficulties when evaluating three dimensional nonlinear terms. On the other hand, the Galerkin spectral method is much easier to analyze, because it does not suffer from aliasing errors. Many works detailing the stability and convergence analysis of spectral methods exist, . In this work, we focus on the analysis of the Fourier pseudo-spectral method applied to (1.4)-(1.5). Despite the aliasing errors that appear in the collocation method, we are able to establish its stability and convergence properties for a fixed final time.

2.1.   Discrete differentiation

Given a collocation approximation to the function $f(x, y, z)$ at the points $x_i, y_j, z_k$ described above,

we can define the discrete differentiation operators ${\mathcal D}_{Nx}$, ${\mathcal D}_{Ny}$, and ${\mathcal D}_{Nz}$ operating on the vector of grid values ${\mathbf f} = f(x_i, y_j, z_k)$. In practice, one may compute the collocation coefficients $(\hat{f_c^N)}_{l, m, n}$ via FFT, and then multiply them by the correct values (given by $2 \pi {\rm i} l$, $2 \pi {\rm i} m$, $2 \pi {\rm i} n$ in the $x$, $y$ and $z$ directions, respectively) and perform the inverse FFT. Alternatively, we can view the differentiation operators ${\mathcal D}_{Nx}$, ${\mathcal D}_{Ny}$, ${\mathcal D}_{Nz}$ as matrices, and the above process can be seen as a matrix-vector multiplication. Once again, we note that the derivative of the interpolation is not the interpolation of the derivative: ${\mathcal D}_{N} {\mathcal I}_N f \ne {\mathcal I}_N {\mathcal D}_{N} f$.

The same process is performed for the second derivatives $\partial_x^2$, $\partial_y^2$, $\partial_z^2$, where this time the collocation coefficients are multiplied by $( - 4 \pi^2 l^2)$, $( - 4 \pi^2 m^2)$ and $( - 4 \pi^2 n^2)$, respectively. Alternatively, the differentiation matrix can be applied twice, i.e. the vector ${\mathbf f}$ is multiplied by ${\mathcal D}^2_{Nx}$ for the $x$-derivative, and so on. In turn, we define the discrete Laplacian, gradient and divergence

in the point-wise sense. It is also straightforward to verify that

2.2.   Norms and inner products

Since the pseudo-spectral differentiation is taken at a point-wise level, we need to introduce a discrete $L^2$ norm and inner product to facilitate the analysis in later sections. Given any periodic grid functions ${\mathbf f}$ and ${\mathbf g}$ (over the 3-D numerical grid), we note that these are simply vectors and define the discrete $L^2$ inner product and norm

This discrete $L^2$ inner product can be computed in Fourier space rather than in physical space, with the help of Parseval's equality:

where $(\hat{f}_c^N)_{l, m, n}$, $(\hat{g}_c^N)_{l, m, n}$ are the Fourier collocation coefficients of the grid functions ${\mathbf f}$ and ${\mathbf g}$ in the expansion as in (2.1).

Furthermore, a detailed calculation shows that the following formulas of integration by parts are also valid at the discrete level [11,27,28,48]:

2.3.   A preliminary estimate in Fourier collocation space

In this section we state a lemma which will be helpful later in bounding the aliasing error for the nonlinear term. Consider a function $\varphi(x, y, z)$ which is in the space ${\mathcal B}^{2N}$. We are interested in the relationship between the function and its interpolation

where $\hat{p}^N_{l, m, n}$ are the collocation coefficients of $\varphi(x, y, z)$ for the $2N+1$ equidistant points in each dimension. Note that, because $\varphi(x, y, z) \in {\mathcal B}^{2N}$, the collocation coefficients $\hat{p}^N_{l, m, n}$ are not typically equal to the Fourier coefficients $\hat{\varphi}_{l, m, n}$. In particular, if $\varphi \in {\mathcal B}^{2N}$ stands for a nonlinear product of two terms, we see that $\varphi \ne {\mathcal I}_N \varphi$, due to the fact that it is not in ${\mathcal B}^N$. The difference between $\varphi$ and ${\mathcal I}_N \varphi$ comes from the aliasing error in the Fourier interpolation; see the detailed derivations in [41]. As a result, a control of the aliasing error associated with the nonlinear product term has always been an essential difficulty in the theoretical analysis of Fourier pseudo-spectral schemes to various nonlinear PDEs.

The following lemma will enable us to obtain an $H^m$ bound of the interpolation of the nonlinear term. In more details, the $\| \cdot \|_{H^k}$ norm of the Fourier interpolation of nonlinear product term is bounded by the same norm of the original nonlinear term, up to a fixed constant. Therefore, the aliasing error in the nonlinear product is controlled, and this lemma is usually referred as an "aliasing error control lemma". The case of $k = 0$ was proven in E's earlier work [21,22], while the case of $k \ge 1$ was analyzed in a more recent work [28]. Such an aliasing error control technique has been extensively applied to various gradient flow models, such as Cahn-Hilliard [9,17], epitaxial thin film growth [6,7,8,10,12,16,14,32,38], phase field crystal [15,46], etc.

Lemma 2.1. [28] For any $\varphi \in {\mathcal B}^{2N}$ in dimension $d$, we have

3.   First order (in time) scheme

For the incompressible NSEs (1.1), we treat the nonlinear convection term explicitly for the sake of numerical convenience, and the diffusion term implicitly to preserve an $L^2$ stability. In addition, a skew-symmetric form is taken for the nonlinear convection, and such a choice assures its $L^2$ orthogonality with the velocity vector at a discrete level. The pressure field, which is determined by the velocity field (as will be discussed later), is also updated explicitly in the same fashion of the nonlinear convection. Such an approach would bring lots of numerical convenience; and also it enforces the divergence-free condition for the numerical velocity. Consequently, the fully discrete scheme is given by:

or, in more detail,

Note that the nonlinear term is a spectral approximation to the skew-symmetric form $\frac12 \left( {{\boldsymbol{u}}} {\!\cdot\!} \nabla {{\boldsymbol{u}}} + \nabla \cdot \left( {{\boldsymbol{u}}} \otimes {{\boldsymbol{u}}} \right) \right)$ at time step $t^n$. Such a form could radically reduce the effect of aliasing error, as will be shown later.

Remark 3.1. In this work, we use a collocation Fourier spectral differentiation other than the Galerkin spectral method. It is well-known that the discrete Fourier expansion (2.1) may contain aliasing errors, while the spectral accuracy is preserved as long as the exact solution is smooth enough. Also note that in the fully discrete scheme (3.1), the gradient is computed in Fourier space by using formulas (2.2), and the multiplication of ${{\boldsymbol{u}}}^n$ and $\nabla_N {{\boldsymbol{u}}}^n$ in the nonlinear convection term is then performed in point-wise physical space. That greatly simplifies the computational efforts in the numerical simulations. Also see the related discussions in [13,27].

This approach is very different from the Galerkin approach, in which the nonlinear convection term has to be expanded in all wave length. There is no simple formula to compute these coefficients, even if it is treated explicitly. Meanwhile, in spite of aliasing errors appearing in the pseudo-spectral method, we are able to establish its local in time stability and convergence property in this work.

3.1.   Spectral solver for the pressure

Note that the pressure gradient is also updated explicitly in scheme (3.1). Such an explicit treatment avoids its coupling with the divergence-free constraint for the velocity vector. The pressure field is solved by a Fourier spectral algorithm for the pressure Poisson equation (1.5). In more detail, the following Poisson equation is formulated in Fourier space:

with the nonlinear force term given by (3.2)-(3.4):

This Poisson equation is solved using a discrete periodic boundary condition:

After the nonlinear convection terms are computed in physical space, we could apply FFT to obtain the Fourier coefficients of the force term in a collocation way. Subsequently, the pressure field can be determined by a combination of backward FFT, with the corresponding Fourier coefficients divided by the eigenvalues of the Poisson operator. An implementation of the discrete periodic boundary condition (3.7) would not effect the spectral accuracy of the numerical scheme.

3.1.1.   Divergence-free property for the velocity in spectral space

Taking a divergence of scheme (3.1) in Fourier space gives

Meanwhile, we observe that the divergence and Laplacian operators are commutative in Fourier space:

In addition, a combination of formula (2.3) and the pressure Poisson equation (3.5) shows that

As a result, its substitution into (3.8) implies that

It is a discrete version of heat equation (1.6) for the divergence field, with implicit Euler in time, Fourier spectral in space.

With the periodic boundary condition and initial incompressibility constraint

we arrive at the divergence-free property for the numerical velocity vector in Fourier spectral space:

3.1.2.   $L^2$ orthogonality between the velocity vector and pressure gradient in Fourier space

The following lemma is needed.

Lemma 3.1. For any vector field ${{\boldsymbol{v}}}$ with $\nabla_N \cdot {{\boldsymbol{v}}} \equiv 0$ and any scalar field $\phi$, we have

Proof. An application of integration by parts formula (2.5) shows that

As a result, due to the divergence-free property (3.13) for the numerical velocity vector ${{\boldsymbol{u}}}^n$, we get

3.1.3.   Further analysis

In fact, a careful analysis shows that the pressure Poisson equation (3.5) is equivalent to the Helmholtz decomposition of the nonlinear convection term in Fourier space, i.e.,

For simplicity, we denote ${{\boldsymbol{v}}}^n = {\mathcal P}_N \left( {{\boldsymbol{u}}}^n {\!\cdot\!} \nabla_N {{\boldsymbol{u}}}^n \right)$, with ${\mathcal P}_N$ a finite dimensional projection operator.

In a more general case, we have the following estimate for the Helmholtz decomposition in finite-dimensional Fourier space.

Lemma 3.2. For any vector function ${\boldsymbol{f}}$ over the 3-D grid, assume its Helmholtz decomposition in Fourier space is given by

we have the identity

Proof. Taking a discrete $L^2$ product of ${\boldsymbol{f}}$ with itself:

in which the $L^2$ orthogonality between ${{\boldsymbol{v}}}$ and $\nabla_N \phi$, as shown by Lemma 1, was used in the last step.

3.2.   Connection with the projection and gauge method

Recall the application of projection method to incompressible NSEs, with a periodic boundary condition and collocation Fourier spectral in space:

The projection method was initiated with the pioneering work of A. Chorin [18], R. Temam [44]. Also see the related works [2,23,36,39], among others.

Note that each stage is equivalent to a Poisson equation in Fourier space. In other words, two Poisson equations need to be solved at each time step, one for the intermediate velocity vector, the other for the pressure. In more detail, stage 2 is exactly a Helmholtz decomposition for ${{\boldsymbol{u}}}^{*n}$:

Moreover, its substitution into stage 1 gives

in which we utilized the commutative property between Fourier operators $\nabla_N$ and $\Delta_N$. In turn, by a comparison between the formulated numerical scheme (3.1) and the projection method (3.24), we see that the two schemes are equivalent if we set

Note that the above equation always has a unique solution in Fourier space, since all the eigenvalues associated with the left hand operator are positive. It is also observed that both sides in equation (3.25) turn out to be the gradient part in the Helmholtz decomposition of the nonlinear convection term ${{\boldsymbol{u}}}^n {\!\cdot\!} \nabla_N {{\boldsymbol{u}}}^n$ in Fourier space.

A more detailed calculation also shows its equivalence to the gauge method; see the related works [24,25,47], etc. In the gauge formulation, a gauge variable $\phi$ is introduced (instead of the pressure) and an auxiliary field is given by ${{\boldsymbol{a}}} = {{\boldsymbol{u}}} + \nabla \phi$. The first order (in time) gauge method with collocation Fourier spectral spatial differentiation is given by

Similarly, both stages lead to a Poisson equation in Fourier space. With a Helmholtz decomposition for ${{\boldsymbol{a}}}^{n+1}$ given by stage 2, we go back to stage 1 and get

Therefore, a comparison between (3.1) and (3.29) implies their equivalence if we set

which is a discrete heat equation for $\phi$. Again, both sides in equation (3.29) are the gradient part in the Helmholtz decomposition of the nonlinear convection term ${{\boldsymbol{u}}}^n {\!\cdot\!} \nabla_N {{\boldsymbol{u}}}^n$ in Fourier space.

4.   Stability and convergence analysis

Note that the numerical solution $({{\boldsymbol{u}}}, p)$ of (3.1) is evaluated at discrete grid points. Before the convergence statement of the scheme, we introduce its continuous extension in space, defined by ${{\boldsymbol{u}}}_{ {\Delta t}, h}^k = {{\boldsymbol{u}}}_N^k$, $p_{ {\Delta t}, h}^k = p_N^k$ in which $({{\boldsymbol{u}}}_N^k, p_N^k) \in {\mathcal B}^N, \forall k$, is the continuous version of the discrete grid function $({{\boldsymbol{u}}}^k, p^k)$, with the interpolation formula given by (2.1).

Our stability and convergence analysis will bound the error between this spatially continuous version of the numerical solution and the exact solution. To bound this function we will be looking at the $\| \cdot \|_{l^\infty (0, T^*; H^2)}$ and $\| \cdot \|_{l^2 (0, T^*; H^3)}$ norms. For a semi-discrete function $w$ (continuous in space and discrete in time), we define the first of these norms by

i.e., we create a discrete time-dependent function by taking the $H^2$ norm of the numerical approximation in space for each time step $t^k$, and then take the maximum of this function over all time steps $0 \leq k \leq N_k$. For the second norm, we perform a similar operation,

Theorem 4.1. For any final time $T^* > 0$, assume the exact solution to the 3-D incompressible NSEs (1.1)-(1.2) has a regularity of ${{\boldsymbol{u}}}_{{ e}} \in H^2 (0, T^*; H^{m+3})$, $p_{{ e}} \in L^\infty (0, T^*; H^{m+2})$, with $m \ge 2$. Denote ${{\boldsymbol{u}}}_{ {\Delta t}, h}$, $p_{ {\Delta t}, h}$ as the continuous (in space) extension of the fully discrete numerical solution given by scheme (3.1). As ${\Delta t}, h \to 0$, we have the following convergence result:

provided that the time step ${\Delta t}$ and the space grid size $h$ are bounded by given constants

Note that the convergence constant in (4.1), (4.2) depend on the exact solution, as well as $T^*$ and $\nu$.

The convergence analysis follows a combination of consistency analysis for the pseudo-spectral scheme (3.1) and a stability analysis for the numerical error function. In the consistency analysis, instead of a direct comparison between the numerical solution and the exact solution, we construct an approximate solution by projecting the exact solution onto ${\mathcal B}^N$. The consistency analysis shows that such an approximate solution satisfies the numerical scheme up to an $O ( {\Delta t})$ accuracy in time and a spectral accuracy in space. In the stability and convergence analysis, we first make an $H^{\frac32 + \delta}$ a-priori assumption for the numerical error function, which overcomes the difficulty in obtaining an $L^\infty$ bound for the numerical solution, with an application of 3-D Sobolev embedding. Based on this a-priori assumption, a detailed error estimate can be performed in both $L^2$ and $H^2$ norms, with the help of Lemma 1 to bound the aliasing errors associated with the nonlinear terms. Afterward, the $H^{\frac32 + \delta}$ a-priori assumption is recovered at the next time step, due to the convergence in the $H^2$ norm, so that induction (in time) can be applied.

This approach avoids an application of the inverse inequality. As a result, an unconditional numerical stability (time step vs. spatial grid size) is obtained, and no scaling law is required between ${\Delta t}$ and $h$, as compared with the classical references [4,1,20,31], among others.

4.1.   Consistency analysis

Let

The following approximation estimate is clear from our discussion from Section 2:

As a result, an application of Sobolev embedding in 3-D gives

at any fixed time. In particular, since ${{\boldsymbol{U}}}_N$ is the Fourier projection of ${{\boldsymbol{u}}}_{{ e}}$, a direct application of projection estimate indicates that

in which the 3-D Sobolev embedding of $H^4 (\Omega)$ into $W^{2, \infty} (\Omega)$, as well as the 1-D embedding of $H^1 (0, T^*)$ into $L^{\infty} (0, T^*)$ (in time), have been applied in the derivation.

In addition, it is also observed that the projected velocity vector ${{\boldsymbol{U}}}_N$ is divergence-free:

due to the fact that the exact solution ${{\boldsymbol{u}}}_{{ e}}$ is.

Looking at the time derivative of the projection operator, we observe that

in other words, $\partial_t^k {{\boldsymbol{U}}}_N$ is the truncation of $\partial_t^k {{\boldsymbol{u}}}_{{ e}}$ for any $k \ge 0$, since projection and differentiation commute. This in turn implies a spectrally accurate approximation of the corresponding temporal derivative:

The bounds on the projection error and its derivatives, namely (4.4), (4.5), (4.6) and (4.10), indicate that

in which Hölder's inequality was utilized to estimate the nonlinear term and the last step is based on the fact that (${{\boldsymbol{u}}}_{{ e}}$, $p_{{ e}}$) is the exact solution. This shows that if the solution (${{\boldsymbol{u}}}_{{ e}}$, $p_{{ e}}$) satisfies the original incompressible NSEs exactly, then its projection (${{\boldsymbol{U}}}_N$, $P_N$) will satisfy the PDE up to spectral accuracy.

Now consider the time discrete version of the equation. The temporal grid is discretized using equidistant points $t^n = n \Delta t$. We denote ${{\boldsymbol{U}}}_N^n ( {\boldsymbol{x}}) = {\mathcal P}_N {{\boldsymbol{u}}}_{{ e}} ( {\boldsymbol{x}}, t^n)$, $P_N^n ( {\boldsymbol{x}}) = {\mathcal P}_N p_{{ e}} ( {\boldsymbol{x}}, t^n)$. Recall the ${{\boldsymbol{U}}}_N^n, P_N^n \in {\mathcal B}^N$, so that it is equal to its interpolation, and its derivatives can be computed exactly by the collocation differentiation operators $\nabla_N$ and $\Delta_N$.

Moreover, for the approximate solution (${{\boldsymbol{U}}}_N^n$, $P_N^n$), we define its grid function (${{\boldsymbol{U}}}^n$, $P^n$) as its interpolation: ${{\boldsymbol{U}}}^n_{i, j, k} = {{\boldsymbol{U}}}_N^n (x_i, y_j, z_k)$, $P^n_{i, j, k} = P_N^n (x_i, y_j, z_k)$. A detailed error estimate indicates that

where $L_h^2$ denotes the discrete $L^2$ norm given by (2.4). For the temporal discretization term, we start from the following estimate

at a point-wise level (in space), in which the derivation is based on an integral form of Taylor's formula. Furthermore, by the observation (4.9), we arrive at

Consequently, a combination of (4.11) and (4.12), (4.14), (4.15), (4.16) implies the consistency of the approximate solution ${{\boldsymbol{U}}}$:

In other words, the projection of the exact solution satisfies the numerical scheme (3.1) up to an $O ( {\Delta t} + h^m)$ truncation error.

In addition, we also observe that the $H^1$ norm of $\tau$ is also bounded at the consistency order, namely

where $\tau_N^k \in {\mathcal B}^N$ is the continuous version of $\tau^k$. Such an estimate is derived based on higher order asymptotic expansions. The details are skipped for simplicity of presentation.

Another key feature of the approximation solution ${{\boldsymbol{U}}}$ is its exact divergence-free property at the spectrally discrete level:

in which the first step is based on the fact that ${{\boldsymbol{U}}}_N^n \in {\mathcal B}_N$, ${{\boldsymbol{U}}}_N^n = {\mathcal I}_N {{\boldsymbol{U}}}^n$, and the second step comes from an earlier estimate (4.8).

4.2.   Stability and convergence analysis

The point-wise numerical error grid function is given by

The difference between scheme (3.1) and the consistency (4.17) shows that

Note that the divergence-free property for the numerical velocity vector is derived from a combination of (4.19) and the same property for the numerical solution as given by (3.13). Such a property assures that the numerical error for the velocity is orthogonal to that one for the pressure gradient in the discrete $L^2$ space, thus avoids a well-known difficulty associated with the Lagrange multiplier.

To facilitate the presentation below, we denote ${{\boldsymbol{u}}}_N^n, {{\boldsymbol{e}}}_N^n, p_N^n, q_N^n \in {\mathcal B}^N$ as the continuous version of the numerical solution ${{\boldsymbol{u}}}^n, {{\boldsymbol{e}}}^n, p^n$ and $q^n$, respectively, with the interpolation formula given by (2.1).

The constructed approximate solution has a $W^{2, \infty}$ bound

which comes from the regularity of the constructed solution.

4.2.1.   An a-priori $H^{\frac32 + \delta}$ assumption

We assume a-priori that the numerical error function for the velocity has an $H^{\frac32 + \delta}$ bound at time step $t^n$:

so that the $L^\infty$ bound for the numerical solution at $t^n$ can be obtained as

It is noticed that the first estimate in (4.25) comes from the triangular inequality, while the second bound is based on the 3-D Sobolev embedding from $H^{\frac32 + \delta}$ into $L^\infty$ (for any $\delta >0$), combined with the $H^{\frac32 + \delta}$ bound derived in the first estimate.

4.2.2.   Leading order error estimate: In $L^\infty (0, T^*; L^2) \cap L^2 (0, T^*; H^1)$ norm for the velocity error

Lemma 4.1. Under the a-priori assumption (4.24), the numerical error function for the first order scheme (3.1) satisfies

Proof. Taking a discrete $L^2$ inner product of (4.21) with $2 {{\boldsymbol{e}}}^{n+1}$ yields

The time marching term and the truncation error term can be handled in a straightforward way:

in which a discrete Cauchy inequality was utilized. A discrete version of the integration by parts formula (2.5) can be applied to analyze the diffusion term:

For the pressure gradient term, we observe that the divergence-free property of ${{\boldsymbol{e}}}^{n+1}$ in the Fourier collocation space (4.25) indicates that

with the help of Lemma 2.

For the numerical error of the nonlinear convection, we see that the first term can be handled by the Cauchy inequality, with the $W^{1, \infty}$ bound of the approximate solution given by (4.23):

Similar analysis can be applied to the second nonlinear error term:

in which the $L^\infty$ a-priori bound (4.25) was used in the second step. The other two nonlinear error terms can be handled with the help of summation by parts:

As a result, a substitution of (4.28), (4.29), (4.30), (4.31), (4.32), (4.33), (4.34) and (4.35) into (4.27) gives

Summing in time gives

with $\tilde{C}_2 = 2 \left( \frac{\tilde{C}_1^2 + (C^*)^2}{\nu} + C^* + 1 \right)$. Note that we used the fact $\left\| {{\boldsymbol{e}}}^0 \right\|_2 \le C h^m$, due to the collocation spectral approximation of the initial data. An application of the discrete Gronwall inequality leads to (4.26), with $M_1 = C e^{\frac12 \tilde{C}_2 T^*}$. Also note that the truncation error estimate (4.17) was used. This in turn gives the $L^\infty (0, T^*; L^2) \cap L^2 (0, T^*; H^1)$ error estimate for the numerical scheme, under the a-priori $H^{\frac32 + \delta}$ assumption (4.24).

It is observed that the $L^2$ convergence (4.26) is based on the a-priori $H^{\frac32 + \delta}$ assumption (4.24) for the numerical solution. We need an $H^2$ error estimate to recover this assumption.

4.2.3.   Higher order error estimate: in $L^\infty (0, T^*; H^2) \cap L^2 (0, T^*; H^3)$ norm for the velocity error

Taking a discrete $L^2$ inner product of (4.21) with $2 \Delta_N^2 {{\boldsymbol{e}}}^{n+1}$ yields

The time marching term, truncation error term and the diffusion term can be analyzed as

Similar to (4.31), with the help of Lemma 2, the error for the pressure gradient is orthogonal to $\Delta_N^2 {{\boldsymbol{e}}}^{n+1}$ in the discrete $L^2$ space:

due to the fact that $\Delta_N^2 {{\boldsymbol{e}}}^{n+1}$ is also divergence-free:

Lemma 4.2. Under the a-priori assumption (4.24), we have the following estimates for the nonlinear error terms

Proof. We start with an application of summation by parts to the first term:

The remaining work is focused on the nonlinear expansion of $\nabla_N \left( {{\boldsymbol{e}}}^n {\!\cdot\!} \nabla_N {{\boldsymbol{U}}}^n \right)$. For simplicity of presentation, we only look at the first row $\nabla_N \left( {{\boldsymbol{e}}}^n {\!\cdot\!} \nabla_N U^n \right)$; the two other rows, $\nabla_N \left( {{\boldsymbol{e}}}^n {\!\cdot\!} \nabla_N V^n \right)$ and $\nabla_N \left( {{\boldsymbol{e}}}^n {\!\cdot\!} \nabla_N W^n \right)$, can be handled in the same way. Recall that ${{\boldsymbol{e}}}_N^n \in {\mathcal B}^N$ is the continuous version of the discrete grid error function ${{\boldsymbol{e}}}^n$, as in (2.1). It is obvious that

The first step is based on the fact that ${{\boldsymbol{e}}}^n$, $\nabla_N U^n$ and ${{\boldsymbol{e}}}_N^n$, $\nabla U_N^n$ have the same interpolation values. Lemma 1 was applied in the second step, due to the fact that ${{\boldsymbol{e}}}_N^n {\!\cdot\!} \nabla U_N^n \in {\mathcal B}^{2N}$. The advantage of this inequality is that the right hand side norm is measured in continuous space. Subsequently, a detailed Sobolev space expansion and applications of Hölder's inequality show that

with the regularity fact (4.23) applied in the last step. Its combination with (4.47) yields

which in turn also gives

Going back to (4.46), we get (4.42), the estimate for the first nonlinear convection error term:

The analysis of the second nonlinear convection error term also starts from the summation by parts:

Similarly, Lemma 1 has to be utilized to deal with the nonlinear expansion of $\nabla_N ( {{\boldsymbol{u}}}^n {\!\cdot\!} \nabla_N {{\boldsymbol{e}}}^n)$ in collocation Fourier space. The detailed estimates are given below.

Note that the a-priori assumption (4.24)-(4.25), combined with the following 3-D Sobolev embedding and elliptic regularity were used in the last step.

This in turn shows that

Going back to (4.48), we arrive at (4.43), the second inequality of Lemma 4.2:

The third and fourth estimates in lemma 5 can be derived in the same manner. The details are skipped for the sake of brevity.

A substitution of (4.37)-(4.40) and Lemma 5 into (4.36) indicates that

in which the $H^1$ consistency (4.18), the leading $L^2$ convergence (4.26) for the numerical scheme, combined with the elliptic regularity: $\left\| \nabla_N {{\boldsymbol{e}}}^n \right\|_2 \le C_3 \left\| \Delta_N {{\boldsymbol{e}}}^n \right\|_2$, was used in the derivation. In turn, applying the discrete Gronwall inequality leads to

with $\tilde{C}_5 = \frac{C \left( \tilde{C}_1^2 + (C^*)^2 \right)}{\nu}$, $M_2^2 = \frac{C \left( 96 (C^*)^2 \tilde{C}_3 + C \right)}{\nu} e^{\tilde{C}_5 T^*}$. Therefore, the $H^2$ error estimate for the numerical scheme is obtained from the following elliptic regularity:

As a direct consequence, the point-wise convergence of the numerical scheme is established by an application of 3-D Sobolev imbedding:

4.2.4.   Recovery of the $H^{\frac32 + \delta}$ a-priori bound (4.24)

With the help of the $H^2$ error estimate (4.49), we see that the a-priori $H^{\frac32 + \delta}$ bound (4.24) is also valid for the numerical error vector ${{\boldsymbol{e}}}^{n+1}$ at time step $t^{n+1}$ provided that

This completes the $L^\infty (0, T^*; H^2)$ convergence analysis for the velocity vector.

4.3.   Error estimate for the pressure gradient

What remains in the proof of Theorem 1 is the analysis for $\nabla_N q^n$. It is observed that equation (4.21) for the error function can be rewritten as

Since ${{\boldsymbol{e}}}^n$, ${{\boldsymbol{e}}}^{n+1}$ and $\Delta_N {{\boldsymbol{e}}}^{n+1}$ are divergence free at the discrete level, we see that $\nabla_N q^n$ is the gradient part of the left hand side of (4.50), in the Helmholtz decomposition. Furthermore, taking a divergence of (4.50) gives

In more detail, the estimate (4.23) for the constructed approximation solution, the recovered a-priori assumption (4.25) for the numerical solution, together with the Sobolev space estimates (4.49) and the established $H^2$ convergence (4.49) for the velocity vector, result in

Therefore, the $L^\infty (0, T^*;H^2)$ convergence (4.2) for the pressure is established with the help of the elliptic regularity

This finishes the proof of Theorem 4.1.

Remark 4.1. The optimal rate convergence analysis and error estimate have been established for the proposed numerical scheme (3.1) in Theorem 4.1. In more details, such a convergence analysis is based on the consistency analysis, presented in section 4.1, combined with the linearized stability analysis for the numerical error functions, provided in section 4.2. In other words, the numerical scheme is treated as a small perturbation of the exact solution and its projection, and the optimal rate error estimate (4.49) would lead to an $H^2$ upper bound of the numerical solution:

with the help of triangular inequality. This $H^2$ upper bound could be viewed as a stability estimate for the numerical solution.

In fact, this stability estimate is valid over a finite time, since the convergence estimate (4.49) is a local-in-time analysis. Because of the fully explicit treatment of the nonlinear convection term in the numerical scheme (3.1), a global-in-time estimate for the numerical solution is not directly available. Meanwhile, if the nonlinear convection is updated in a semi-implicit pattern

a global-in-time $L^2$ estimate for the numerical solution could be theoretically justified.

On the other hand, for two-dimensional (2-D) incompressible Navier-Stokes equation in vorticity-stream function formulation, a global-in-time $L^2$ and $H^1$ bound for the numerical solution is theoretically valid, even with a fully explicit treatment of the nonlinear convection term; see a few recent works [13,27], etc.

5.   Higher order (in time) schemes

To derive a second order scheme, we use a similar semi-implicit approach. As before, the nonlinear term is updated explicitly. For the convection term we use a standard second order Adams-Bashforth extrapolation formula, which involves the numerical solutions at time node points $t^n$, $t^{n-1}$, with well-known coefficients $3/2$ and $-1/2$, respectively. The same extrapolation formula can be applied to the pressure gradient term. The diffusion term is treated implicitly, using a second order Adams-Moulton interpolation. However, we do not use the standard second order formula, as this leads to difficulties in the stability analysis. Instead, we look for an Adams-Moulton interpolation such that the diffusion term is more focused on the time step $t^{n+1}$, i.e., the coefficient at time step $t^{n+1}$ dominates the sum of all other diffusion coefficients. It was discovered in [28] that the Adams-Moulton interpolation which involves the time node points $t^{n+1}$ and $t^{n-1}$ gives the corresponding coefficients as $3/4$, $1/4$, respectively, which satisfies the unconditional stability condition. Therefore, with the following notation for the skew-symmetric nonlinear term:

we formulate the fully discrete scheme:

in which the pressure field at time steps $t^n$, $t^{n-1}$ are determined by the Poisson equation as in (3.5). Since the velocity profile at these two time steps are explicit, the pressure Poisson equation is fully decoupled from the momentum equation, which greatly simplifies the computational effort.

Similar ideas can be applied to derive third and fourth order in time schemes for (1.1)-(1.2). The nonlinear convection term and the pressure gradient are updated by an explicit Adams-Bashforth extrapolation formula, with the time node points $t^n$, $t^{n-1}$, ..., $t^{n-k+1}$ involved and an order of accuracy $k$. The diffusion term is computed by an implicit Adams-Moulton interpolation with the given accuracy order in time. To ensure unconditional numerical stability for a fixed time, we have to derive an Adams-Moulton formula such that the coefficient at time step $t^{n+1}$ dominates the sum of the other diffusion coefficients. In more detail, a $k$-th order (in time) scheme takes the form of (notice that ${\mathcal NL} ({{\boldsymbol{u}}})$ is the skew-symmetric nonlinear term):

in which $B_i \mid_{i = 0}^{k-1}$ are the standard Adams-Bashforth coefficients with extrapolation points $t^n$, $t^{n-1}$, ..., $t^{n-k+1}$, $j (i) \mid_{i = 0}^{k-1}$ are a set of (distinct) indices with $j (i) \ge 0$, and $D_0$, $D_{j (i)} \mid_{i = 0}^{k-1}$ correspond to the Adams-Moulton coefficients to achieve the $k$-th order accuracy. Moreover, a necessary condition for unconditional numerical stability is given by

To derive an Adams-Moulton formula for the diffusion term, whose coefficients satisfy the condition (5.4), we require a stretched stencil. In particular, for the third order scheme, it can be shown that a stencil comprised of the node points $t^{n+1}$, $t^{n-1}$ and $t^{n-3}$ is adequate. The fully discrete scheme can be formulated as:

Similarly, the pressure field at time steps $t^n$, $t^{n-1}$, $t^{n-2}$ are determined by the Poisson equation as in (3.5); the explicit treatment of the pressure gradient decouples the pressure Poisson equation and the momentum equation.

For the fourth order scheme, we use an Adams-Moulton interpolation at node points $t^{n+1}$, $t^{n-1}$, $t^{n-5}$ and $t^{n-7}$ for the diffusion term. Combined with the Adams-Bashforth extrapolation for the nonlinear convection term, the scheme is given by

in which the $p^n$, $p^{n-1}$, $p^{n-2}$ and $p^{n-3}$ are determined by the pressure Poisson equation as in (3.5).

5.1.   Divergence-free property of the velocity in the multi-step schemes

It is a well-known numerical difficulty for incompressible fluid flow to ensure the divergence-free condition for the computed velocity vector, especially for the higher order (in time) schemes. We will prove that all the multi-step schemes proposed above satisfy this condition at the discrete level.

For simplicity, we only present the analysis for the third order scheme (5.5). The divergence-free properties of (5.2) and (5.6) can be derived in the same manner. Taking a divergence of (5.5) gives

Meanwhile, using the Poisson equation (3.5) (applied at $t^n$, $t^{n-1}$, $t^{n-2}$), we observe the Laplacian of the pressure cancels the divergence of the nonlinear term:

since $\Delta_N p^i = - \nabla_N \cdot {\mathcal NL} ({{\boldsymbol{u}}}^i)$. As a result, we obtain

With the periodic boundary condition and initial incompressibility constraint $\nabla_N \cdot {{\boldsymbol{u}}}^0 \equiv 0$, we arrive at $\nabla_N \cdot {{\boldsymbol{u}}}^n \equiv 0$, $\forall n$. the divergence-free property of the third order scheme (5.5) is proven.

5.2.   Convergence analysis for a fixed final time

Numerical stability and full order convergence (at a fixed final time) are valid for all these multi-step schemes, provided that condition (5.4) is satisfied. The result is given in the following theorem. For simplicity of presentation, we denote $K$ as the order of accuracy in time for each multi-step scheme, such as (5.2), (5.5) and (5.6).

Theorem 5.1. For any final time $T^* > 0$, assume the exact solution to the 3-D incompressible NSEs (1.1)-(1.2) has a regularity of ${{\boldsymbol{u}}}_{{ e}} \in H^K (0, T^*; H^{m+3} )$, $p_{{ e}} \in H^K (0, T^*; H^{m+2})$, with $m \ge 2$. Denote ${{\boldsymbol{u}}}_{ {\Delta t}, h}$ as the continuous (in space) extension of the fully discrete numerical solution given by the $K$-th order accurate (in time) scheme, (5.2) with $K = 2$, (5.5) with $K = 3$, (5.6) with $K = 4$. As ${\Delta t}, h \to 0$, we have the following convergence result:

provided that the time step ${\Delta t}$ and grid size $h$ are bounded by given constants

Note that the convergence constants in (5.10), (5.11) also depend on the exact solution, as well as $T^*$ and $\nu$.

We only present analysis for the third order scheme (5.5) with $K = 3$. The two other cases can be handled in the same way.

Proof. We look at the approximate solution ${{\boldsymbol{U}}}_N ( {\boldsymbol{x}}, t)$, $P_N ( {\boldsymbol{x}}, t)$, given by (4.3), and denote $({{\boldsymbol{U}}}_{i, j, k}^n, P_{i, j, k}^n)$ as the numerical value of $({{\boldsymbol{U}}}_N, P_N)$ at grid point $(x_i, y_j, z_k, t^n)$. Again, (4.11) is satisfied. The following truncation error estimates can be derived using a high order Taylor expansion in time:

That in turn leads to the third order consistency of the approximate solution ${{\boldsymbol{U}}}$:

with $\left\| \tau \right\|_{l^2 (0, T^*; L_h^2)} \le C \left( {\Delta t}^3 + h^m \right)$. The $H^1$ bound of $\tau$ can also be derived:

Subsequently, with the numerical error function denoted by (4.20) (while $({{\boldsymbol{e}}}_N^n, q_N^n)$ represents the continuous version of $({{\boldsymbol{e}}}^n, q^n)$), the difference between (5.5) and (5.12) gives

Similarly, the constructed solution ${{\boldsymbol{U}}}$ satisfies the $W^{2, \infty}$ regularity condition (4.23). To deal with a multi-step method, the $H^{\frac32 + \delta}$ a-priori bound is assumed for the numerical error function at all previous time steps:

for any $1 \le k \le n$, with $\tilde{C}_0, \tilde{C}_1$ given by (4.24)-(4.25).

Taking a discrete $L^2$ inner product of (5.13) with $2 {{\boldsymbol{e}}}^{n+1}$ yields

Since the diffusion coefficient at $t^{n+1}$ dominates the other ones, the viscosity term can be analyzed by

The estimates for the nonlinear convection terms are given below. The details are skipped.

Going back to (5.15) results in

Consequently, summing in time, applying the discrete Gronwall inequality and using a simple fact that $\frac34 - \frac{2}{24} - \frac{1}{12} - \frac{11}{24} = \frac{1}{8} > 0$, we get a fixed time $O( {\Delta t}^3 + h^m)$ convergence for the third order scheme (5.5) in $L^\infty (0, T^* ; L^2)$ norm:

under the a-priori $H^{\frac32 +\delta}$ assumption (5.14).

Similar convergence analysis in $L^\infty (0, T^*; H^2) \cap L^2 (0, T^*; H^3)$ can be performed by taking a discrete $L^2$ inner product of (5.13) with $2 \Delta_N^2 {{\boldsymbol{e}}}^{n+1}$. The details are skipped.

hence,

In turn, the a-priori $H^{\frac32 + \delta}$ bound (5.14) for ${{\boldsymbol{e}}}^{n+1}$ at time step $t^{n+1}$ is assured, provided that

Finally, the error estimate for the pressure can be carried out in the same way as in section 4.3. The technical details are skipped for the sake of brevity.

Remark 5.1. For the second order in time method, treating the diffusion term with a standard second order Adams-Moulton formula leads to a Crank-Nicolson scheme. In the third order case, we observe that a direct application of the Adams-Moulton formula at the nodes $t^{n+1}$, $t^n$ and $t^{n-1}$ (corresponding to $j (i) = i$ in the general form (5.3)) does not give a formula with the stated stability property. For example, a "naive" combination of 3rd order Adams-Bashforth for the nonlinear convection and Adams-Moulton for the diffusion term results in the following scheme

which violates the stability condition (5.4). Numerical experiments also showed that this method is unstable in time. This case highlights the need to choose an appropriate time-discretization to couple with the pseudospectral method.

Remark 5.2. There have been extensive numerical simulation works of high-order multi-step, semi-implicit schemes applied to nonlinear time-dependent PDEs. Among them, it is worth mentioning the $AB/BDI$ (Adams-Bashforth/Implicit Backward Differentiation) schemes, introduced by Crouzeix [19], which update the nonlinear convection term using the standard Adams-Bashforth extrapolation and the time marching term with implicit backward differentiation for the diffusion term. In particular, the third order scheme of this family has been applied to the incompressible Navier-Stokes equations with various spatial approximations; see the relevant works [3,35]. However, only the linear stability analysis is covered in these existing works; see the relevant discussions in [40]. To the authors' knowledge, this article is the first to provide the stability and convergence analysis for a third order scheme applied to incompressible Navier-Stokes equations.

6.   The numerical results

In this section we perform a numerical accuracy check for the proposed numerical schemes (3.1), (5.2), (5.5), (5.6), in the first, second, third and fourth order temporal accuracy orders, respectively. The computational domain is chosen as $\Omega = (0, 1)^2$, and the exact profiles for ${{\boldsymbol{u}}}_exac$ and $p_{{ e}}$ are set to be

The viscosity parameter is taken as $\nu = 0.5$. To make $({{\boldsymbol{u}}}_e, p_e)$ satisfy the original PDE system (1.1)-(1.2), we have to add an artificial, time-dependent forcing term. The proposed schemes (3.1), (5.2), (5.5), (5.6), can be very efficiently implemented using the FFT solver.

In the accuracy check for the temporal accuracy, we fix the spatial resolution as $N = 256$ (with $h = \frac{1}{256}$), so that the spatial numerical error is negligible, because of the spectral accuracy order in space. The final time is set as $T = 1$. Naturally, a sequence of time step sizes are taken as ${\Delta t} = \frac T{N_T}$, with $N_T = 100:100:1000$. The expected temporal numerical accuracy assumption $e = C {\Delta t}^k$, with $k$ the temporal accuracy order, indicates that $\ln |e| = \ln (C T^k) - k \ln N_T$, so that we plot $\ln |e|$ vs. $\ln N_T$ to demonstrate the temporal convergence order. The fitted line displayed in Figure 1 shows approximate slopes of -0.9993, - 2.0018, -2.9972, -4.01522, respectively. These numerical results have verified perfect temporal convergence orders of the proposed numerical schemes, for the physical variable of velocity component $u$, in both the discrete $\ell^2$ and $\ell^\infty$ norms. The results for the velocity component $v$ and the pressure variable $p$ are very similar.

7.   Concluding remarks

In this paper we develop stable time-stepping methods of order up to four for the numerical solution of the three dimensional incompressible Navier-Stokes equation, with Fourier pseudo-spectral spatial approximation, and present stability and convergence analyses for these fully discrete methods.

In the first order (in time) scheme, a semi-implicit approach is applied, which updates the nonlinear convection term and the pressure gradient term explicitly, and treats the diffusion term implicitly. The pressure variable is solved by a pressure Poisson equation, instead of being teated as a Lagrange multiplier. As a result, the computed velocity vector is proved to be divergence-free at a discrete level, so that it is $L^2$ orthogonal to the pressure gradient part, which would greatly facilitate the convergence estimate. Since the pseu-dospectral method is evaluated at the interpolation grid points, a discrete $L^2$ and $H^m$ inner product is needed to carry out the analysis, in which an estimate of the nonlinear convection term in the $L^2$ and $H^1$ norms has to be derived. An aliasing error control technique, originally developed in [28] to deal with viscous Burgers' equation, is applied. With the help of such a technique, all the energy estimates in the Sobolev norms can be performed almost in the same way as in the Galerkin approach, so that stability and convergence over a fixed finite time is demonstrated. In particular, the $L^\infty$ bound of the numerical solution is the key difficulty in the stability and convergence analysis. We provide an $H^2$ error estimate so that the $L^\infty$ norm of the numerical solution is automatically bounded, using a 3-D Sobolev imbedding. This approach avoids a time-step restriction in terms of the spatial grid size.

Similar ideas could be applied to develop higher order in time schemes using a multi-step approach, with numerical stability preserved. For the sake of numerical efficiency, the semi-implicit pattern is kept, with a standard Adams-Bashforth extrapolation formula (with the given accuracy) applied to the nonlinear term and the pressure gradient term, while the diffusion term is treated using an implicit Adams-Moulton interpolation formula on certain time node points. On the other hand, the numerical stability requires specialized stretched stencil in the Adams-Moulton approximation. We present these multi-step schemes with accuracy up to fourth order and show that, coupled with the Fourier pseudo-spectral spatial approximation, these multi-step numerical schemes are stable and convergent. Some numerical results have also been presented, which demonstrate the perfect convergence rates of the proposed schemes.

Acknowledgments

This work is supported in part by NSF DMS-2012669.