Research article

Influence of cholesterol on human calcitonin channel formation. Possible role of sterol as molecular chaperone

  • Received: 31 October 2018 Accepted: 24 January 2019 Published: 01 February 2019
  • The interplay between lipids and embedded proteins in plasma membrane is complex. Membrane proteins affect the stretching or disorder of lipid chains, transbilayer movement and lateral organization of lipids, thus influencing biological processes such as fusion or fission. Membrane lipids can regulate some protein functions by modulating their structure and organization. Cholesterol is a lipid of cell membranes that has been intensively investigated and found to be associated with some membrane proteins and to play an important role in diseases. Human calcitonin (hCt), an amyloid-forming peptide, is a small peptide hormone. The oligomerization and fibrillation processes of hCt can be modulated by different factors such as pH, solvent, peptide concentration, and chaperones. In this work, we investigated the role of cholesterol in hCt incorporation and channel formation in planar lipid membranes made up of palmitoyl-oleoyl-phosphatidylcholine in which no channel activity had been found. The results obtained in this study indicate that cholesterol promotes hCt incorporation and channel formation in planar lipid membranes, suggesting a possible role of sterol as a lipid target for hCt.

    Citation: Daniela Meleleo, Cesare Sblano. Influence of cholesterol on human calcitonin channel formation. Possible role of sterol as molecular chaperone[J]. AIMS Biophysics, 2019, 6(1): 23-38. doi: 10.3934/biophy.2019.1.23

    Related Papers:

    [1] Thuy Hien Nguyen, Catherine C. Moore, Preston B. Moore, Zhiwei Liu . Molecular dynamics study of homo-oligomeric ion channels: Structures of the surrounding lipids and dynamics of water movement. AIMS Biophysics, 2018, 5(1): 50-76. doi: 10.3934/biophy.2018.1.50
    [2] Carl S. Helrich . Studies of cholesterol structures in phospholipid bilayers. AIMS Biophysics, 2017, 4(3): 415-437. doi: 10.3934/biophy.2017.3.415
    [3] Chia-Wen Wang, Meng-Han Lin, Wolfgang B. Fischer . Cholesterol affected dynamics of lipids in tailor-made vesicles by ArcVes software during multi micro second coarse grained molecular dynamics simulations. AIMS Biophysics, 2023, 10(4): 482-502. doi: 10.3934/biophy.2023027
    [4] Iva Valkova, Petar Todorov, Victoria Vitkova . VV-hemorphin-5 association to lipid bilayers and alterations of membrane bending rigidity. AIMS Biophysics, 2022, 9(4): 294-307. doi: 10.3934/biophy.2022025
    [5] Jacques Fantini, Francisco J. Barrantes . How membrane lipids control the 3D structure and function of receptors. AIMS Biophysics, 2018, 5(1): 22-35. doi: 10.3934/biophy.2018.1.22
    [6] Rianne Bartelds, Jonathan Barnoud, Arnold J. Boersma, Siewert J. Marrink, Bert Poolman . Lipid phase separation in the presence of hydrocarbons in giant unilamellar vesicles. AIMS Biophysics, 2017, 4(4): 528-542. doi: 10.3934/biophy.2017.4.528
    [7] Sudarat Tharad, Chontida Tangsongcharoen, Panadda Boonserm, José L. Toca-Herrera, Kanokporn Srisucharitpanit . Local conformations affect the histidine tag-Ni2+ binding affinity of BinA and BinB proteins. AIMS Biophysics, 2020, 7(3): 133-143. doi: 10.3934/biophy.2020011
    [8] Shan-Yang Lin . Salmon calcitonin: conformational changes and stabilizer effects. AIMS Biophysics, 2015, 2(4): 695-723. doi: 10.3934/biophy.2015.4.695
    [9] Mathieu F. M. Cellier . Evolutionary analysis of Slc11 mechanism of proton-coupled metal-ion transmembrane import. AIMS Biophysics, 2016, 3(2): 286-318. doi: 10.3934/biophy.2016.2.286
    [10] Nicola Gaetano Gatta, Gaetano Cammarota, Vittorio Gentile . Possible roles of transglutaminases in molecular mechanisms responsible for human neurodegenerative diseases. AIMS Biophysics, 2016, 3(4): 529-545. doi: 10.3934/biophy.2016.4.529
  • The interplay between lipids and embedded proteins in plasma membrane is complex. Membrane proteins affect the stretching or disorder of lipid chains, transbilayer movement and lateral organization of lipids, thus influencing biological processes such as fusion or fission. Membrane lipids can regulate some protein functions by modulating their structure and organization. Cholesterol is a lipid of cell membranes that has been intensively investigated and found to be associated with some membrane proteins and to play an important role in diseases. Human calcitonin (hCt), an amyloid-forming peptide, is a small peptide hormone. The oligomerization and fibrillation processes of hCt can be modulated by different factors such as pH, solvent, peptide concentration, and chaperones. In this work, we investigated the role of cholesterol in hCt incorporation and channel formation in planar lipid membranes made up of palmitoyl-oleoyl-phosphatidylcholine in which no channel activity had been found. The results obtained in this study indicate that cholesterol promotes hCt incorporation and channel formation in planar lipid membranes, suggesting a possible role of sterol as a lipid target for hCt.


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

    $ tu+uu+p=νΔu,
    $
    (1.1)
    $ u=0,
    $
    (1.2)

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

    $ Dαu(0,y,z)=Dαu(1,y,z),Dαu(x,0,z)=Dαu(x,1,z),Dαu(x,y,0)=Dαu(x,y,1).
    $
    (1.3)

    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

    $ tu+uu+p=νΔu,
    $
    (1.4)
    $ Δp=(u)u,
    $
    (1.5)

    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

    $ t(u)=Δ(u).
    $
    (1.6)

    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

    $ u0.
    $
    (1.7)

    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.

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

    $ f(x,y,z)=l,m,n=ˆfl,m,ne2πi(lx+my+nz),withˆfl,m,n=Ωf(x,y,z)e2πi(lx+my+nz)dxdydz.
    $

    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

    $ {\mathcal P}_N f(x, y, z) = \sum\limits_{l, m, n = -N}^{N} \hat{f}_{l, m, n} {\rm e}^ {2 \pi {\rm i} (l x + m y + n z) }. $

    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

    $ \| f(x, y, z) - {\mathcal P}_N f(x, y, z) \| \leq C M N^{-m}, $

    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

    $ \left({\mathcal I}_N f \right)(x, y, z) = \sum\limits_{l, m, n = -N}^{N} (\hat{f}_c^N)_{l, m, n} {\rm e}^ {2 \pi {\rm i} (l x + m y + n z) }, $

    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

    $ \|\partial^k f(x, y, z) - \partial^k {\mathcal P}_N f(x, y, z) \| \leq \| f^{(m)} \| N^{k-m} , \quad \mbox{for} \, \, \, 0 \le k \le m , $
    $ \|\partial^k f(x, y, z) - \partial^k {\mathcal I}_N f(x, y, z) \| \leq \| f \|_{H^m} N^{k-m} , \quad \mbox{for} \, \, \, 0 \le k \le m , \, m > \frac{d}{2} . $

    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.

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

    $ f(xi,yj,zk)=(INf)i,j,k=Nl,m,n=N(ˆfNc)l,m,ne2πi(lxi+myj+nzk),
    $
    (2.1)

    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

    $ ΔNf=(D2Nx+D2Ny+D2Nz)f,Nf=(DNxfDNyfDNzf),N(fgh)=DNxf+DNyg+DNzh,
    $
    (2.2)

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

    $ NNf=ΔNf.
    $
    (2.3)

    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

    $ f2=f,f,withf,g=(12N+1)32Ni,j,k=0fi,j,kgi,j,k.
    $
    (2.4)

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

    $ f,g=Nl,m,n=N(ˆfNc)l,m,n¯(ˆgNc)l,m,n=Nl,m,n=N(ˆgNc)l,m,n¯(ˆfNc)l,m,n,
    $

    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]:

    $ f,N(g1g2g3)=Nf,(g1g2g3),f,ΔNg=Nf,Ng.
    $
    (2.5)

    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

    $ {\mathcal I}_N \varphi(x, y, z) = \sum\limits_{l, m, n = -N}^N \hat{p}^N_{l, m, n} {\rm e}^{2 \pi {\rm i} ( l x + m y + n z)}, $

    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

    $ INφHk(2)dφHk.
    $
    (2.6)

    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:

    $ un+1unΔt+12(unNun+N(unun))+Npn=νΔNun+1,
    $
    (3.1)

    or, in more detail,

    $ un+1unΔt+12(unNun+N(unun))+DNxpn=νΔNun+1,
    $
    (3.2)
    $ vn+1vnΔt+12(unNvn+N(vnun))+DNypn=νΔNvn+1,
    $
    (3.3)
    $ wn+1wnΔt+12(unNwn+N(wnun))+DNzpn=νΔNwn+1.
    $
    (3.4)

    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.

    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:

    $ ΔNpn=12N(unNun+N(unun))=N(Fn1Fn2Fn3),
    $
    (3.5)

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

    $ Fn1=12(unDNxun+vnDNyun+wnDNzun+DNx(unun)+DNy(unvn)+DNz(unwn)),Fn2=12(unDNxvn+vnDNyvn+wnDNzvn+DNx(vnun)+DNy(vnvn)+DNz(vnwn)),Fn3=12(unDNxwn+vnDNywn+wnDNzwn+DNx(wnun)+DNy(wnvn)+DNz(wnwn)).
    $
    (3.6)

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

    $ pn(0,j,k)=pn(2N+1,j,k),pn(i,0,k)=pn(i,2N+1,k),pn(i,j,0)=pn(i,j,2N+1).
    $
    (3.7)

    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.

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

    $ (Nun+1)(Nun)Δt+N(unNun)+NNpn=νNΔNun+1.
    $
    (3.8)

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

    $ NΔNun+1=ΔN(Nun+1).
    $
    (3.9)

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

    $ N(unNun)+NNpn=N(unNun)+ΔNpn0.
    $
    (3.10)

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

    $ (Nun+1)(Nun)Δt=νΔN(Nun+1).
    $
    (3.11)

    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

    $ Nu00,
    $
    (3.12)

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

    $ Nun0,n0.
    $
    (3.13)

    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

    $ Nϕ,v=0.
    $
    (3.14)

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

    $ Nϕ,v=ϕ,Nv=0.
    $
    (3.15)

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

    $ Npn,un=0.
    $
    (3.16)

    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.,

    $ unNun=vn+Npn,withNvn0.
    $
    (3.17)

    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

    $ f=v+Nϕ,withv=PNf,Nv0,
    $
    (3.18)

    we have the identity

    $ f22=v22+Nϕ22.
    $
    (3.19)

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

    $ f22=f,f=v+Nϕ,v+Nϕ=v22+Nϕ22+2Nϕ,v=v22+Nϕ22,
    $
    (3.20)

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

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

    $ Stage 1:ununΔt+12(unNun+N(unun))=νΔNun,
    $
    (3.21)
    $ Stage 2:{un+1unΔt+NΦn+1=0,Nun+1=0.
    $
    (3.22)

    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} $:

    $ un=un+1+N(ΔtΦn+1),withNun+1=0.
    $
    (3.23)

    Moreover, its substitution into stage 1 gives

    $ un+1unΔt+12(unNun+N(unun))+N(Φn+1νΔtΔNΦn+1)=νΔNun+1,
    $
    (3.24)

    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

    $ (IνΔtΔN)Φn+1=pn.
    $
    (3.25)

    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

    $ Stage 1:an+1anΔt+12(unNun+N(unun))=νΔNan+1,
    $
    (3.26)
    $ Stage 2:{un+1an+1+Nϕn+1=0,Nun+1=0.
    $
    (3.27)

    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

    $ un+1unΔt+12(unNun+N(unun))+N(ϕn+1ϕnΔtνΔNϕn+1)=νΔNun+1.
    $
    (3.28)

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

    $ ϕn+1ϕnΔtνΔNϕn+1=pn,
    $
    (3.29)

    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.

    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

    $ \| w \|_{l^\infty (0, T^*; H^2)} = \max\limits_{0 \leq k \leq N_k} \| w^k ( {\boldsymbol{x}})\|_{H^2}, \quad N_k = \left[ \frac{T^*}{ {\Delta t}} \right] , $

    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,

    $ \| w \|_{l^2 (0, T^*; H^3)} = \sqrt{ {\Delta t} \sum\limits_{k = 0}^{N_k} \| w^k (x) \|_{H^3}^2 }. $

    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:

    $ uΔt,huel(0,T;H2)+νuΔt,huel2(0,T;H3)C(Δt+hm),
    $
    (4.1)
    $ pΔt,hpel(0,T;H2)C(Δt+hm),
    $
    (4.2)

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

    $ ΔtL1(T,ν),hL2(T,ν).
    $

    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.

    Let

    $ UN(x,t)=PNue(x,t),PN(x,t)=PNpe(x,t).
    $
    (4.3)

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

    $ UNueL(0,T;Hr)ChmueL(0,T;Hm+r),forr0,
    $
    (4.4)
    $ PNpeL(0,T;Hr)ChmpeL(0,T;Hm+r),forr0.
    $
    (4.5)

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

    $ UNueL(Ω)CUNueH2ChmueL(0,T;Hm+2),
    $
    (4.6)

    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

    $ UNL(0,T;Hr)ueL(0,T;Hr),forr0,so thatUNL(0,T;W2,)CUNL(0,T;H4)CueL(0,T;H4)CueH1(0,T;H4),
    $
    (4.7)

    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:

    $ uN=(PNu)=PN(u)0,
    $
    (4.8)

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

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

    $ ktkUN(x,t)=ktkPNue(x,t)=PNkue(x,t)tk;
    $
    (4.9)

    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:

    $ kt(UNue)ChmktueHm,for0k2.
    $
    (4.10)

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

    $ tUN+12(UNUN+(UNUN))PNνΔUN=tue+ueuepeνΔue+τ1=τ1,withτ1L2Chm,
    $
    (4.11)

    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

    $ UnNUn=UnNUnN,ΔNUn+1=ΔUn+1N,NPn=PnNsinceUnN=INUn,PnN=INPn
    $
    (4.12)
    $ N(UnUn)=(UnNUnN)+τn2,τn2L2hChmUnN2Hm+3,
    $
    (4.13)
    $ ΔUn+1N=ΔUnN+τn3,τ3l2(0,T;L2h)CΔtUNH1(0,T;H2)CΔt,
    $
    (4.14)

    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

    $ Un+1NUnNΔt=tUN(,tn)+τn4(),τ4()l2(0,T)CΔtUN()H2(0,T),
    $
    (4.15)

    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

    $ τ4l2(0,T;L2h(Ω))CΔt2tUNL2(0,T;H2(Ω))CΔt2tueL2(0,T;H2(Ω))CΔt.
    $
    (4.16)

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

    $ Un+1UnΔt+12(UnNUn+N(UnUn))NPnνΔNUn+1=τn,withτ=τ1+τ2+τ3+τ4,τl2(0,T;L2h)C(Δt+hm).
    $
    (4.17)

    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

    $ τl2(0,T;H1h):=ΔtNkk=0τkN(x)2H1C(Δt+hm),
    $
    (4.18)

    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:

    $ NUn=UnN0,n,
    $
    (4.19)

    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).

    The point-wise numerical error grid function is given by

    $ eni,j,k=Uni,j,kuni,j,k,qni,j,k=Pni,j,kpni,j,k.
    $
    (4.20)

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

    $ en+1enΔt+12(enNUn+unNen+N(en(Un+un)))+Nqn=νΔNen+1+τn,
    $
    (4.21)
    $ Nen+10.
    $
    (4.22)

    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

    $ UNL(0,T;W2,)C,i.e.UnNLC,UnNLC,(UnN)LC,n0,
    $
    (4.23)

    which comes from the regularity of the constructed solution.

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

    $ enNH32+δ1,with enN=INen,
    $
    (4.24)

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

    $ unNH32+δ=UnNenNH32+δUnH32+δ+enNH3/2+δC+1:=˜C0,ununNLCδunNH32+δCδ˜C0:=˜C1.
    $
    (4.25)

    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.

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

    $ en+12+νΔtn+1k=1Nek22M1(Δt+hm).
    $
    (4.26)

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

    $ en+1en,2en+12νΔtΔNen+1,en+1+2ΔtNqn,en+1=ΔtenNUn,en+1ΔtunNen,en+1ΔtN(enUn),en+1ΔtN(enun),en+1+2Δtτn,en+1.
    $
    (4.27)

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

    $ en+1en,2en+1=en+122en22+en+1en22,
    $
    (4.28)
    $ 2τn,en+1τn22+en+122,
    $
    (4.29)

    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:

    $ 2ΔNen+1,en+1=2Nen+1,Nen+1=2Nen+122.
    $
    (4.30)

    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

    $ Nqn,en+1=0,
    $
    (4.31)

    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):

    $ enNUn,en+1NUnen2en+12Cen2en+1212Cen22+12Cen+122.
    $
    (4.32)

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

    $ unNen,en+1unNen2en+12˜C1Nen2en+1212νNen22+˜C212νen+122,
    $
    (4.33)

    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:

    $ N(enUn),en+1=enUn,Nen+114νNen+122+(C)2νen22,
    $
    (4.34)
    $ N(enun),en+1=enun,Nen+114νNen+122+˜C21νen22.
    $
    (4.35)

    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

    $ en+122en22+en+1en22+32νΔtNen+12234νΔtNen22+(˜C21+(C)2ν+C+1)Δt(en+122+en22)+Δtτn22.
    $

    Summing in time gives

    $ en+122+nk=0ek+1ek22+34νΔtn+1k=1Nek22˜C2Δtn+1k=0ek22+Δtnk=0τn22+Ch2m,
    $

    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.

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

    $ 2en+1en,Δ2Nen+12νΔtΔNen+1,Δ2Nen+1+2ΔtNqn,Δ2Nen+1=Δt(enNUn,Δ2Nen+1+unNen,Δ2Nen+1)ΔtN(en(Un+un)),Δ2Nen+1+τn,Δ2Nen+1.
    $
    (4.36)

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

    $ 2en+1en,Δ2Nen+1=ΔNen+122ΔNen22+ΔN(en+1en)22,2τn,Δ2Nen+1=2Nτn,NΔNen+1
    $
    (4.37)
    $ 4νNτn22+14νNΔNen+122,
    $
    (4.38)
    $ 2ΔNen+1,Δ2Nen+1=2NΔNen+1,NΔNen+1=2NΔNen+122.
    $
    (4.39)

    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:

    $ Nqn,Δ2Nen+1=0,
    $
    (4.40)

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

    $ N(Δ2Nen+1)=ΔN(Nen+1)0.
    $
    (4.41)

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

    $ enNUn,Δ2Nen+114νNΔNen+122+48(C)2ν(Nen22+en22),
    $
    (4.42)
    $ unNen,Δ2Nen+114νNΔNen+122+C˜C21νΔNen22,
    $
    (4.43)
    $ N(enUn),Δ2Nen+114νNΔNen+122+48(C)2ν(Nen22+en22),
    $
    (4.44)
    $ N(enun),Δ2Nen+114νNΔNen+122+C˜C21νΔNen22.
    $
    (4.45)

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

    $ enNUn,Δ2Nen+1=N(enNUn),NΔNen+1.
    $
    (4.46)

    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

    $ N(enNUn)2=(IN(enNUnN))L222(enNUnN)L2.
    $
    (4.47)

    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

    $ (enNUnN)L2=UnNenN+enN(UnN)L2UnNenNL2+enN(UnN)L2UnNLenNL2+enNL2(UnN)LC(enNL2+enNL2)=C(Nen2+en2),
    $

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

    $ N(enNUn)222C(Nen2+en2),
    $

    which in turn also gives

    $ N(enNUn)226C(Nen2+en2).
    $

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

    $ enNUn,Δ2Nen+126C(Nen2+en2)NΔNen+1214νNΔNen+122+48(C)2ν(Nen22+en22).
    $

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

    $ unNen,Δ2Nen+1=N(unNen),NΔNen+1.
    $
    (4.48)

    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.

    $ N(unNen)2=(IN(unNenN))L222(unNenN)L2,(unNenN)L2=enNunN+unN(enN)L2enNunNL2+unN(enN)L2unNL3enNL6+unNL(enN)L2C˜C1ΔenNL2=C˜C1ΔNen2.
    $

    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.

    $ unNL3CunNH32C˜C1,enNL6CΔenNL2,unNL˜C1,(enN)L2CΔenNL2.
    $

    This in turn shows that

    $ N(unNen)222(unNenN)L2C˜C1ΔNen2.
    $

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

    $ unNen,Δ2Nen+1C˜C1ΔNen2NΔNen+1212νNΔNen+122+C˜C21νΔNen22.
    $

    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

    $ ΔNen+122ΔNen22+ΔN(en+1en)22+34νΔtNΔNen+122C˜C21ΔtνΔNen22+96(C)2Δtν(Nen22+en22)+4ΔtνNτn22,C(˜C21+(C)2)ΔtνΔNen22+(96(C)2˜C3+C)Δtν(Δt+hm)2,
    $

    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

    $ ΔNen+122+nk=0ΔN(ek+1ek)22+34νΔtn+1k=1NΔNek22C(96(C)2˜C3+C)νe˜C5T(Δt+hm)2M22(Δt+hm)2,
    $

    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:

    $ en+1NH2C(en+1NL2+Δen+1NL2)˜C6(Δt+hm),˜C6=C(M1+M2).
    $
    (4.49)

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

    $ en+1en+1NLCen+1NH2C(M1+M2)(Δt+hm).
    $

    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

    $ Δt(˜C6)1,h(˜C6)1m,with ˜C6 dependent on ν and exp(T).
    $

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

    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

    $ 12(enNUn+unNen+N(en(Un+un)))τn=(en+1enΔtνΔNen+1)Nqn.
    $
    (4.50)

    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

    $ ΔNqn=12N(enNUn+unNen+N(en(Un+un)))+Nτn.
    $
    (4.51)

    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

    $ ΔNqn212N(enNUn+unNen+N(en(Un+un)))2+Nτn26(enNUnN+unNenN+(enN(UnN+unN)))L2+τnH1C(UnNW2,enNH2+unNLenNH2+unNL3enNL6+unNH2enNL)+τnH1C(C+˜C1)enNH2+τnH1˜C7(Δt+hm).
    $
    (4.52)

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

    $ qnNH2CΔqnN2=CΔNqn2C˜C7(Δt+hm).
    $
    (4.53)

    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:

    $ ukNH2UkNekNH2UkNH2+ekNH2C+˜C6(Δt+hm)C+1,
    $
    (4.54)

    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

    $ \frac12 \left( {{\boldsymbol{u}}}^n {\!\cdot\!} \nabla_N {{\boldsymbol{u}}}^{n+1} + \nabla_N \cdot ( {{\boldsymbol{u}}}^n \otimes {{\boldsymbol{u}}}^{n+1} ) \right) , $

    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.

    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:

    $ NL(uk)=12(ukNuk+N(ukuk)),
    $
    (5.1)

    we formulate the fully discrete scheme:

    $ un+1unΔt+32NL(un)12NL(un1)+N(32pn12pn1)=νΔN(34un+1+14un1),
    $
    (5.2)

    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):

    $ un+1unΔt+k1i=0Bi(NL(uni)+Npni)=νΔN(D0un+1+k1i=0Dj(i)unj(i)).
    $
    (5.3)

    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

    $ D0>k1i=0|Dj(i)|.
    $
    (5.4)

    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:

    $ un+1unΔt+2312NL(un)43NL(un1)+512NL(un2)+N(2312pn43pn1+512pn2)=νΔN(23un+1+512un1112un3).
    $
    (5.5)

    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

    $ un+1unΔt+5524NL(un)5924NL(un1)+3724NL(un2)38NL(un3)+N(5524pn5924pn1+3724pn238pn3)=νΔN(7571152un+1+4701152un11181152un5+431152un7),
    $
    (5.6)

    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).

    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

    $ (Nun+1)(Nun)Δt+N(2312NL(un)43NL(un1)+512NL(un2))+ΔN(2312pn43pn1+512pn2)=νNΔN(23un+1+512un1112un3).
    $
    (5.7)

    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:

    $ ΔN(2312pn43pn1+512pn2)=N(2312NL(un)43NL(un1)+512NL(un2)),
    $
    (5.8)

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

    $ (Nun+1)(Nun)Δt=νΔNN(23un+1+512un1112un3).
    $
    (5.9)

    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.

    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:

    $ uΔt,huel(0,T;H2)+νuΔt,huel2(0,T;H3)C(ΔtK+hm),
    $
    (5.10)
    $ pΔt,hpel(0,T;H2)C(ΔtK+hm),
    $
    (5.11)

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

    $ ΔtL3(T,ν),hL4(T,ν).
    $

    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:

    $ Un+1NUnN=tn+1tntUN(,t)dt,2312UnNUnN43Un1NUn1N+512Un2NUn2N=1Δttn+1tnUNUN(,t)dt+τn2,withτ2l2(0,T;L2h)CΔt3UN2H3(0,T;H3)CΔt3,(2312UnNUnN43Un1NUn1N+512Un2NUn2N)=1Δttn+1tn(UNUN)dt+τn3,
    $
    $ withτ3l2(0,T;L2h)CΔt3UNUNH3(0,T;H1)CΔt3,(2312PnN43Pn1N+512Pn2N)=1Δttn+1tnPN(,t)dt+τn4,withτ4l2(0,T;L2h)CΔt3PNH3(0,T;H1)CΔt3,Δ(23Un+1N+512Un1N112Un3N)=1Δttn+1tnΔUN(,t)dt+τn5,withτ5l2(0,T;L2h)CΔt3UNH3(0,T;H2)CΔt3.
    $

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

    $ Un+1UnΔt+2324(UnNUn+(UnUn))23(Un1NUn1+(Un1Un1))+524(Un2NUn2+(Un2Un2))+N(2312Pn43Pn1+512Pn2)νΔN(23Un+1+512Un1112Un3)=τn,
    $
    (5.12)

    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:

    $ τl2(0,T;H1h)C(Δt3+hm).
    $

    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

    $ en+1enΔt+2324(enNUn+unNen+N(en(Un+un)))23(en1NUn1+un1Nen1+N(en1(Un1+un1)))+524(en2NUn2+un2Nen2+N(en2(Un2+un2)))+N(2312qn43qn1+512qn2)=νΔN(23en+1+512en1112en3)+τn.
    $
    (5.13)

    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:

    $ ekNH32+δ1,so thatukNH32+δ˜C0,ukukNL˜C1,
    $
    (5.14)

    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

    $ 1Δt(en+122en22+en+1en22)+νN(43en+1+56en116en3),Nen+1=2312enNUn+43en1NUn1512en2NUn2,en+1+2312unNen+43un1Nen1512un2Nen2,en+1+N(2312enUn+43en1Un1512en2Un2),en+1+N(2312enun+43en1un1512en2un2),en+1+2τn,en+1.
    $
    (5.15)

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

    $ N(43en+1+56en116en3),Nen+1=43Nen+122+56Nen1,Nen+116Nen3,Nen+156Nen+122512Nen122112Nen322.
    $

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

    $ 2312enNUn+43en1NUn1512en2NUn2,en+1116C(en+122+en22+en122+en222),2312unNen+43un1Nen1512un2Nen2,en+1124ν(Nen22+Nen122+Nen222)+C˜C21νen+122,N(2312enUn+43en1Un1512en2Un2),en+1124νNen+122+C(C)2ν(en22+en122+en222),N(2312enun+43en1un1512en2un2),en+1124νNen+122+C˜C21ν(en22+en122+en222).
    $

    Going back to (5.15) results in

    $ en+122en22+34νΔtNen+122124νΔtNen221124νΔtNen122124νΔtNen222112νΔtNen322(C˜C21+C(C)2ν+C)Δt(en+122+en22+en122+en222)+Δtτn22.
    $

    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:

    $ en+122+18νΔtn+1k=1Nek22˜C10(Δt3+hm)2,
    $

    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.

    $ ΔNen+122+nk=0ΔN(ek+1ek)22+124νΔtn+1k=1NΔNek22˜C11(Δt3+hm)2;
    $

    hence,

    $ en+1NH2C(en+1N+Δen+1N)˜C12(Δt3+hm),˜C12=C˜C10+˜C11.
    $

    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

    $ Δt(˜C12)13,h(˜C12)1m,with ˜C12 dependent on ν and exp(T).
    $

    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

    $ un+1unΔt+2312NL(un)43NL(un1)+512NL(un2)+N(2312pn43pn1+512pn2)=νΔN(512un+1+23un112un1),
    $

    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.

    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

    $ ue(x,y,t)=(sin(2πx)cos(2πy),cos(2πx)sin(2πy))Tcos(t),pe(x,y,t)=cos(2πx)cos(2πy)cos(t).
    $
    (6.1)

    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.

    Figure 1.  The discrete $ \ell^2 $ and $ \ell^\infty $ numerical errors vs. temporal resolution $ N_T $ for $ N_T = 100:100:1000 $, with a spatial resolution $ N = 256 $. The numerical results are obtained by the computation using the proposed scheme (3.1), (5.2), (5.5), (5.6), in the first, second, third and fourth order temporal accuracy orders, respectively. The viscosity parameter is taken as $ \nu = 0.5 $. The numerical errors for the velocity variable $ u $ are displayed. The data lie roughly on curves $ CN_T^{-1} $, $ C N_T^{-2} $, $ C N_T^{-3} $ and $ C N_T^{-4} $, respectively, for appropriate choices of $ C $, confirming the full temporal accuracy orders of the proposed schemes.

    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.

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


    Acknowledgments



    The authors acknowledge Anthony Green for proofreading and providing linguistic advice.

    Conflict of interest



    All authors declare no conflicts of interest in this paper

    [1] Cho W, Stahelin RV (2005) Membrane-protein interactions in cell signaling and membrane trafficking. Annu Rev Biophys Biomol Struct 34: 119–151. doi: 10.1146/annurev.biophys.33.110502.133337
    [2] McLaughlin S, Wang J, Gambhir A, et al. (2002) PIP(2) and proteins: interactions, organization, and information flow. Annu Rev Biophys Biomol Struct 31: 151–175. doi: 10.1146/annurev.biophys.31.082901.134259
    [3] Moffett S, Brown DA, Linder ME (2000) Lipid-dependent targeting of G proteins into rafts. J Biol Chem 275: 2191–2198. doi: 10.1074/jbc.275.3.2191
    [4] Lee SY, MacKinnon R (2004) A membrane-access mechanism of ion channel inhibition by voltage sensor toxins from spider venom. Nature 430: 232–235. doi: 10.1038/nature02632
    [5] Gura T (2001) Ancient system gets new respect. Science 291: 2068–2071. doi: 10.1126/science.291.5511.2068
    [6] Simons K, Toomre D (2000) Lipid rafts and signal transduction. Nat Rev Mol Cell Biol 1: 31–39.
    [7] Grouleff J, Irudayam SJ, Skeby KK, et al. (2015) The influence of cholesterol on membrane protein structure, function, and dynamics studied by molecular dynamics simulations. Biochim Biophys Acta 1848: 1783–1795. doi: 10.1016/j.bbamem.2015.03.029
    [8] Fantini J, Barrantes FJ (2013) How cholesterol interacts with membrane proteins: an exploration of cholesterol-binding sites including CRAC, CARC, and tilted domains. Front Physiol 4: 31.
    [9] Segre GV, Goldring SR (1993) Receptors for secretin, calcitonin, parathyroid hormone (PTH)/PTH-related peptide, vasoactive intestinal peptide, glucagonlike peptide 1, growth hormone-releasing hormone, and glucagon belong to a newly discovered G-protein-linked receptor family. Trends Endocrinol Metab 4: 309–314. doi: 10.1016/1043-2760(93)90071-L
    [10] Rymer DL, Good TA (2001) The role of G protein activation in the toxicity of amyloidogenic Abeta-(1–40), Abeta-(25–35), and bovine calcitonin. J Biol Chem 276: 2523–2530. doi: 10.1074/jbc.M005800200
    [11] Kamihira M, Naito A, Tuzi S, et al. (2000) Conformational transitions and fibrillation mechanism of human calcitonin as studied by high-resolution solid-state 13C NMR. Protein Sci 9: 867–877. doi: 10.1110/ps.9.5.867
    [12] Reches M, Porat Y, Gazit E (2002) Amyloid fibril formation by pentapeptide and tetrapeptide fragments of human calcitonin. J Biol Chem 277: 35475–35480. doi: 10.1074/jbc.M206039200
    [13] Wang SS, Good TA, Rymer DL (2005) The influence of phospholipid membranes on bovine calcitonin peptide's secondary structure and induced neurotoxic effects. Int J Biochem Cell Biol 37: 1656–1669. doi: 10.1016/j.biocel.2005.02.006
    [14] Ji S, Wu Y, Sui S (2002) Cholesterol is an important factor affecting the membrane insertion of beta-amyloid peptide (A beta 1–40), which may potentially inhibit the fibril formation. J Biol Chem 277: 6273–6279.
    [15] Micelli S, Meleleo D, Picciarelli V, et al. (2004) Effect of nanomolar concentrations of sodium dodecyl sulfate, a catalytic inductor of alpha-helices, on human calcitonin incorporation and channel formation in planar lipid membranes. Biophys J 87: 1065–1075.
    [16] Micelli S, Meleleo D, Picciarelli V, et al. (2006) Effect of pH-variation on insertion and ion channel formation of human calcitonin into planar lipid bilayers. Front Biosci 11: 2035–2044. doi: 10.2741/1945
    [17] Meleleo D, Micelli S, Toma K, et al. (2006) Effect of eel calcitonin glycosylation on incorporation and channel formation in planar phospholipid membranes. Peptides 27: 805–811. doi: 10.1016/j.peptides.2005.09.011
    [18] Meleleo D, Picciarelli V (2016) Effect of calcium ions on human calcitonin. Possible implications for bone resorption by osteoclasts. Biometals 29: 61–79.
    [19] Diociaiuti M, Polzi LZ, Valvo L, et al. (2006) Calcitonin forms oligomeric pore-like structures in lipid membranes. Biophys J 91: 2275–2281.
    [20] Micelli S, Meleleo D, Picciarelli V, et al. (2004) Effect of sterols on beta-amyloid peptide (AbetaP 1–40) channel formation and their properties in planar lipid membranes. Biophys J 86: 2231–2237.
    [21] Micelli S, Gallucci E, Meleleo D, et al. (2002) Mitochondrial porin incorporation into black lipid membranes: ionic and gating contribution to the total current. Bioelectrochemistry 57: 97–106. doi: 10.1016/S1567-5394(02)00003-8
    [22] Hodge T, Colombini M (1997) Regulation of metabolite flux through voltage-gating of VDAC channels. J Membr Biol 157: 271–279. doi: 10.1007/s002329900235
    [23] Bogdanov M, Dowhan W (1998) Phospholipid-assisted protein folding: phosphatidylethanolamine is required at a late step of the conformational maturation of the polytopic membrane protein lactose permease. EMBO J 17: 5255–5264. doi: 10.1093/emboj/17.18.5255
    [24] Bogdanov M, Dowhan W (1999) Lipid-assisted protein folding. J Biol Chem 274: 36827–36830. doi: 10.1074/jbc.274.52.36827
    [25] Bogdanov M, Umeda M, Dowhan W (1999) Phospholipid-assisted refolding of an integral membrane protein. Minimum structural features for phosphatidylethanolamine to act as a molecular chaperone. J Biol Chem 274: 12339–12345.
    [26] Zhang W, Bogdanov M, Pi J, et al. (2003) Reversible topological organization within a polytopic membrane protein is governed by a change in membrane phospholipid composition. J Biol Chem 278: 50128–50135. doi: 10.1074/jbc.M309840200
    [27] Zhang W, Campbell HA, King SC, et al. (2005) Phospholipids as determinants of membrane protein topology. Phosphatidylethanolamine is required for the proper topological organization of the gamma-aminobutyric acid permease (GabP) of Escherichia coli. J Biol Chem 280: 26032–26038.
    [28] Chaparro Sosa AF, Kienle DF, Falatach RM, et al. (2018) Stabilization of Immobilized Enzymes via the Chaperone-Like Activity of Mixed Lipid Bilayers. ACS Appl Mater Interfaces 10: 19504–19513. doi: 10.1021/acsami.8b05523
    [29] Genheden S, Essex JW, Lee AG (2017) G protein coupled receptor interactions with cholesterol deep in the membrane. Biochim Biophys Acta Biomembr 1859: 268–281. doi: 10.1016/j.bbamem.2016.12.001
    [30] Guixà-González R, Albasanz JL, Rodriguez-Espigares I, et al. (2017) Membrane cholesterol access into a G-protein-coupled receptor. Nat Commun 8: 14505. doi: 10.1038/ncomms14505
    [31] O'Connor JW, Klauda JB (2011) Lipid membranes with a majority of cholesterol: applications to the ocular lens and aquaporin 0. J Phys Chem B 115: 6455–6464. doi: 10.1021/jp108650u
    [32] Bradshaw JP (1997) Phosphatydylglycerol promotes bilayer insertion of salmon calcitonin. Biophys J 72: 2180–2186.
    [33] Stipani V, Gallucci E, Micelli S, et al. (2001) Channel formation by salmon and human calcitonin in black lipid membranes. Biophys J 81: 3332–3338. doi: 10.1016/S0006-3495(01)75966-8
    [34] Westerhoff HV, Hendler RW, Zasloff M, et al. (1989) Interactions between a new class of eukaryotic antimicrobial agents and isolated rat liver mitochondria. Biochim Biophys Acta 975: 361–369. doi: 10.1016/S0005-2728(89)80344-5
    [35] Westerhoff HV, Juretić D, Hendler RW, et al. (1989) Magainins and the disruption of membrane-linked free-energy transduction. Proc Natl Acad Sci USA 86: 6597–6601.
    [36] Matsuzaki K, Mitani Y, Akada K, et al. (1998) Mechanism of synergism between antimicrobial peptides magainin 2 and PGLa. Biochemistry 37: 15144–15153. doi: 10.1021/bi9811617
    [37] Cruciani RA, Barker JL, Zasloff M, et al. (1991) Antibiotic magainins exert cytolytic activity against transformed cell lines through channel formation. Proc Natl Acad Sci USA 88: 3792–3796.
    [38] Gallucci E, Meleleo D, Micelli S, et al. (2003) Magainin 2 channel formation in planar lipid membranes: the role of lipid polar groups and ergosterol. Eur Biophys J 32: 22–32.
    [39] Ashley R, Harroun T, Hauss T, et al. (2006) Autoinsertion of soluble oligomers of Alzheimer's Abeta(1–42) peptide into cholesterol-containing membranes is accompanied by relocation of the sterol towards the bilayer surface. BMC Struct Biol 6: 21. doi: 10.1186/1472-6807-6-21
    [40] Qiu L, Buie C, Reay A, et al. (2011) Molecular dynamics simulations reveal the protective role of cholesterol in β-amyloid protein-induced membrane disruptions in neuronal membrane mimics. J Phys Chem B 115: 9795–9812. doi: 10.1021/jp2012842
    [41] Smart OS, Breed J, Smith GR, et al. (1997) A novel method for structure-based prediction of ion channel conductance properties. Biophys J 72: 1109–1126. doi: 10.1016/S0006-3495(97)78760-5
    [42] Hiller S, Garces RG, Malia TJ, et al. (2008) Solution structure of the integral human membrane protein VDAC-1 in detergent micelles. Science 321: 1206–1210. doi: 10.1126/science.1161302
    [43] Gallucci E, Micelli S, Monticelli G (1996) Pore formation in lipid bilayer membranes made of phosphatidylinositol and oxidized cholesterol followed by means of alternating current. Biophys J 71: 824–831. doi: 10.1016/S0006-3495(96)79283-4
    [44] Micelli S, Gallucci E, Picciarelli V (2000) Studies of mitochondrial porin incorporation parameters and voltage-gated mechanism with different black lipid membranes. Bioelectrochemistry 52: 63–75. doi: 10.1016/S0302-4598(00)00085-4
  • This article has been cited by:

    1. Marco Diociaiuti, Cecilia Bombelli, Laura Zanetti-Polzi, Marcello Belfiore, Raoul Fioravanti, Gianfranco Macchia, Cristiano Giordani, The Interaction between Amyloid Prefibrillar Oligomers of Salmon Calcitonin and a Lipid-Raft Model: Molecular Mechanisms Leading to Membrane Damage, Ca2+-Influx and Neurotoxicity, 2019, 10, 2218-273X, 58, 10.3390/biom10010058
    2. Richard Lantz, Brian Busbee, Ewa P. Wojcikiewicz, Deguo Du, Effects of disulfide bond and cholesterol derivatives on human calcitonin amyloid formation, 2020, 111, 0006-3525, 10.1002/bip.23343
    3. Igor Krivoi, Alexey Petrov, Cholesterol and the Safety Factor for Neuromuscular Transmission, 2019, 20, 1422-0067, 1046, 10.3390/ijms20051046
    4. Lipika Mirdha, Aggregation Behavior of Amyloid Beta Peptide Depends Upon the Membrane Lipid Composition, 2024, 257, 0022-2631, 151, 10.1007/s00232-024-00314-3
  • Reader Comments
  • © 2019 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(5960) PDF downloads(1812) Cited by(4)

Figures and Tables

Figures(4)  /  Tables(2)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog