A VARIABLE TIME-STEP IMEX-BDF2 SAV SCHEME AND ITS SHARP ERROR ESTIMATE FOR THE NAVIER–STOKES EQUATIONS

. We generalize the implicit-explicit (IMEX) second-order backward difference (BDF2) scalar auxiliary variable (SAV) scheme for Navier–Stokes equation with periodic boundary conditions (Huang and Shen, SIAM J. Numer. Anal. 59 (2021) 2926–2954) to a variable time-step IMEX-BDF2 SAV scheme, and carry out a rigorous stability and convergence analysis. The key ingredients of our analysis are a new modified discrete Gr¨onwall inequality, exploration of the discrete orthogonal convolution (DOC) kernels, and the unconditional stability of the proposed scheme. We derive global and local optimal 𝐻 1 error estimates in 2D and 3D, respectively. Our analysis provides a theoretical support for solving Navier–Stokes equations using variable time-step IMEX-BDF2 SAV schemes. We also design an adaptive time-stepping strategy, and provide ample numerical examples to confirm the effectiveness and efficiency of our proposed methods.


Introduction
We here consider numerical approximation of the following incompressible Navier-Stokes (N-S) equations with periodic boundary conditions, where Ω = (−, )  ( = 2, 3) and  > 0 represents the viscosity.How to efficiently and accurately solve the N-S equations has been a research focus for many decades, see [4,[7][8][9]21] and the references therein.While most of the work are concerned with non periodic boundary conditions, the N-S equations with periodic boundary conditions, which retain the basic mathematical properties of N-S equations with non-periodic boundary conditions but can be more efficiently solved with a Fourier-spectral method, are also of important theoretical and practical interests, particularly in the study of well posedness [25] and of homogeneous turbulence [19,20,24].Recently, a high-order IMEX SAV scheme was developed for solving the N-S equations [11], and its numerical solution is shown to be stable without any constraint on time-step size.The results in [11] are established for schemes with constant time-step size.However, in order to efficiently capture the dynamics at different stages of the problem, it is highly beneficial to use variable time-step schemes [6,10,12,13].In fact, the variable time-step BDF2 scheme plays an important role in constructing efficient and accurate algorithms for solving stiff problems, and has been used frequently for solving the parabolic-type equations [3, 14-16, 26, 30].Thus, it is natural to ask if the stability and convergence properties proved in [11] still hold with variable time-step scheme.
The main contributions of this paper are the stability and convergence analysis for the variable time-step IMEX-BDF2 SAV scheme for N-S equations, more precisely, its unconditional stability and optimal convergence order in  1 -norm.Our theoretical results are achieved under the following mild condition on the adjacent timestep ratio A1 : 0 <   ≤  max −  for any small constant 0 <  <  max ≈ 4.8645 and 2 ≤  ≤  , where  max is a root of  3 = 1 + 2.
We point out that the main difficulties of analyzing the variable time-step IMEX-BDF2 SAV scheme are two-fold.On the one hand, since a first-order scheme is used in the first time step, a direct use of the standard discrete Grönwall inequality for the error estimate in the  1 -norm would lead to order reduction.This order reduction of ( 1.5 ) has been observed for linear parabolic equations [27,31] and N-S equations [28].Recently, Ma et al. proved unconditional optimal ( 2 ) convergence in  1 -norm by constructing a modified discrete Grönwall inequality [18].Inspired by the idea in [18], which is used to obtain the optimal convergence order in  1 -norm in [18], we further generalize the discrete Grönwall inequality to make it applicable to the theoretical analysis of IMEX-BDF2 SAV scheme for N-S equations.On the other hand, the SAV and IMEX approaches in the variable time-step IMEX-BDF2 SAV scheme makes the analysis more difficult than a typical semi-implicit numerical scheme.The error estimate of SAV requires the use of the stability of numerical solutions in the  2 -norm, which has not been studied in the framework of the discrete orthogonal convolution (DOC) kernels, and the existing theories require that the time-step ratio satisfies a more strict assumption than   ≤ 4.8645 [28].Therefore, in order to establish a rigorous theories under   ≤ 4.8645, we need to discover new properties of DOC kernels to circumvent the difficulties arisen from applying DOC kernels to IMEX-BDF2 SAV scheme.We point out that for the Newton linearized variable time-step BDF2 scheme, Zhao et al. studied the unconditionally optimal convergence in  2 -norm for general semi-linear equations [31].
The remainder of this paper is organized as follows.In Section 2, we present some preliminary theories that will be used in the paper, including important properties of the N-S equations and of IMEX-BDF2 SAV schemes 2.2.In Section 2.3 we present some important lemmas for the DOC kernels.In Section 3, we present our main results: the global optimal second-order  1 -error estimate of the IMEX-BDF2 SAV scheme in 2D, and a local optimal second-order  1 -error estimate in 3D, and defer their proofs to Section 5.In Section 4, we propose an adaptive time-stepping strategy, and present numerical results to validate our theoretical findings.

Some basic functional settings and the trilinear form
The N-S equations considered here satisfy a periodic boundary condition and will be solved by Fourier spectral methods.We now introduce some relevant functional spaces related to the periodic boundary conditions.If we assume ∫︀ Ω  0 d = 0, then it is easy to see that the solution  of (1.1) also satisfies ∫︀ Ω  d = 0. Hence, in this paper, we assume ∫︀ Ω  0 d = 0 so that ∫︀ Ω  d = 0 at all times.
By using these inequalities, it is straightforward to verify (2.4)

The SAV scheme
Following [11], the form of N-S equations (1.1) can be simplified to a system that only involves .To do so, we take the divergence on both sides of (1.1) and derive By applying the equality (2.2), one has Thus, we equivalently reformulate the N-S equation into Applying the properties of trilinear form (•, •, •), it is easy to check the solution to N-S equations (1.1) satisfies the following energy dissipation law: Based on the energy law, an IMEX-BDFk SAV scheme is developed in [11] by introducing a SAV () = ℰ() + 1, and expand (1.1) as (2.7) We construct below a variable time step IMEX-BDF2 scheme for the above system.We partition the time interval [0,  ] into a general nonuniform time mesh, i.e., 0 =  0 <  1 < • • • <   =  , with a given integer  .Denote by   :=   −  −1 the th time-step size, by  := max 1≤≤   the maximum time-step size, and by   =   / −1 (2 ≤  ≤  ), the adjacent time-step ratio.Denote by Ū  the approximation of the exact solution (  ), and by ∇  Ū  := Ū  − Ū −1 the difference operator.The BDF2 with variable timestep is defined as and two-step Adams-Bashforth extrapolation is defined as We point out that the start of BDF2 scheme needs two steps' information.Differing from the initial setting in [11], we here use BDF1 (i.e., set  1 = 0), and one-step Adams-Bashforth extrapolation to compute first step value  1 .By setting  1 = 0, /(  (1 +   )) and  ()  = 0 for 2 ≤  ≤  − 1, the BDF1 and BDF2 can be reformulated into a unified convolution form of The semi-discrete variable time-step IMEX-BDF2 SAV scheme for (2.6) and (2.7) can be written as: with Next, we construct a Fourier-spectral method for (2.9).For the sake of brevity, we only consider the threedimensional case with Ω = (−, ) 3 , the two-dimensional case can be dealt with similarly.We define the Fourier approximation space as In remainder of this paper, we fix   =   =   =  for simplicity.It is easy to notice that   is a subspace that contains low frequency functions.
It is easy to see that the following properties hold: -Given the initial condition  0 ∈  , then the IMEX-BDF2 SAV scheme (2.9) (resp.(2.10)) admits a unique solution satisfies    (resp.   ∈   ∪  ).-Whenever the pressure is needed in (2.10), it can be computed from Since the mathematical properties of semi-discrete scheme (2.9) and fully discrete scheme (2.10) are similar, so in this paper we only focus on the fully discrete scheme (2.10).
The proof of this theorem is exactly the same to the one in [11], and we omit it here.Note that since Ū   ∈  2 0 (Ω), it also holds ‖   ‖ 1 ≤  for  = 2 thanks to the Poincaré inequality.

Some useful lemmas and their proofs
We now present several Lemmas used in our proof later.The technique of DOC kernels plays a key role in the proof of our main results, and their definition and properties will be presented in this subsection. Set where   := /280, and (, ) :=  • .
Proof.By defining matrix  and Θ in the same way in [30] as where Let us first consider the following equation We claim that   +  is a Hermitian matrix.In fact, according to Lemma 2.3, we see  is a real matrix which has the following estimate: So, by Hermitian matrix's properties, it is easy to check the following estimate holds in complex space: Applying the above inequality to (2.20), we have Here  Θ   2 Θ is a positive definite real matrix, and its inverse matrix is  −1  −2    −1 , which is a positive definite matrix that could be written explicitly as: where (for simplicity, we set  +1 = 0 here)

2𝐷
() This lemma is naturally proved if we can find an uniform up bound for the largest eigenvalue of ̃︀ , i.e., In fact, it is easy to find a proper Λ which is only decided by  max : where Λ : e., Λ ≈ 13.39).By the property of positive definite matrix, it is obvious that  Θ   2 Θ has a smallest eigenvalue, and it is bigger than Λ −1 , which together with (2.21) derive that Taking the definition of Θ,  into the above inequality, the proof is completed.
Lemma 2.6 (Discrete Grönwall inequality).Assume  > 0 and the sequences Lemma 2.6 can be proved by the standard induction hypothesis and is omitted here.To prove the optimal  1 -error estimate, we further introduce a modified Grönwall's inequality, see also [18].

Lemma 2.7 (A modified Grönwall's inequality). Assume a constant 𝐶 and sequences {𝑥
then it holds (2.24) Proof.Denote by Using this property, we can get the estimate of  0 immediately: For the rest of the   , it holds , which together with the inequality above, has (2.27) This lemma is directly proved by noticing where (2.29) One has the following estimate One can refer to Lemma 2.1 and Theorem 3.4 of [30] for the proof of Lemma 2.8 in details.We now deliver a lemma that will only be used in the analysis of local error estimates for 3D case.It is a variant of the lemma in [17], so we omit its proof here.Lemma 2.9.Let  : (0, ∞) → (0, ∞) be continuous and increasing, and let Then there exists a constant  * > 0, which is independent of   > 0 but dependent of  * , satisfies

Main results of variable time-step IMEX-BDF2 SAV scheme
In this section, we state the main results on the global optimal error estimate for the 2D variable time-step IMEX-BDF2 SAV scheme in Theorem 3.1, and the local optimal error estimate for the 3D case in Theorem 3.2.The proofs of these results are based on the stability properties given in Theorem 2.2 and will be deferred to Section 5.
As pointed out in [11], it is no longer possible to obtain a global error estimate in the 3D case.But, the related local error estimate of 3D case in [11] can still be established for variable time-step scheme.Here we again use the DOC kernels to overcome the analysis difficulties raised by variable time-step.
∫︀ ∞ ℱ d/() ( will be given later in the proof ), denote by   :=   −  −1 the th time-step size, by  := max 1≤≤   the maximum time-step size, and ℱ is a positive constant that only depends on real solution .
where the constants  0 ,  Π and  are independent of ,  .
The error analysis for the pressure  is exactly the same as the one in [11].We present the theorem and omit its proof here.Theorem 3.3.Under the same assumptions of Theorem 3.1 in 2D and Theorem 3.2 in 3D, there hold and where  +1  is computed from (2.12),  is a constant independent of  and  ,  * is the biggest integer satisfies ∑︀ * =1   ≤  * , and  * is defined in Theorem 3.2.

Numerical results
We now present three examples to verify the effectiveness and convergence of variable time-step IMEX-BDF2 SAV scheme (2.10).In Example 1, we use the benchmark problems given in [5,11] to investigate the convergence order of our proposed scheme.In Example 2, we design an adaptive time-stepping strategy, and then construct a benchmark problem with different time scales to investigate the effectiveness of our variable time-step scheme, and also present several simulations to investigate the influence of the parameters in the adaptive time-stepping strategy on the CPU-time,  1 -error and the minimum time step-size.In Example 3, we use our adaptive timestepping strategy to deal with the double shear layer problem in [5,11].The experimental results show that, in the same total CPU time, we can successfully solve some problems by our variable time-step IMEX-BDF2 SAV scheme, but, which fails to be solved by the constant time-step IMEX-BDF2 SAV scheme.
Example 1.We first consider the computation of the N-S equation (1.1) in Ω = (−, ) 2 with periodic boundary condition, a benchmark problem in [5], by constructing an exact solution satisfying In the simulations, we take  = 1,  = 16, and divide the time interval [0, 1] by a variable time mesh with  = (1/ ).The maximum  1 errors of   ℎ are listed in Table 1, which shows the second-order convergence of the numerical simulations.According to Theorem 3.1, the optimal error estimate in  1 -norm is at least , which is consistent with the results of our experiments.
In the simulations, we take  = 1,  = 40, and partition the time interval [0, 1] by a variable-time-mesh with  = (1/ ).The maximum  1 errors of   ℎ are listed in Table 2, which again shows the second-order convergence in time.
Overall, both the results in Tables 1 and 2 confirm the analysis results in Theorem 3.1.
Example 2. The main advantage of the variable time-step scheme is its ability to efficiently capture the dynamics of numerical solutions in different time scales.This enables us to obtain smaller errors using less CPU time than the constant time-step scheme.Therefore, how to design an appropriate adaptive time-stepping strategy plays an important role in for the variable time-step scheme.Several common strategies have been used to construct adaptive time stepping schemes.For example, adaptive time-stepping strategies using energy have been well studied in literatures [22,23].However some mathematical indicators are not suitable for designing adaptive time-stepping strategies for N-S equations.In fact, when spatial scalar is fixed and  → 0, it will reduces to energy conservation for inviscid flows, and the energy dissipation is almost negligible.For the 2D fluid problem we consider here (i.e.,  = 2 and  = 0), it will satisfy the Euler equation, the vorticity of the fluid system is also conserved.Therefore, it is not ideal to apply the changes of energy and vorticity directly to the adaptive time-stepping strategy for N-S equations.The adaptive time-stepping strategy in [1] uses a posteriori time error indicator that related to the changes in the adjacent time-step  1 -norm of velocity, and has successfully simulated the unsteady N-S problems.
In this paper, we will design an adaptive time-stepping strategy based on the changes of velocity.A similar strategy is discussed in [1].Our strategy is motivated by an idea that the errors at two adjacent steps should be roughly equal to each other, and an observation that the changes of velocity can effectively reflect the choice of the adaptive time steps.Specifically, our adaptive time strategy is given by where the parameters  and  are the given positive constant to adjust the change of time steps. Noting ℰ( −1

𝑡
) .For the practical simulations, the first step-size  1 is set to be small, but our strategy can guarantee that   will evolve to a proper size in a few steps based on the dynamics of the solution itself.
As a benchmark example, we construct an exact solution having different time scales to verify the superiority of the variable time-step IMEX-BDF2 SAV scheme.In fact, by considering the relative speed of change ‖  ‖/‖‖ =  (), we can choose  () to be a function that is extremely large at some time, such as where  is a large number.We now design the benchmark problem over Ω = (−1, 1) 2 with the following exact solution In the simulations, we take  = 1 and  = 16.The problem will be solved by using the variable and constant time-step strategies, respectively.In the adaptive time-stepping strategy (4.1), we choose  1 = 1.25e−7,  = 0.2 and  max = 5.0e−03.For variable time-step scheme, we choose  = 1.0e−05, 5.0e−06 and 1.0e−06 respectively to solve the problem (4.2) until  = 4.The step sizes at different time levels and the evolution of   in the experiments above are plotted in Figures 1a and 1b respectively.Notes.In Table 3, the first 3 lines are experiment results using variable time-step and the last 3 lines are results using constant time-step.We re-emphasize that  := max 1≤≤   here.As a contrast, we solve the same problem using constant-step sizes  = 1/1250, 1/2500 and 1/5000, respectively.The numerical results are listed in Table 3. Table 3 shows that the variable time-step scheme gives a lower  1 -error with lower computational costs, while the constant time-step scheme gives a large  1 -error at the almost same computational cost.For example, the  1 -error of variable time-step scheme with  = 1.0e−06 and the constant time-step scheme with  = 2.0e−04 is similar, but the CPU-time of constant time-step scheme is almost five times than the variable time-step one.This can also be noticed by observing the evolution of   over time in Figure 1b, and it is noticed that   is much closer to 1 for variable time-step scheme, while the   is much farther than 1 for constant time-step scheme.Figure 1a shows that the variable time-step IMEX-BDF2 SAV scheme captures changes over different time scales very well, and by setting a lower , one can get a more accurate numerical result.
We further investigate how the parameters in the adaptive time-stepping strategy influence the CPU-time,  1 -error and the minimum time step-size.Table 4 shows the CPU-time,  1 -error and the minimum time stepsize (since the initial step-size is very small, we only consider steps after the 100th step here) for different .In this experiment, all numerical solutions are obtained by fixing  = 16,  = 0.2 and  max = 5.0e−03.One can observe in Table 4 that as  becomes smaller, the numerical error and the minimum time-step size also become smaller, which implies one can get more accurate numerical results by lowering .However, smaller  may result in longer CPU-time, and how to select appropriate  in actual computation is worthy of further study.
Table 5 shows the CPU-time and  1 -error for different parameter  max in adaptive time-stepping strategy.Numerical solutions were obtained by setting  = 16,  = 0.2 and  = 4.0e−05.One can clearly observe that when the  max gets smaller, one can have a more accurate result.By setting an appropriate  max , one can get a much smaller error at a negligible cost in CPU-time than with no  max .But, as the  max becomes smaller, the variable time-step scheme under this strategy will gradually degenerate into a constant time-step scheme, and we will spend much more CPU-time to obtain a smaller  1 -error.
From the above results, one can see that our adaptive time-stepping strategy (4.1) is suitable for the variable time-step IMEX-BDF2 SAV scheme, which significantly improves the computational efficiency for this benchmark problem.
Example 3. We also use our adaptive time-stepping strategy to carry out experiments on the well known e double shear layer problem [5,11].Consider the N-S equation (1.1) in Ω = (−0.5,0.5) 2 with periodic boundary conditions and the initial conditions given by where  determines the slope of the shear layer and  represents the size of the perturbation.As the same setting as [11], we fix  = 0.05,  = 100 and  = 0.00005.
In [11] it has been shown that the constant time-step IMEX-BDF2 SAV scheme can compute a reasonably correct solution at  = 1.2 by taking  = 2.5e−04 and  = 128, and when the step-size is  = 3.0e−04, the solution at  = 1.2 will blow up.However, through a large number of experiments and observations, we find that it may not be sufficient to use only the vorticity contours as in [11], or to use the magnitude of  deviation from 1 to determine whether the numerical results are correctly computed.In fact, if we examine in the computation of the solution, when using step-size  = 2.5e−04, we can see that the velocity solved by the constant time-step scheme is not accurate enough at  = 1.2.However, by using the variable time-step IMEX-BDF2 SAV scheme, we may compute correct and accurate results with less computational cost.
In the simulations, we take  = 128.The problem is solved by using the variable and constant time-step schemes, respectively.With the adaptive time-stepping strategy (4.1), we choose  1 = 1.0e−04 and  max = ℰ(  ) . Figure 2a shows that the variable time-step scheme can compute the correct relative change of velocity on [0, 1.2] with less computational cost, while the relative change of velocity computed by the constant time-step scheme with the same or even higher computational cost will blow up.For example, the numerical result of constant time-step scheme using  = 2.8e−04 blows up at around 0.8.Meanwhile, for the variable time-step scheme, the average step-size used in the simulation from 0 to 1.2 is 2.803e−04, but it can compute the correct relative change on [0, 1.2].We further consider the experiment in [11], i.e., using stepsize  = 2.5e−04 to solve this problem.Although its   does not deviate much from 1 [11], and the vorticity contour is correctly simulated, we can still see the relative change of velocity blows up after 1.05.However, the variable time-step scheme still works.The result of variable time-step scheme with less computational cost behaves better, and the relative change of velocity maintains correct over the time interval [0, 1.2].As shown in Figure 2b, it can be seen that the variable time-step IMEX-BDF2 SAV scheme well captures the changes of different time scales in the double shear layer problem.This is also the reason why it can obtain correct numerical results at a lower computational cost than using the constant time-step scheme.
As a contrast, we solve the same problem using constant-step sizes   = 1.0e−04, 2.5e−04 and 2.8e−04, respectively, and the relative changes of velocity at different time level, i.e.,

ℰ(𝑈
ℰ(  ) , are also plotted in Figure 2a.Table 6 shows the results of constant time-step IMEX-BDF2 SAV scheme by using  = 1.0e−04 as the reference solution, and  1 -error of velocity.It can be seen that the error of variable step-size scheme is smaller.
In order to make it easier for the readers to realize the difference in results, we compute the vorticity contours of the difference between numerical results according to the following equation Here    are results using constant step-size 2.5e−04 and 2.8e−04 and variable step-size when setting  = 2.4e−05, and   ref is the reference solution, i.e., result using constant step-size 1.0e−04.The Figure 3 shows the vorticity contours of the difference between numerical solution    and reference solution   ref at  = 1.2.

Proofs of Theorems 3.1 and 3.2
We now present the proofs of Theorems 3.1 and 3.2.We carry out a complete rigorous analysis of the error estimate in the 2D case, and only present the part of the proof for the 3D case which is different from the 2D case.For brevity, the proof for the following error estimate is divided into three steps, and the technique of DOC kernels is mainly used for the theoretical analysis in the second step.

Proof of Theorem 3.1
Proof.Similarly as in [11], the main task is to establish by induction , and we want to show |1 −  +1 | ≤  0  +  Π  2− where  0 ,  Π will be determined later.
We point out that  4 is bounded thanks to (5.4).By applying the modified Grönwall's inequality in Lemma 2.7 to (5.14), we can prove: (5.15) The summation terms on the righthand side can be bounded as follows.
Hence, there exists a constan  5 > 0 such that )︂ ≤  5 . (5.17) On the other hand, thanks to Lemma 2.8, we have Similarly, Using (5.17)-(5.19),we derive from (5.15) that where  8 ,  9 are constants independent of  0 ,  Π .By simple computation, one also finds where  10 is a constant depending on the exact solution () only.The proof here can be easily generalized to prove the boundedness of The bound here is independent of  0 ,  Π and only decided by the parameters , ,  max , Ω and the exact solution .For simplicity, we denote all the up bound here as  11 .

Proof of Theorem 3.2
Proof.A main difficulty in the three dimensional case is that Theorem 2.2 only provides  2 -stability for    , while in the two dimensional case,  1 -stability for    is provided in Theorem 2.2.Thus, our main task here is to prove the stability of    in  1 -norm.The rest of the proof is similar to the one for Theorem 3.1.We proceed by induction.Assuming |1 −   | ≤  0  +  Π  2− , ∀ ≤ , we shall prove |1 −  +1 | ≤  0  +  Π  2− where  0 ,  Π will be determined below, but we may need them greater than 1 without loss of generality.
Clearly, it is satisfied for  = 0. Step Multiplying the first equation in (2.10) by DOC kernels, and summing up from  = 0 to  − 1, we get )︁ .
Thanks to the Poincaré's inequality, there holds the following inequality: )︂ +   .
Here   is a constant decided by  0 .By applying lemma 2.9, and choosing () =  2 +  3 , there exist 0 <  * < ∫︀ ∞ ℱ d () and  * > 0 such that With the above bound, we can then prove the desired result by following similar procedures in Steps 2 and 3 in the proof of Theorem 3.1.

Concluding remarks
We considered in this paper an energy stable variable time-step IMEX-BDF2 SAV scheme for the Navier-Stokes equations.Based on the energy stability, we proved the sharp global optimal error estimates by using DOC kernels through rigorous mathematical induction process.We point out that our results are obtained without any restriction on the ratio of time step and spatial grid sizes.
We introduced a suitable adaptive time-stepping strategy, and provided numerical examples to verify that the variable time-step size IMEX-BDF2 SAV scheme can effectively obtain correct numerical solutions with lower computational cost compared with a corresponding constant time-step scheme.
We only considered the analysis for the Navier-Stokes equations with periodic boundary conditions, although a similar variable time-step IMEX-BDF2 SAV scheme can be developed for the Navier-Stokes equation with non-periodic boundary conditions based on the constant step size scheme in [29].However, the analysis for the non-periodic case is much more difficult due to the lack of a priori bound on the pressure approximation.But the many techniques that we developed in this paper will be useful in the future study of the variable time-step IMEX-BDF2 SAV scheme for the Navier-Stokes equation with non-periodic boundary conditions.

Table 3 .
Errors and parameters of experiments for solving (4.2).

Table 5 .
Numerical results by choosing different  max for solving (4.2).

Table 6 .
Errors and parameters of experiments for solving double shear layer problem.