Convergence Analysis of a Fully Discrete Family of Iterated Deconvolution Methods for Turbulence Modeling with Time Relaxation

Abstract

We present a general theory for regularization models of the Navier-Stokes equations based on the Leray deconvolution model with a general deconvolution operator designed to fit a few important key properties. We provide examples of this type of operator, such as the (modified) Tikhonov-Lavrentiev and (modified) Iterated Tikhonov-Lavrentiev operators, and study their mathematical properties. An existence theory is derived for the family of models and a rigorous convergence theory is derived for the resulting algorithms. Our theoretical results are supported by numerical testing with the Taylor-Green vortex problem, presented for the special operator cases mentioned above.

1. Approximate Deconvolution for Turbulence Modeling

Numerical simulations of complex flows present many challenges. The Navier-Stokes (NS) equations (NSE), given by the following𝐮_𝑡+𝐮⋅∇𝐮−𝜈Δ𝐮+∇𝑝=𝐟,∇⋅𝐮=0,inΩ×(0,𝑇) (1.1) is an exact model for the flow of a viscous, incompressible fluid, [1]. For turbulent flows (characterized by Reynold’s number Re≫1), it is infeasible to properly resolve all significant scales above the Kolmogorov length scale 𝑂(Re^−3/4) by direct numerical simulation. Thus, numerical simulations are often based on various regularizations of NSE, rather than NSE themselves. Accordingly, regularization methods provide a (computationally) efficient and (algorithmically) simple family of turbulence models. Several of the most commonly applied regularization methods include: Leray𝐰_𝑡+𝐰⋅∇𝐰−Re⁻¹()𝐰Δ𝐰+∇𝑝=𝐟,∇⋅𝐰=0,NS-𝛼_𝑡+𝐰×(∇×𝐰)−Re⁻¹(Δ𝐰+∇𝑃=𝐟,NS-𝜔)𝐰_𝑡+𝐰×(∇×𝐰)−Re⁻¹()𝐰Δ𝐰+∇𝑃=𝐟,timerelaxation_𝑡+𝐰⋅∇𝐰−Re⁻¹Δ𝐰+∇𝑝+𝜒(𝐰−𝐰)=𝐟, (1.2)(1.3)(1.4)(1.5) where ∇⋅𝐰=0 in each case and 𝐰 is an averaged velocity field 𝐰, 𝑝 is pressure, and 𝑃 is the Bernoulli pressure. More details about these models can be found for instance in [2–5] and references therein. Although these regularization methods achieve high theoretical accuracy and perform well in select practical tests, those models do not provide a fully-developed numerical solution for decoupling the scales in a turbulent flow. In fact, results show that only time relaxation regularization truncates scales sufficiently for practical computations. Indeed, it is shown that time relaxation term 𝜒(𝐰−𝐰) for 𝜒>0 damps unresolved fluctuations over time [5, 6]. Note that the choice of 𝜒 is an active area of research and that solutions are very sensitive to variations in 𝜒.

Deconvolution-based regularization is also an active area of research obtained, for example, by replacing 𝐰 by 𝐷(𝐰) in each (1.2)–(1.5) for some deconvolution operator 𝐷. In [7], Dunca proposed the general Leray-deconvolution problem (𝐷(𝐰) instead of 𝐰) as a more accurate extension to Leray’s model [8]. Leray used the Gaussian filter as the smoothing (averaging) filter 𝐺, denoted above by overbar. In [9], Germano proposed the differential filter (approximate-Gaussian) 𝐺=(−𝛿²Δ+𝐼)⁻¹ where 𝛿>0 is the filter length. The differential filter is easily modeled in the variational framework of the finite element (FE method (FEM)). We provide a brief overview of continuous and discrete operators (Section 3).

The deconvolution-based models have proven themselves to be very promising. However, among the very many known approximate deconvolution operators from image processing, for example, [10], so far only few have been studied for turbulence modeling, for example, the van Cittert and the modified Tikhonov-Lavrentiev deconvolution operators. Their success suggests that it is time to develop a general theory for regularization models of the NSE as a guide to development of models based on other, possibly better, deconvolution operators and refinement of existing ones.

Herein, we present a general theory for regularization models of the NSE based on the Leray deconvolution model with a general deconvolution operator. We prove energetic stability (and hence existence) and convergence of an FE (in space) and Crank-Nicolson (CN) (in time) discretization of the following family of Leray deconvolution regularization models with time relaxation: find 𝐰∶Ω×(0,𝑇]→ℝ³ and 𝜋∶Ω×(0,𝑇]→ℝ satisfying the following 𝐰_𝑡+𝐷(𝐰)⋅∇𝐰−Re⁻¹Δ𝐰+∇𝜋+𝜒(𝐰−𝐷(𝐰))=𝐟,∇⋅𝐰=0, (1.6) for some appropriate boundary and initial conditions (the fully discrete model is presented in Problem 1 with energetic stability and well posedness proved in Theorem 4.5).

1.1. Improving Accuracy of Approximate Deconvolution Methods

The fundamental difficulties corresponding to regularization methods applied as a viable turbulence model include ensuring that:(i)scales are appropriately truncated (model microscale = filter radius = mesh width)(ii)smooth parts of the solution are accurately approximated (𝐷(𝐰)≈𝐰 for smooth 𝐰)(iii)physical fidelity of flow is preserved.

Due to the nonlinearity in (1.6), different choices of the filter and deconvolution operator yield significant changes in the solution of the corresponding model. Implementation concerns for deconvolution methods, for example, Tikhonov-Lavrentiev regularization given by 𝐷=(𝐺+𝛼𝐼)⁻¹, include selection of deconvolution parameter 𝛼>0. Iterated deconvolution methods reduce approximation sensitivity relative to 𝛼-selection and, hence, allow a conservatively large 𝛼-selection for stability with updates (fixed number of iterations) used to recover higher accuracy (Section 3.3). For example we prove, under usual conditions, (Proposition 5.3), ModifiedTikhonov𝑗=0error(𝐰−𝐷(𝐰))≤𝒪𝛼𝛿²,IteratedModifiedTikhonov𝑗>0error(𝐰−𝐷(𝐰))≤𝒪𝛼𝛿^2𝑗+1, (1.7) so that iterated modified Tikhonov regularization gives geometric convergence with respect to the update number 𝑗. In either case, Tikhonov-Lavrentiev or iterated Tikhonov-Lavrentiev regularization, we prove, under usual conditions, (Proposition 5.4) error𝐷^ℎ𝐰^ℎ−𝐷(𝐰)≤𝒪ℎ^𝑘, (1.8) where 𝐷^ℎ and 𝐰^ℎ represent the discrete deconvolution operator and discrete filter, respectively, and 𝑘 is order of FE polynomial space. We propose minimal properties for a general family of deconvolution operators 𝐷 and filters 𝐺 (Section 3.1) satisfying (Assumptions 3.4, 3.5), for example,(i)||𝐷𝐺||≤1, forces spectrum of 𝐷𝐺 in [0,1](ii)||∇𝐷𝐺𝐯||≤𝑑₁||∇𝐯||, controls size of ∇𝐷𝐺(iii)||(𝐼−𝐷𝐺)𝐯||≤𝑐₁(𝛼,𝛿)||Δ𝐯||→0 as 𝛼, 𝛿→0 (for smooth 𝐯), ensures convergence of method(iv)||(𝐷𝐺−𝐷^ℎ𝐺^ℎ)𝐯||≤𝑐₂(ℎ,𝛼,𝛿)→0 as ℎ→0 (for smooth 𝐯), ensures convergence of method.

In fact, the updates {𝝎_𝑗}_𝑗≥0 satisfying 𝝎₀∶=𝐷(𝐰),𝝎_𝑗−𝝎_𝑗−1∶=𝐷𝐰−𝝎_𝑗−1, (1.9) inherit the properties assumed for the base operator 𝐷 in Assumptions 3.4, 3.5 (Propositions 3.10, 3.11). These iterates represent defect correction generalization of iterated Tikhonov regularization operator [11]. We prove that the FE-CN approximation 𝐰^ℎ of the general deconvolution turbulence model (Problem 1) satisfies (Theorem 4.6) error𝐰^ℎ−𝐮_NSE≤𝐶ℎ^𝑘+Δ𝑡²+𝑐₁𝛼,𝛿,𝑗+𝑐₂ℎ,𝛼,𝛿,𝑗→0,asℎ,𝛿,𝛼→0 (1.10) for smooth enough solutions 𝐮_NSE of the NSE (see variational formulation (2.3)–(2.5)). We show that 𝑐₁=𝐶(𝛼𝛿²)^𝛽 (for any 0≤𝛽≤𝑗+1 and Δ^𝛽𝐮_NSE∈𝐿²(Ω)) (Proposition 5.3) and 𝑐₂=𝐶ℎ^𝑘 (for 𝐷(𝐰)∈𝐻^𝑘+1(Ω)) (Proposition 5.4) for modified iterated Tikhonov-Lavrentiev regularization (see Corollary 5.5 for corresponding error estimate). We conclude with a numerical test that verifies the theoretical convergence rate predicted in Theorem 4.6 (Section 5.2).

1.2. Background and Overview

One of the most interesting approaches to generate turbulence models is via approximate deconvolution or approximate/asymptotic inverse of the filtering operator. Examples of such models include: Approximate Deconvolution Models (ADM) and Leray-Tikhonov Deconvolution Models. Layton and Rebholz compiled a comprehensive overview and detailed analysis of ADM [12] (see also references therein). Previous analysis of the ADM with and without the time-relaxation term used van Cittert deconvolution operators [5, 6]; although easily programmed, van Cittert schemes can be computationally expensive [5]. Tikhonov-Lavrentiev regularization is another popular regularization scheme [13]. Determining the appropriate value of 𝛼 to ensure stability while preserving accuracy is challenging, see for example, [14–19]. Alternatively, iterated Tikhonov regularization is well known to decouple stability and accuracy from the selection of regularization parameter 𝛼, see for example, [11, 20–22]. Iterated Tikhonov regularization is one special case of the general deconvolution operator we propose herein.

2. Function Spaces and Approximations

Let the flow domain Ω∈ℝ^𝑑 for 𝑑=2,3 be a regular and bounded polyhedral. We use standard notation for Lebesgue and Sobolev spaces and their norms. Let ||⋅|| and (⋅,⋅) be the 𝐿²-norm and inner product, respectively. Let ||⋅||_𝑝,𝑘∶=||⋅||_𝑊^𝑘_𝑝(Ω) represent the 𝑊^𝑘_𝑝(Ω)-norm. We write 𝐻^𝑘(Ω)∶=𝑊^𝑘₂(Ω) and ||⋅||_𝑘 for the corresponding norm. Let the context determine whether 𝑊^𝑘_𝑝(Ω) denotes a scalar, vector, or tensor function space. For example let 𝐯∶Ω→ℝ^𝑑. Then, 𝐯∈𝐻¹(Ω) implies that 𝐯∈𝐻¹(Ω)^𝑑 and ∇𝐯∈𝐻¹(Ω) implies that ∇𝐯∈𝐻¹(Ω)^𝑑×𝑑. Write 𝑊^𝑚_𝑞(𝑊^𝑘_𝑝(Ω))∶=𝑊^𝑚_𝑞(0,𝑇;𝑊^𝑘_𝑝(Ω)) equipped with the standard norm. For example, ‖𝐯‖_𝐿^𝑞_(𝑊^𝑘_𝑝)∫∶=^𝑇0‖𝐯(⋅,𝑡)‖^𝑞_𝑝,𝑘𝑑𝑡^1/𝑞,if1≤𝑞<∞,esssup_0<𝑡<𝑇‖𝐯(⋅,𝑡)‖_𝑝,𝑘,if𝑞=∞. (2.1)

Denote the pressure and velocity spaces by 𝑄∶=𝐿²₀(Ω)={𝑞∈𝐿²(Ω)∶∫_Ω𝑞=0} and 𝑋∶=𝐻¹₀(Ω)={𝐯∈𝐻¹(Ω)∶𝐯|_𝜕Ω=0}, respectively. Moreover, the dual space of 𝑋 is denoted 𝑋^′∶=𝑊₂⁻¹(Ω) and equipped with the norm ‖𝐟‖₋₁∶=sup_{0≠𝐯∈𝑋}⟨𝐟,𝐯⟩_{𝑋′×𝑋}|𝐯|₁. (2.2)

Fix 𝐟∈𝑋^′ and Re>0. In this setting, we consider strong NS solutions: find 𝐮∈𝐿²(𝑋)∩𝐿^∞(𝐿²(Ω)) and 𝑝∈𝑊^−1,∞(𝑄) satisfying 𝐮_𝑡,𝐯+(𝐮⋅∇𝐮,𝐯)+Re⁻¹()−=()(](]𝐮()∇𝐮,∇𝐯𝑝,∇⋅𝐯𝐟,𝐯,a.e.𝑡∈0,𝑇,∀𝐯∈𝑋,𝑞,∇⋅𝐮=0,a.e.𝑡∈0,𝑇,∀𝑞∈𝑄,⋅,0=𝐮₀in𝐿²(Ω). (2.3)(2.4)(2.5) Let 𝑉∶={𝐯∈𝑋∶∇⋅𝐯=0}. Restricting test functions 𝐯∈𝑉 reduces (2.3)–(2.5) to find 𝐮∶(0,𝑇]→𝑉 satisfying 𝐮_𝑡,𝐯+(𝐮⋅∇𝐮,𝐯)+Re⁻¹(∇𝐮,∇𝐯)=(𝐟,𝐯),a.e.𝑡∈(0,𝑇],∀𝐯∈𝑉 (2.6) and (2.5). For smooth enough solutions, solving the problem associated with (2.6), (2.5) is equivalent to (2.3)–(2.5).

Control of the nonlinear term is essential for establishing a priori estimates and convergence estimates. We state a selection of inequalities here that will be utilized later: ||(Ω)‖𝐮‖𝐮⋅∇𝐯,𝐰≤𝐶¹|𝐯|₁‖𝐰‖₁‖𝐮‖‖𝐯‖₂‖𝐰‖₁∀𝐯∈𝐻²(Ω)‖𝐮‖₂‖𝐯‖‖𝐰‖₁∀𝐮∈𝐻²(Ω)∩𝑉. (2.7)

2.1. Discrete Function Setting

Fix ℎ>0. Let 𝒯^ℎ be a family of subdivisions (e.g., triangulation) of Ω⊂ℝ^𝑑 satisfying Ω=⋃_{𝐸∈𝒯ℎ}𝐸 so that 𝑑𝑖𝑎𝑚𝑒𝑡𝑒𝑟(𝐸)≤ℎ and any two closed elements 𝐸₁, 𝐸₂∈𝒯^ℎ are either disjoint or share exactly one face, side, or vertex. See Chapter II, Appendix A in [23] for more on this subject in context of Stokes problem and [24] for a more general treatment. For example, 𝒯^ℎ consists of triangles for 𝑑=2 or tetrahedra for 𝑑=3 that are nondegenerate as ℎ→0.

Let 𝑋^ℎ⊂𝑋 and 𝑄^ℎ⊂𝑄 be a conforming velocity-pressure mixed FE space. For example, let 𝑋^ℎ and 𝑄^ℎ be continuous, piecewise (on each 𝐸∈𝒯^ℎ) polynomial spaces. The discretely divergence-free space is given by 𝑉^ℎ=𝐯^ℎ∈𝑋^ℎ∶𝑞^ℎ,∇⋅𝐯^ℎ=0∀𝑞^ℎ∈𝑄^ℎ. (2.8) Note that in general 𝑉^ℎ⊄𝑉 (e.g., for Taylor-Hood elements). In order to avoid stability issues arising when FE solutions are not exactly divergence free (i.e., when 𝑉^ℎ⊄𝑉), we introduce the explicitly skew-symmetric convective term 𝑏^ℎ()1𝐮,𝐯,𝐰∶=2((𝐮⋅∇𝐯,𝐰)−(𝐮⋅∇𝐰,𝐯)), (2.9) so that 𝑏^ℎ(𝐮,𝐯,𝐯)=0. (2.10) Note that 𝑏^ℎ(𝐮,𝐯,𝐰)=(𝐮⋅∇𝐯,𝐰) when 𝐮∈𝑉. Moreover, the trilinear from 𝑏^ℎ(⋅,⋅,⋅) is continuous and skew-symmetric on 𝑋×𝑋×𝑋.

Lemma 2.1. If 𝐮,𝐯,𝐰∈𝑋, 𝑏^ℎ((Ω‖𝐮𝐮‖𝐮,𝐯,𝐰)≤𝐶)‖‖₁^1/2‖𝐯‖₁‖𝐰‖₁‖𝐮‖‖𝐯‖₂‖𝐰‖₁∀𝐯∈𝐻²(Ω). (2.11)

Proof. The proof of the first inequality can be found in [25]. The second follows from Hölder’s and Poincaré’s inequalities.

For the time discretization, let 0=𝑡₀<𝑡₁<⋯<𝑡_𝑀−1=𝑇<∞ be a discretization of the time interval [0,𝑇] for a constant time step Δ𝑡=𝑡_𝑛+1−𝑡_𝑛. Write 𝑡_𝑛+1/2∶=(𝑡_𝑛+1+𝑡_𝑛)/2,𝑧_𝑛=𝑧(𝑡_𝑛) and, if 𝑧∈𝐶⁰([𝑡_𝑛,𝑡_𝑛+1]), 𝑧_𝑛+1/2=(𝑧_𝑛+1+𝑧_𝑛)/2. Define ‖𝐮‖_𝑙^𝑞_{(𝑚1,𝑚2;𝑊}^𝑘_𝑝(Ω))∶=Δ𝑡∑^𝑚2_𝑛=𝑚1𝐮_𝑛^𝑞_𝑘,𝑝^1/𝑞,𝑞∈[1,∞),max_{𝑚1≤𝑛≤𝑚2}𝐮_{𝑛𝑘,𝑝},𝑞=∞, (2.12) for any 0≤𝑛=𝑚₁,𝑚₁+1,…,𝑚₂≤𝑀. Write ||𝐮||_𝑙^𝑞_(𝑊^𝑘_𝑝(Ω))=||𝐮||_𝑙^𝑞_(0,𝑀;𝑊^𝑘_𝑝(Ω)). We say that 𝐮∈𝑙^𝑞(𝑚₁,𝑚₂;𝑊^𝑘_𝑝(Ω)) if the associated norm defined above stays finite as Δ𝑡→0.

The discrete Gronwall inequality is essential to the convergence analysis in Section 4.2.

Lemma 2.2. Let 𝐷≥0 and 𝜅_𝑛, 𝐴_𝑛, 𝐵_𝑛, 𝐶_𝑛≥0 for any integer 𝑛≥0 and satisfy 𝐴_𝑀+Δ𝑡^𝑀_𝑛=0𝐵_𝑛≤Δ𝑡^𝑀_𝑛=0𝜅_𝑛𝐴_𝑛+Δ𝑡^𝑀_𝑛=0𝐶_𝑛+𝐷,∀𝑀≥0. (2.13) Suppose that for all 𝑛, Δ𝑡𝜅_𝑛<1 and set 𝑔_𝑛=(1−Δ𝑡𝜅_𝑛)⁻¹. Then, 𝐴_𝑀+Δ𝑡^𝑀_𝑛=0𝐵_𝑛≤expΔ𝑡^𝑀_𝑛=0𝑔_𝑛𝜅_𝑛Δ𝑡^𝑀_𝑛=0𝐶_𝑛+𝐷,∀𝑀≥0. (2.14)

Proof. The proof follows from [26].

2.2. Approximation Theory

Let 𝐶>0 be a generic constant independent of ℎ→0⁺. Preserving an abstract framework for the FE spaces, we assume that 𝑋^ℎ×𝑄^ℎ inherit several fundamental approximation properties.

Assumption 2.3. The FE spaces 𝑋^ℎ×𝑄^ℎ satisfy:

Uniform inf-sup (LBB) condition inf_{𝑞ℎ∈𝑄ℎ}sup_{𝐯ℎ∈𝑋ℎ}𝑞^ℎ,∇⋅𝐯^ℎ𝐯^ℎ₁𝑞≥𝐶>0. (2.15) FE-approximation inf_{𝐯ℎ∈𝑋ℎ}𝐮−𝐯^ℎ₁≤𝐶ℎ^𝑘‖𝐮‖_𝑘+1,inf_{𝑞ℎ∈𝑄ℎ}𝑝−𝑞^ℎ≤𝐶ℎ^𝑠+1𝑝_𝑠+1for𝑘≥0,𝑠≥−1when𝐮∈𝐻^𝑘+1(Ω)∩𝑋,𝑝∈𝐻^𝑠+1(Ω)∩𝑄. (2.16) Inverse-estimate 𝐯^ℎ₁≤𝐶ℎ⁻¹𝐯^ℎ,∀𝐯^ℎ∈𝑋^ℎ. (2.17) The well-known Taylor-Hood mixed FE is one such example satisfying Assumption 2.3.

Estimates in (2.18)–(2.20) stated below are used in proving error estimates for time-dependent problems: for any 𝑛=0,1,…,𝑀−1, 𝜽_𝑛+1−𝜽_𝑛Δ𝑡²_𝑘≤𝐶Δ𝑡^{−1𝑡𝑛+1}_𝑡𝑛𝜽_𝑡(𝑡)²_𝑘𝜽𝑑𝑡,_𝑛+1/2−𝜽𝑡_𝑛+1/2²_𝑘≤𝐶Δ𝑡^3𝑡𝑛+1_𝑡𝑛𝜽_𝑡𝑡(𝑡)²_𝑘1𝑑𝑡,𝜽Δ𝑡^𝑛+1−𝜽_𝑛−𝜽_𝑡n+1/2²_𝑘≤𝐶Δ𝑡^3𝑡𝑛+1_𝑡𝑛𝜽_𝑡𝑡𝑡(𝑡)²_𝑘𝑑𝑡, (2.18)(2.19)(2.20) where 𝜽∈𝐻¹(𝐻^𝑘(Ω)), 𝜽∈𝐻²(𝐻^𝑘(Ω)), and 𝜽∈𝐻³(𝐻^𝑘(Ω)) is required, respectively, for some 𝑘≥−1. Each estimate (2.18)–(2.20) is a result of a Taylor expansion with integral remainder.

These higher-order spatial (𝑘≥2 or 𝑠≥1) and temporal estimates (2.18)–(2.20) require that the nonlocal compatibility condition addressed by Heywood and Rannacher in [26, 27] (and more recently, for example, by He in [28, 29] and He and Li in [30, 31]) is satisfied. Suppose, for example, that 𝑝₀ is the solution of the (well-posed) Neumann problem Δ𝑝₀=∇⋅𝐟₀−𝐮₀⋅∇𝐮₀,inΩ,∇𝑝₀⋅𝐧_𝜕Ω=Δ𝐮₀+𝐟₀−𝐮₀⋅∇𝐮₀⋅𝐧_𝜕Ω. (2.21) In order to avoid the accompanying factor min{𝑡⁻¹,1} in the error estimates contained herein, the following compatibility condition is necessarily required (e.g., see [27, Corollary 2.1]): ∇𝑝_0𝜕Ω=Δ𝐮₀+𝐟₀−𝐮₀⋅∇𝐮_0𝜕Ω. (2.22) Replacing (2.21) with (2.21)(a), (2.22) defines an overdetermined Neumann-type problem. Condition (2.22) is a nonlocal condition relating 𝐮₀ and 𝐟₀. Condition (2.22) is satisfied for several practical applications including start up from rest with zero force, 𝐮₀=0, 𝐟₀=0. In general, however, condition (2.22) cannot be verified. In this case, it is shown that, for example, |𝐮_𝐭(⋅,𝑡)|₁,||𝐮(⋅,𝑡)||₃→∞as𝑡→0⁺.

We finish with an approximation property of the 𝐿²-projection. Indeed, Assumption 2.4 holds for smooth enough Ω.

Assumption 2.4. Fix 𝐰∈𝑉 and let 𝐰^ℎ∈𝑉^ℎ be the unique solution satisfying (𝐰−𝐰^ℎ,𝐯^ℎ) for all 𝐯^ℎ∈𝑉^ℎ. Then 𝐰−𝐰^ℎ_−𝑚≤𝐶ℎ^𝑚+1inf_{𝐯ℎ∈𝑋ℎ}𝐰−𝐯^ℎ₁ (2.23) for 𝑚=−1,0,1.

Note that the infimum in (2.23) is over all 𝑋^ℎ (see intermediate estimate (1.16) of Theorem II.1.1 in [23] for the corresponding estimate relating the spaces 𝑉^ℎ and 𝑋^ℎ).

3. Filters and Deconvolution

We prescribe the essential properties our filter 𝐺 and deconvolution operator 𝐷 in this section.

Definition 3.1. Let 𝑌 be a Hilbert space and 𝑇∶𝑌→𝑌. Write 𝑇≥0 if 𝑇 is self-adjoint 𝑇=𝑇^′ and (𝑇𝐯,𝐯)_𝑌≥0 for all 𝐯∈𝑌 and call 𝑇 symmetric nonnegative (snn). Write 𝑇>0 if 𝑇 is self-adjoint 𝑇=𝑇^′ and (𝑇𝐯,𝐯)_𝑌>0 for all 0≠𝐯∈𝑌 and call 𝑇 symmetric positive definite (spd).
Let 𝐺=𝐺(𝛿)>0 be a linear, bounded, compact operator on 𝑋 representing a generic smoothing filter with filter radius 𝛿>0: 𝐺∶𝐿²(Ω)→𝐿²(Ω),𝜙∶=𝐺𝜙. (3.1) One example of this operator is the continuous differential filter 𝐺=𝐴⁻¹=(−𝛿²Δ+𝐼)⁻¹ (Definition 3.2), which is used, together with its discrete counterpart (𝐴^ℎ)⁻¹ (Definition 3.3), for implementation of our numerical scheme (Section 5.2).

Definition 3.2 (continuous differential filter). Fix 𝜙∈𝐿²(Ω). Then 𝜙∈𝑋 is the unique solution of −𝛿²Δ𝜙+𝜙=𝜙 with corresponding weak formulation 𝛿²∇𝜙,∇𝐯+𝜙,𝐯=𝜙,𝐯,∀𝐯∈𝑋. (3.2) Set 𝐴=−𝛿²Δ+𝐼 so that 𝐴⁻¹∶𝐿²(Ω)→𝑋, defined by 𝜙∶=𝐴⁻¹𝜙, is well defined.

Definition 3.3 (discrete differential filter). Fix 𝜙∈𝐿²(Ω). Then 𝜙^ℎ∈𝑋^ℎ is the unique solution of the following 𝛿²∇𝜙^ℎ,∇𝐯^ℎ+𝜙^ℎ,𝐯^ℎ=𝜙,𝐯^ℎ,∀𝐯^ℎ∈𝑋^ℎ. (3.3) Set 𝐴^ℎ=−𝛿²Δ^ℎ+Π^ℎ so that (𝐴^ℎ)⁻¹∶𝐿²(Ω)→𝑋^ℎ, defined by 𝜙^ℎ∶=(𝐴^ℎ)⁻¹𝜙, is well defined. Here, Π^ℎ∶𝐿²(Ω)→𝑋^ℎ is the 𝐿² projection and Δ^ℎ∶𝑋→𝑋^ℎ the discrete Laplace operator satisfying the following Π^ℎ𝜙−𝜙,𝐯^ℎ=0,Δ^ℎ𝜙,𝐯^ℎ=−∇𝜙,∇𝐯^ℎ∀𝐯^ℎ∈𝑋^ℎ. (3.4)
It is well known that 𝐴⁻¹ and (𝐴^ℎ)⁻¹ are each linear and bounded, 𝐴⁻¹ is compact, and the spectrum of 𝐴 and 𝐴^ℎ (on 𝑋 and 𝑋^ℎ, resp.) is contained in [1,∞) and spectrum of 𝐴⁻¹ and (𝐴^ℎ)⁻¹ (on 𝑋 and 𝑋^ℎ, resp.) is contained in (0,1] so that 𝐴⁻¹>0on𝑋,𝐴^ℎ−1>0on𝑋^ℎ. (3.5) For more detailed exposition on these operators, see [13].

3.1. A Family of Deconvolution Operators

We analyze (1.6) for stable, accurate deconvolution 𝐷 of the smoothing filter 𝐺 introduced in Section 3 so that 𝐷𝐺𝐮 accurately approximates the smooth parts of 𝐮.

Assumption 3.4 (continuous deconvolution operator). Suppose that 𝐷∶𝑋→𝑋 is linear, bounded, spd, and commutes with 𝐺 so that ‖𝐷𝐺‖≤1,|𝐷𝐯|₁≤𝑑₁|𝐯|₁∀𝐯∈𝑋, (3.6) for some constant 𝑑₁>0. Moreover, suppose that 𝐷 is parametrized by 𝛼>0, 𝛿>0 so that ‖(𝐼−𝐷𝐺)𝐯‖≤𝑐₁(𝛿,𝛼)‖Δ𝐯‖→0,as𝛼,𝛿→0 (3.7) for smooth enough 𝐯∈𝑋.
Note that the first estimate in (3.6) is required so that the spectral radius satisfies 𝜌(𝐷𝐺)≤1. The second estimate in (3.6) (which controls the 𝐻¹-seminorm of 𝐷𝐺) is required for the convergence analysis in Section 4.2.
Assumption 3.5 prescribes properties of the discrete analogue 𝐷^ℎ∶𝑋^ℎ→𝑋^ℎ corresponding to the continuous deconvolution operator 𝐷∶𝑋→𝑋 (Assumptions 3.4).

Assumption 3.5 (discrete deconvolution operator). Let 𝐷 satisfy Assumption 3.4. Let 𝐺^ℎ∶𝑋^ℎ→𝑋^ℎ be a discrete analogue of 𝐺 that is linear, bounded, spd. Suppose that 𝐷^ℎ∶𝑋^ℎ→𝑋^ℎ is linear, bounded, spd, and commutes with 𝐺^ℎ such that 𝐷^ℎ𝐺^ℎ≤1,𝐷^ℎ𝐯^ℎ₁≤𝑑₁|𝐯|₁∀𝐯∈𝑋, (3.8) for some constant 𝑑₁>0. Moreover, suppose that 𝐷^ℎ is parametrized by 𝛼>0,𝛿>0,ℎ such that 𝐷=𝐷(ℎ,𝛿,𝛼) and 𝐷𝐺−𝐷^ℎ𝐺^ℎ𝐯≤𝑐₂(ℎ,𝛿,𝛼)‖𝐯‖_𝑘+1→0,asℎ,𝛿,𝛼→0, (3.9) for all 𝐯∈𝑋∩𝐻^𝑘+1(Ω) for some 𝑘≥0.

The estimates in (3.8) are motivated by the continuous case of (3.6). The approximation (3.9) is required for the convergence analysis in Section 4.2 (see Theorem 4.6, Corollary 4.7).

Remark 3.6. If 𝐷=𝑓(𝐺) for some continuous map 𝑓∶ℝ→ℝ, then commutativity is satisfied 𝐷𝐺=𝐺𝐷. Tikhonov-Lavrentiev (modified) regularization with 𝐺=𝐴⁻¹, 𝐺^ℎ=(𝐴^ℎ)⁻¹ given by 𝐷=((1−𝛼)𝐴⁻¹+𝛼𝐼)⁻¹, 𝐷^ℎ=((1−𝛼)(𝐴^ℎ)⁻¹+𝛼Π^ℎ)⁻¹ is one such example with 𝑓(𝑥)=((1−𝛼)𝑥+𝛼)⁻¹ and 𝑑₁=1, 𝑐₁=𝛼𝛿², 𝑐₂=𝛼𝛿²ℎ^𝑘+ℎ^𝑘+1, see [13].

Remark 3.7. Letting 𝜆_𝑘(⋅) denote the 𝑘th (ordered) eigenvalue of a given operator, commutativity of 𝐷 and 𝐺 provides 𝜆_𝑘(𝐷𝐺)=𝜆_𝑘(𝐷)𝜆_𝑘(𝐺) and similarly for the discrete operator 𝐷^ℎ𝐺^ℎ.

We next derive several important consequences of 𝐷 and 𝐷^ℎ under Assumptions 3.4, 3.5 required in the forthcoming analysis.

Lemma 3.8. Suppose that 𝐺, 𝐺^ℎ, 𝐷, 𝐷^ℎ satisfy Assumptions 3.4, 3.5. Then, ‖𝐷𝐯‖≤‖𝐯‖,𝐷^ℎ𝐯^ℎ≤‖𝐯‖∀𝐯∈𝐿²(Ω). (3.10)

Proof. For the continuous operator, ‖𝐷𝐯‖≤‖𝐷𝐺‖‖𝐯‖. (3.11) Then, (3.10)(a) follows from Assumption 3.4, and (3.10)(b) is derived similarly applying Assumption 3.5 instead.

Lemma 3.9. Suppose that 𝐺, 𝐺^ℎ, 𝐷, 𝐷^ℎ satisfy Assumptions 3.4, 3.5. Then, the spectrum of both 𝐷𝐺 and 𝐷^ℎ𝐺^ℎ are contained in [0,1] so that ‖𝐼−𝐷𝐺‖≤1,𝐼−𝐷^ℎ𝐺^ℎ≤1. (3.12) As a consequence, ‖𝐯‖²_⋆∶=(𝐯−𝐷)‖𝐯‖𝐯,𝐯_2⋆ℎ∶=𝐯−𝐷^ℎ𝐯^ℎ,𝐯≥0,∀𝐯∈𝐿²(Ω). (3.13)

Proof. Assumptions 3.4, 3.5 guarantee that the spectral radius 𝜌(𝐷𝐺)≤1 and 𝜌(𝐷^ℎ𝐺^ℎ)≤1. Also, 𝐷>0 and 𝐺>0 and commute so that 𝐷𝐺≥0. Similarly, 𝐷^ℎ𝐺^ℎ≥0. Therefore, the spectrum of 𝐷𝐺, 𝐷^ℎ𝐺^ℎ≥0 are each contained in [0,1]. So, 𝐼−𝐷𝐺, 𝐼−𝐷^ℎ𝐺^ℎ≥0 have spectrum contained in [0,1] which ensures the non-negativity of both ||⋅||_⋆ and ||⋅||_⋆ℎ.

3.2. Iterated Deconvolution

One can show, by eliminating intermediate steps in the definition of the iterated regularization operator 𝐷_𝑗 in (1.9) with base operator 𝐷 satisfying Assumption 3.4, that 𝐷_𝑗=𝐷^𝑗_𝑖=0(𝐹𝐷)^𝑖,𝐹∶=𝐷⁻¹−𝐺. (3.14) Similarly, the discrete iterated regularization operator 𝐷^ℎ_𝑗 with discrete base operator 𝐷^ℎ satisfying Assumption 3.5, is given by the following 𝐷^ℎ_𝑗=𝐷^ℎ𝑗_𝑖=0𝐹^ℎ𝐷^ℎ𝑖,𝐹^ℎ∶=𝐷^ℎ−1−𝐺^ℎ. (3.15) We next show that 𝐷_𝑗 and 𝐷^ℎ_𝑗 for 𝑗>0 inherit several important properties from 𝐷 and 𝐷^ℎ, respectively, via Assumption 3.5.

Proposition 3.10. Fix 𝑗∈ℕ. Then 𝐷_𝑗∶𝑋→𝑋 defined by (3.14) satisfies Assumption 3.4. In particular, 𝐷_𝑗>0 is linear, bounded, commutes with 𝐺 and satisfies (3.6)(a). Estimate (3.6)(b) is replaced by the following 𝐷_𝑗𝐯₁≤𝑑_1,𝑗|𝐯|₁∀𝐯∈𝑋 (3.16) for some constant 𝑑_1,𝑗>0. Estimate (3.7) is replaced by the following 𝐼−𝐷_𝑗𝐺𝐯≤𝑐_1,𝑗(𝛿,𝛼)‖Δ𝐯‖²→0,as𝛼,𝛿→0. (3.17) Moreover, 𝑑_1,𝑗≤∑^𝑗_𝑖=0𝑑^𝑖₁ and 𝑐_1,𝑗≤∑^𝑗_𝑖=0𝑐^𝑖₁.

Proof. First notice that 𝐷_𝑗 is linear and bounded since it is a linear combination of linear and bounded operators 𝐷(𝐹𝐷)^𝑖=𝐷(𝐼−𝐷𝐺)^𝑖, for 𝑖=0,1,…,𝑗. Moreover, since 𝐺 commutes with 𝐷, it follows that 𝐺 commutes with 𝐷(𝐼−𝐷𝐺)^𝑖 and hence with 𝐷_𝑗. Next, 𝐷_𝑗 is a sum of spd and snn operators 𝐷>0, 𝐷(𝐼−𝐷𝐺)^𝑖≥0. Hence, 𝐷_𝑗>0. Next, notice that 𝐷_𝑗𝐺=^𝑗_𝑖=0(𝐼−𝐷𝐺)^𝑖𝐷𝐺=𝐼+(𝐼−𝐷𝐺)+⋯+(𝐼−𝐷𝐺)^𝑗𝐷𝐺. (3.18) Letting 𝜆_𝑘(⋅) denote the 𝑘th (ordered) eigenvalue of a given operator, we can characterize the spectrum of 𝐷_𝑗 by summing the resulting finite geometric series (3.18) to get 𝜆_𝑘𝐷_𝑗𝐺=𝜆_𝑘(𝐷)𝜆_𝑘(𝐺)^𝑗_𝑖=01−𝜆_𝑘(𝐷𝐺)^𝑖=1−1−𝜆_𝑘(𝐷𝐺)^𝑗+1. (3.19) Then under Assumption 3.4, Lemma 3.9 with (3.19) implies that 0≤𝜆_𝑘(𝐷_𝑗𝐺)≤||𝐷_𝑗𝐺||≤1. Hence, 𝐷_𝑗 satisfies (3.6)(a). Expanding the terms in (3.18) as powers of 𝐷𝐺, we see that (3.18) can be written as a polynomial (with coefficients 𝑎_𝑖) in 𝐷𝐺, so that ∇𝐷_𝑗𝐯=^𝑗_𝑖=0𝑎_𝑖∇(𝐷𝐺)^𝑖𝐯,𝐷_𝑗𝐯₁≤^𝑗_𝑖=0𝑑^𝑖₁|𝐯|₁, (3.20) since |𝐷𝐯|₁≤𝑑₁|𝐯|₁ can be applied successfully. Therefore (3.16) follows with 𝑑_1,𝑗=∑^𝑗_𝑖=0𝑑^𝑖₁. Next, start with (3.18) to get 𝐼−𝐷_𝑗𝐺𝐯=(𝐼−𝐷𝐺)𝐯+𝐷𝐺(𝐼−𝐷𝐺)𝐯+⋯+𝐷𝐺(𝐼−𝐷𝐺)^𝑗𝐯≤‖(𝐼−𝐷𝐺)𝐯‖+‖𝐷𝐺‖‖(𝐼−𝐷𝐺)𝐯‖+⋯+‖𝐷𝐺‖‖𝐼−𝐷𝐺‖^𝑗−1‖(𝐼−𝐷𝐺)𝐯‖. (3.21) Estimate (3.17) follows by noting ||𝐷𝐺||≤1, ||𝐼−𝐷𝐺||≤1, and by Assumption 3.5, ||(𝐼−𝐷𝐺)𝐯||≤𝑐₁||Δ𝐯||.

Proposition 3.11. Fix 𝑗∈ℕ. Then 𝐷^ℎ_𝑗∶𝑋^ℎ→𝑋^ℎ defined by (3.15) satisfies Assumption 3.5. In particular, 𝐷^ℎ_𝑗>0 is linear, bounded, commutes with 𝐺^ℎ and satisfies (3.8)(a). Estimate (3.8)(b) is replaced by the following 𝐷^ℎ_𝑗𝐯^ℎ₁≤𝑑_1,𝑗|𝐯|₁∀𝐯∈𝑋 (3.22) for some constant 𝑑_1,𝑗>0. Estimate (3.9) is replaced by the following 𝐷_j𝐺−𝐷^ℎ_𝑗𝐺^ℎ𝐯≤𝑐_2,𝑗(ℎ,𝛿,𝛼)‖𝐯‖_𝑘+1→0,asℎ,𝛿,𝛼→0 (3.23) for any 𝐯∈𝑋∩𝐻^𝑘+1(Ω) for some 𝑘≥0. Moreover, 𝑐_2,𝑗≤𝛽(𝑗)𝑐₂ for some constant 𝛽=𝛽(𝑗)>0.

Proof. The first two assertions follow similarly as in the previous proof of Proposition 3.10. To prove (3.23), we start by writing 𝐷^ℎ_𝑗𝐺^ℎ=^𝑗_𝑖=0𝐼−𝐷^ℎ𝐺^ℎ𝑖𝐷^ℎ𝐺^ℎ, (3.24) and then subtract (3.24) from (3.18) to get 𝐷_𝑗𝐺−𝐷^ℎ_𝑗𝐺^ℎ=Λ_𝑗𝐷𝐺−𝐷^ℎ𝐺^ℎ+Λ_𝑗−Λ^ℎ_𝑗𝐷^ℎ𝐺^ℎ, (3.25) where Λ_𝑗=^𝑗_𝑖=0(𝐼−𝐷𝐺)^𝑖,Λ^ℎ_𝑗=^𝑗_𝑖=0𝐼−𝐷^ℎ𝐺^ℎ𝑖. (3.26) Then taking norms across (3.25), we get 𝐷_𝑗𝐺−𝐷^ℎ_𝑗𝐺^ℎ𝐯=Λ_𝑗𝐷𝐺−𝐷^ℎ𝐺^ℎ𝐯+𝐷^ℎ𝐺^ℎΛ_𝑗−Λ^ℎ_𝑗𝐯. (3.27) Notice that ||𝐼−𝐷𝐺||≤1 so that ||Λ_𝑗||≤𝑗+1. Moreover, ||(𝐷𝐺−𝐷^ℎ𝐺^ℎ)𝐯||≤𝑐₂||𝐯||_𝑘+1 via Assumption 3.5. Next, using the binomial theorem and factoring, we get Λ_𝑗−Λ^ℎ_𝑗𝐯=^𝑗_𝑖=0𝑗!!()𝑖!𝑗−𝑖−1^𝑖(𝐷𝐺)^𝑖−𝐷^ℎ𝐺^ℎ𝑖𝐯=^𝑗_𝑖=0𝑗!!()𝑖!𝑗−𝑖−1^𝑖𝑖_𝑛=0(𝐷𝐺)^𝑛𝐷^ℎ𝐺^ℎ_𝑛−𝑖𝐷𝐺−𝐷^ℎ𝐺^ℎ𝐯. (3.28) Then, applying ||𝐷𝐺||≤1, ||𝐷^ℎ𝐺^ℎ||≤1 to (3.28) provides Λ_𝑗−Λ^ℎ_𝑗𝐯=^𝑗_𝑖=0𝑗!𝑖(𝑖)!𝑗−𝑖!𝐷𝐺−𝐷^ℎ𝐺^ℎ𝐯. (3.29) Again, ||(𝐷𝐺−𝐷^ℎ𝐺^ℎ)||≤𝑐₂||𝐯||_𝑘+1 via Assumption 3.5. So, we combine these above results to conclude (3.23) with 𝛽(𝑗)=∑^𝑗_𝑖=0𝑗!𝑖/(𝑖)!(𝑗−𝑖)!.

3.3. Tikhonov-Lavrentiev Regularization

We provide two examples of discrete deconvolution operators 𝐷^ℎ to make the abstract formulation in the previous section more concrete. The Tikhonov-Lavrentiev and modified Tikhonov-Lavrentiev operator (for linear, compact 𝐺>0) is given by the following Tikhonov-Lavrentiev⇒𝐷_𝛼,0=(𝐺+𝛼𝐼)⁻¹modifiedTikhonov-Lavrentiev⇒𝐷_𝛼,0=((1−𝛼)𝐺+𝛼𝐼)⁻¹. (3.30)

Definition 3.12 ((weak) modified Tikhonov-Lavrentiev deconvolution). Fix 𝛼>0. Let 𝐺=𝐴⁻¹. For any 𝐰∈𝑋, let 𝝎₀∶=𝐷_𝛼,0𝐰∈𝑋 be the unique solution of 𝛼𝛿²∇𝝎₀,∇𝐯+𝝎₀,𝐯=(𝐰,𝐯),∀𝐯∈𝑋. (3.31)

Definition 3.13 ((discrete) modified Tikhonov-Lavrentiev deconvolution). Fix 𝛼>0. Let 𝐺^ℎ=(𝐴^ℎ)⁻¹ and 𝐷^ℎ_𝛼,0=((1−𝛼)(𝐴^ℎ)⁻¹+𝛼Π^ℎ)⁻¹. For any 𝐰∈𝑋, let 𝝎^ℎ₀∶=𝐷^ℎ_𝛼,0𝐰^ℎ∈𝑋^ℎ be the unique solution of 𝛼𝛿²∇𝝎^ℎ₀,∇𝐯^ℎ+𝝎^ℎ₀,𝐯^ℎ=𝐰,𝐯^ℎ,∀𝐯^ℎ∈𝑋^ℎ. (3.32)
The iterated modified Tikhonov-Lavrentiev operator (for linear, compact 𝐺>0) is obtained from the Tikhonov-Lavrentiev operator with updates via (1.9): IteratedTikhonov-Lavrentiev⇒𝐷_𝛼,𝑗=𝐷^𝑗_𝛼,0𝑖=0𝛼𝐷_𝛼,0^𝑖,IteratedmodifiedTikhonov-Lavrentiev⇒𝐷_𝛼,𝑗=𝐷^𝑗_𝛼,0𝑖=0𝛼(𝐼−𝐺)𝐷_𝛼,0^𝑖. (3.33)

Definition 3.14 (iterated modified Tikhonov-Lavrentiev deconvolution (weak)). Fix 𝛼>0 and 𝐽∈ℕ. Let 𝐺=𝐴⁻¹. Define 𝝎₋₁=0, then for any 𝐰∈𝑋 and 𝑗=0,1,…,𝐽, let 𝝎_𝑗∶=𝐷_𝛼,𝑗𝐰∈𝑋 be the unique solution of 𝛼𝛿²∇𝝎_𝑗,∇𝐯+𝝎_𝑗,𝐯=(𝐰,𝐯)+𝛼𝛿²∇𝝎_𝑗−1,∇𝐯,∀𝐯∈𝑋. (3.34)

Definition 3.15 (iterated modified Tikhonov-Lavrentiev deconvolution (discrete)). Fix 𝛼>0 and 𝐽∈ℕ. Let 𝐺^ℎ=(𝐴^ℎ)⁻¹, and 𝐷^ℎ_𝛼,𝑗=𝐷^ℎ_𝛼,0∑^𝑗_𝑖=0(𝛼(Π^ℎ−(𝐴^ℎ)⁻¹)𝐷^ℎ_𝛼,0)^𝑖. Define 𝝎^ℎ₋₁=0, then for any 𝐰∈𝑋 and 𝑗=0,1,…,𝐽, let 𝝎^ℎ_𝑗∶=𝐷^ℎ_𝛼,𝑗𝐰^ℎ∈𝑋^ℎ be the unique solution of 𝛼𝛿²∇𝝎^ℎ_𝑗,∇𝐯^ℎ+𝝎^ℎ_𝑗,𝐯^ℎ=𝐰,𝐯^ℎ+𝛼𝛿²∇𝝎^ℎ_𝑗−1,∇𝐯^ℎ,∀𝐯^ℎ∈𝑋^ℎ. (3.35)

4. Well Posedness of the Fully Discrete Model

We now state the proposed algorithm.

Problem 1 (CNFE for Leray-deconvolution). Let (𝐰₀,𝜋₀)∈(𝑋^ℎ,𝑄^ℎ). Then, for each 𝑛=0,1,…,𝑀−1, find (𝐰^ℎ_𝑛+1,𝜋^ℎ_𝑛+1)∈(𝑋^ℎ,𝑄^ℎ) satisfying 1𝐰Δ𝑡^ℎ𝑛+1−𝐰^ℎ_𝑛,𝐯^ℎ+𝑏^ℎ𝜑^ℎ𝐰^ℎ_𝑛+1/2,𝐰^ℎ_𝑛+1/2,𝐯^ℎ−𝜋^ℎ_𝑛+1/2,∇⋅𝐯^ℎ+Re⁻¹∇𝐰^ℎ_𝑛+1/2,∇𝐯^ℎ+𝜒𝐰^ℎ_𝑛+1/2−𝜑^ℎ𝐰^ℎ_𝑛+1/2,𝐯^ℎ=𝐟_𝑛+1/2,𝐯^ℎ,∀𝐯^ℎ∈𝑋^ℎ,∇⋅𝐰^ℎ_𝑛+1,𝑞^ℎ=0,∀𝑞^ℎ∈𝑄^ℎ, (4.1)(4.2) where 𝜑^ℎ(𝐰^ℎ_𝑛+1/2)=𝐷^ℎ𝐰^ℎℎ_𝑛+1/2.

Notice that (𝑞^ℎ_𝑛+1/2,∇⋅𝐯^ℎ)=0 when 𝐯^ℎ∈𝐕^ℎ so that the problem of finding 𝐰^ℎ_𝑛+1∈𝑉^ℎ satisfying 1𝐰Δ𝑡^ℎ𝑛+1−𝐰^ℎ_𝑛,𝐯^ℎ+𝑏^ℎ𝜑^ℎ𝐰^ℎ_𝑛+1/2,𝐰^ℎ_𝑛+1/2,𝐯^ℎ+Re⁻¹∇𝐰^ℎ_𝑛+1/2,∇𝐯^ℎ+𝜒𝐰^ℎ_𝑛+1/2−𝜑^ℎ𝐰^ℎ_𝑛+1/2,𝐯^ℎ=𝐟_𝑛+1/2,𝐯^ℎ,∀𝐯^ℎ∈𝑉^ℎ. (4.3)

4.1. Well Posedness

We establish existence of 𝐰 at each time step of (4.3) by Leray-Schauder’s fixed-point theorem.

Lemma 4.1. Let 𝑎𝜽^ℎ,𝐯^ℎ=Δ𝑡2Re∇𝜽^ℎ,∇𝐯^ℎ+𝜒Δ𝑡2𝜽^ℎ−𝜑^ℎ𝜽^ℎ,𝐯^ℎ,𝑙_𝑦𝐯^ℎ=𝐲,𝐯^ℎ. (4.4) for any 𝐲∈𝑋^′ and 𝜽^ℎ, 𝐯^ℎ∈𝑉^ℎ. Suppose that 𝐷^ℎ satisfies Assumption 3.5. Then 𝑎(⋅,⋅)∶𝑉^ℎ×𝑉^ℎ→ℝ is a continuous and coercive bilinear form and 𝑙_𝑦(⋅)∶𝑉^ℎ→ℝ is a linear, continuous functional.

Proof. Linearity for 𝑙_𝑦(⋅) is obvious, and continuity follows from an application of Hölder’s inequality. Continuity for 𝑎(⋅,⋅) also follows from Hölder’s inequality and Assumption 3.5. Coercivity is proven by application of (3.13).

Lemma 4.2. Let 𝑇∶𝑋^′→𝑉^ℎ be such that, for any 𝐲∈𝑋^′, 𝜽^ℎ∶=𝑇(𝐲) solves 𝑎𝜽^ℎ,𝐯^ℎ=𝑙_𝑦𝐯^ℎ,∀𝐯^ℎ∈𝑉^ℎ. (4.5) Then 𝑇 is a well-defined, linear, bounded operator.

Proof. Linearity is clear. The results of Lemma 4.1, and the Lax-Milgram theorem prove the rest.

Lemma 4.3. Fix 𝑛=0,1,…,𝑀−1. Let 𝐰^ℎ_𝑛 be a solution of Problem 1 and let 𝑁∶𝑉^ℎ→𝑋^′ satisfy, for any 𝜽^ℎ∈𝑉^ℎ, 𝑁𝜽^ℎ,𝐯^ℎ=−𝜽^ℎ−2𝐰^ℎ_𝑛,𝐯^ℎ−Δ𝑡4𝑏^ℎ𝜑^ℎ𝜽^ℎ,𝜽^ℎ,𝐯^ℎ+Δ𝑡𝐟_𝑛+1/2,𝐯^ℎ=∶𝑐𝜽^ℎ,𝐯^ℎ,∀𝐯^ℎ∈𝑉^ℎ. (4.6) Then 𝑁(𝜽^ℎ) is well-defined, bounded, and continuous.

Proof. For each 𝜽^ℎ∈𝑉^ℎ, the map 𝐯^ℎ∈𝑉^ℎ↦𝑐(𝜽^ℎ,𝐯^ℎ) is a bounded, linear functional (apply Hölder’s inequality and (2.11)). Since 𝑉^ℎ is a Hilbert space, we conclude that 𝑁(𝜽^ℎ) is well defined, by the Riesz-Representation theorem. Moreover, 𝑁(𝜽^ℎ) is bounded on 𝑉^ℎ and since the underlying function space is finite dimensional, continuity follows.

Lemma 4.4. Fix 𝑛∈ℕ. Let 𝐹∶𝑉^ℎ→𝑉^ℎ be defined such that 𝐹(𝜽^ℎ)=(𝑇∘𝑁)(𝜽^ℎ). Then, 𝐹 is a compact operator.

Proof. 𝑁(⋅) is a compact operator (continuous on a finite dimensional function space). Thus, 𝐹 is a continuous composition of a compact operator and hence compact itself.

Theorem 4.5 (well posedness). Fix 𝑛=0,1,2,…,𝑀−1<∞. There exists (𝐰^ℎ_𝑛,𝜋^ℎ_𝑛)∈𝑋^ℎ×𝑄^ℎ satisfying Problem 1. Moreover, 𝐰^ℎ_𝑚²+12ReΔ𝑡^𝑚−1_𝑛=0𝐰^ℎ_𝑛+1/2²₁+𝜒Δ𝑡^𝑚−1_𝑛=0𝐰^ℎ_𝑛+1/2²_⋆ℎ≤𝐰^ℎ₀²+Δ𝑡Re2^𝑚−1_𝑛=0𝐟_𝑛+1/2²₋₁, (4.7) for all integers 1≤𝑚≤𝑀, independent of Δ𝑡>0.

Proof. First, assume that (𝐰^ℎ_𝑛+1,𝑞^ℎ_𝑛+1) is a solution to (4.1), (4.2). Set 𝐯^ℎ=𝐰^ℎ_𝑛+1/2 in (4.1) so that skew-symmetry of the nonlinear term provides 1𝐰2Δ𝑡^ℎ𝑛+12−𝐰^ℎ_𝑚²+Re⁻¹𝐰^ℎ_𝑛+1/2²₁+𝜒𝐰^ℎ_𝑛+1/2−𝜑^ℎ𝐰^ℎ_𝑛+1/2,𝐰^ℎ_𝑛+1/2=𝐟_𝑛+1/2,𝐰^ℎ_𝑛+1/2. (4.8) Duality of 𝑋^′×𝑋 with Young’s inequality implies 𝐟_𝑛+1/2,𝐰^ℎ_𝑛+1/2≤Re2𝐟_𝑛+1/2²₋₁+Re2𝐰^ℎ_𝑛+1/2²₁. (4.9) From (3.13), we have 𝐰^ℎ_𝑛+1/2²_⋆ℎ=𝐰^ℎ_𝑛+1/2−𝜑^ℎ𝐰^ℎ_𝑛+1/2,𝐰^ℎ_𝑛+1/2≥0. (4.10) Then applying (4.9), (4.10) to (4.8), combining-like terms and simplifying provides 1𝐰2Δ𝑡^ℎ𝑛+12−𝐰^ℎ_𝑛²+Re2𝐰^ℎ_𝑛+1/2²₁+𝜒𝐰^ℎ_𝑛+1/2²_⋆ℎ≤Re2𝐟_𝑛+1/2²₋₁. (4.11) Summing from 𝑛=0 to 𝑚−1, we get the desired bound.
Next, let 𝐖^ℎ_𝑛=𝐰^ℎ_𝑛+1+𝐰^ℎ_𝑛. Showing that 𝐖^ℎ_𝑛=𝐹(𝐖^ℎ_𝑛) has a fixed point will ensure existence of solutions to (4.3). Indeed, if we can show that 𝐖^ℎ₀=𝐹(𝐖^ℎ₀), then since 𝐰^ℎ₀ is given initial data, existence of 𝐰^ℎ₁ is immediate. Induction can be applied to prove existence of (𝐰^ℎ_𝑛)_{1≤𝑛≤𝑀}. To this end, since 𝐹 is compact, it is enough to show (via Leray Schauder) that any solution 𝐖^ℎ_𝑛,𝜆 of the fixed-point problem 𝐖^ℎ_𝑛,𝜆=𝜆𝐹(𝐖^ℎ_𝑛,𝜆) is uniformly bounded with respect to 0≤𝜆≤1. Hence, we consider 𝑎𝐖^ℎ_𝑛,𝜆,𝐯^ℎ=𝜆𝑁𝐖^ℎ_𝑛,𝜆,𝐯^ℎ. (4.12)
Test with 𝐯^ℎ=𝐖^ℎ_𝑛,𝜆, use skew-symmetry of the trilinear form and properties of 𝐷^ℎ given in Assumption 3.5 and (3.13) to get 𝜆𝐖^ℎ_𝑛,𝜆²+Δ𝑡𝐖2Re^{ℎ𝑛,𝜆21}+𝜒Δ𝑡2𝐖^ℎ_𝑛,𝜆²_⋆ℎ≤2𝜆𝐰^ℎ_𝑛,𝐖^ℎ_𝑛,𝜆+𝜆Δ𝑡𝐟_𝑛+1/2,𝐖^ℎ_𝑛,𝜆. (4.13) Duality of 𝑋^′×𝑋 followed by Young’s inequality implies 𝜆Δ𝑡𝐟_𝑛+1/2,𝐖^ℎ_𝑛,𝜆≤Δ𝑡Re𝐟_𝑛+1/2²₋₁+Δ𝑡𝐖4Re^{ℎ𝑛,𝜆21}. (4.14) Since 𝐰^ℎ_𝑛∈𝐿²(Ω) from the a priori estimate (4.7), we apply Hölder’s and Young’s inequalities to get 2𝜆𝐰^ℎ_𝑛,𝐖^ℎ_𝑛,𝜆≤2𝐰^ℎ_𝑛²+𝜆2𝐖^ℎ_𝑛,𝜆². (4.15) Applying estimates (4.14), (4.15) to (4.13) we get that |𝐖^ℎ_𝑛,𝜆|₁≤𝐶<∞ independent of 𝜆. By the Leray-Schauder fixed point theorem, given 𝐰^ℎ_𝑛, there exists a solution to the fixed-point theorem 𝐖^ℎ_𝑛=𝐹(𝐖^ℎ_𝑛). By the induction argument noted above, there exists a solution 𝐰^ℎ_𝑛 for each 𝑛=0,1,2,…,𝑀−1 to (4.3). Existence of an associated discrete pressure follows by a classical argument, since the pair (𝑋^ℎ,𝑄^ℎ) satisfies the discrete inf-sup condition (2.15).

4.2. Convergence Analysis

Under usual regularity assumptions, we summarize the main convergence estimate in Theorem 4.6. Suppose that 𝐷 represents deconvolution with 𝐽-updates.

Theorem 4.6. Suppose that (𝐮,𝑝) are strong solutions to (2.3), (2.4), (2.5) and that 𝐺, 𝐺^ℎ, 𝐷, 𝐷^ℎ satisfy Assumptions 3.4, 3.5. Suppose further that 𝐮∈𝑙²(𝐻²(Ω)∩𝑉)∩𝑙^∞(𝑋), 𝐮_𝑡∈𝑙²(𝑋^′)∩𝐿²(𝑋^′), Δ^𝛽+1𝐮∈𝐿²(Ω) for some 0≤𝛽≤𝐽+1, 𝑝∈𝑙²(𝑄). If 𝐶ReΔ𝑡𝐮_𝑛²₂<1,∀𝑛=0,1,…,𝑀 (4.16) then, 𝐮_𝑀−𝐰^ℎ_𝑀²+Re⁻¹Δ𝑡^𝑀−1_𝑛=0𝐮_𝑛+1/2−𝐰^ℎ_𝑛+1/2²₁≤𝐶𝐮₀−𝐰^ℎ₀²+𝐸+𝑝²_𝑙²_(𝐿2(Ω))+‖∇𝐮‖²_𝑙^∞_(𝐿2(Ω))+𝜒²‖𝐮‖²_{𝑙2(𝐻1(Ω))}+‖𝐮‖²_𝑙^∞_(𝐿2(Ω))+‖𝐮‖²_{𝑙2(𝐿2(Ω))}+𝜒²+‖𝐮‖²_𝑙^∞_(𝐻2(Ω))𝑐²_1,𝐽Δ^𝛽𝐮²+𝑐²_2,𝑗‖∇𝐮‖²_{𝑙2(𝐿2(Ω))}, (4.17) where 𝐸>0 is given in (4.52).

Corollary 4.7 (convergence estimate). Under the assumptions of Theorem 4.6, suppose further that (𝐮,𝑝) satisfy the assumptions for (2.16) for some 𝑘≥1 and 𝑠≥0, 𝐮∈𝑙^∞(𝐻^𝑘(Ω))∩𝑙²(𝐻^𝑘+1(Ω)), 𝐮_𝑡∈𝐿²(𝐻^𝑘−1(Ω)∩𝐻¹(Ω))∩𝑙^∞(𝐻¹(Ω)), 𝐮_𝑡𝑡∈𝐿²(𝐿²(Ω)), 𝐮_𝑡𝑡𝑡∈𝐿²(𝑋^′), and 𝑝∈𝑙²(𝐻^𝑠+1(Ω)). If ||𝐮₀−𝐰^ℎ₀||≤𝐶(ℎ^𝑘+(𝛼𝛿²)^𝛽) then 𝐮_𝑀−𝐰^ℎ_𝑀²+Re⁻¹Δ𝑡^𝑀−1_𝑛=0𝐮_𝑛+1/2−𝐰^ℎ_𝑛+1/2²₁≤𝐶ℎ^2𝑘+ℎ^2𝑠+2+Δ𝑡⁴+𝑐_1,𝐽(𝛿,𝛼)²+𝑐_2,𝐽(ℎ,𝛿,𝛼)². (4.18)

Proof of Theorem 4.6, Corollary 4.7. Suppose that 𝐮, 𝑝 satisfying (2.3), (2.4) also satisfy 𝐮∈𝐶⁰(𝑉), 𝐮_𝑡∈𝐶⁰(𝑋^′), and 𝑝∈𝐶⁰(𝑄) so that, for each 𝑛=𝑛₀,𝑛₀+1,…,𝑀−1, 𝐮_𝑡𝑛+1/2+1,𝐯2𝐮_𝑛+1⋅∇𝐮_𝑛+1+1,𝐯2𝐮_𝑛⋅∇𝐮_𝑛,𝐯+Re⁻¹∇𝐮_𝑛+1/2,∇𝐯−𝑝_𝑛+1/2=𝐟,∇⋅𝐯_𝑛+1/2,𝐯,∀𝐯∈𝑋,∇⋅𝐮_𝑛+1=0. (4.19) The consistency error for the time-discretization 𝜏_𝑛⁽¹⁾(𝐮,𝑝;𝐯^ℎ) and regularization/time-relaxation error 𝜏_𝑛⁽²⁾(𝐮,𝑝;𝐯^ℎ) are given by, for 𝑛=0,1,…,𝑀−1, 𝜏_𝑛⁽¹⁾𝐮,𝑝;𝐯^ℎ𝐮∶=^𝑛+1−𝐮_𝑛−𝐮Δ𝑡^𝑡𝑛+1/2,𝐯^ℎ+𝐮_𝑛+1/2⋅∇𝐮_𝑛+1/2,𝐯^ℎ−12𝐮_𝑛+1⋅∇𝐮_𝑛+1,𝐯^ℎ−12𝐮_𝑛⋅∇𝐮_𝑛,𝜏,𝐯_𝑛(2)𝐮,𝑝;𝐯^ℎ∶=−𝑏^ℎ𝐮_𝑛+1/2−𝜑^ℎ𝐮_𝑛+1/2,𝐮_𝑛+1/2,𝐯^ℎ+𝜒𝐮_𝑛+1/2−𝜑^ℎ𝐮_𝑛+1/2,𝐯^ℎ, (4.20) where 𝐯^ℎ∈𝑋^ℎ. Write 𝜏_𝑛∶=𝜏_𝑛⁽¹⁾+𝜏_𝑛⁽²⁾. Using (4.20), rewrite (4.19) in a form conducive to analyzing the error between the continuous and discrete models: 𝐮_𝑛+1−𝐮_𝑛Δ𝑡,𝐯^ℎ+𝑏^ℎ𝜑^ℎ𝐮_𝑛+1/2,𝐮_𝑛+1/2,𝐯^ℎ+Re⁻¹∇𝐮_𝑛+1/2,∇𝐯^ℎ−𝑝_𝑛+1/2,∇⋅𝐯^ℎ+𝜒𝐮_𝑛+1/2−𝜑^ℎ𝐮_𝑛+1/2,𝐯^ℎ=𝐟_𝑛+1/2,𝐯^ℎ+𝜏_𝑛𝐮,𝑝;𝐯^ℎ. (4.21) Let 𝐯^ℎ=𝐔^ℎ_𝑛 be the 𝐿²-projection of 𝐮(⋅,𝑡_𝑛) so that (𝜼_𝑛+1−𝜼_𝑛,𝜙^ℎ_𝑛+1/2)=0. Decompose the velocity error 𝐞_𝑛=𝐰^ℎ_𝑛−𝐮_𝑛=𝜙^ℎ_𝑛−𝜼_𝑛,𝜙^ℎ_𝑛=𝐰^ℎ_𝑛−𝐔^ℎ_𝑛,𝜼_𝑛=𝐮_𝑛−𝐔^ℎ_𝑛, (4.22) where 𝐔^ℎ_𝑛∈𝑉^ℎ. Fix 𝑞^ℎ_𝑛+1/2∈𝑄^ℎ. Note that (𝑞^ℎ_𝑛+1/2,∇⋅𝐯^ℎ)=0 for any 𝐯^ℎ∈𝑉^ℎ. Subtract (4.21) from (4.3), apply (3.13)(b), and test with 𝐯^ℎ=𝜙^ℎ_𝑛+1/2 to get 1𝜙2Δ𝑡^ℎ𝑛+12−𝜙^ℎ_𝑛²+Re⁻¹𝜙^ℎ_𝑛+1/2²₁+𝜒𝜙^ℎ_𝑛+1/2²_⋆ℎ=−𝑞^ℎ_𝑛+1/2−𝑝_𝑛+1/2,∇⋅𝜙^ℎ_𝑛+1/2−𝑏^ℎ𝜑^ℎ𝜙^ℎ_𝑛+1/2,𝐮_𝑛+1/2,𝜙^ℎ_𝑛+1/2+𝑏^ℎ𝜑^ℎ𝜼_𝑛+1/2,𝐮_𝑛+1/2,𝜙^ℎ_𝑛+1/2+𝑏^ℎ𝜑^ℎ𝐰^ℎ_𝑛+1/2,𝜼_𝑛+1/2,𝜙^ℎ_𝑛+1/2+Re⁻¹∇𝜼_𝑛+1/2,∇𝜙^ℎ_𝑛+1/2+𝜒𝜼_𝑛+1/2−𝜑^ℎ𝜼_𝑛+1/2,𝜙^ℎ_𝑛+1/2−𝜏_𝑛𝐮,𝑝;𝜙^ℎ_𝑛+1/2. (4.23) (Spatial discretization error): Fix 𝜀>0. First, apply Hölder’s and Young’s inequalities to get Re⁻¹∇𝜼_𝑛+1/2,∇𝜙^ℎ_𝑛+1/2+𝑝_𝑛+1/2−𝑞^ℎ_𝑛+1/2,∇⋅𝜙^ℎ_𝑛+1/2≤𝐶Re⁻¹𝜼_𝑛+1/2²₁+𝐶Re𝑝_𝑛+1/2−𝑞^ℎ_𝑛+1/2²+1𝜙𝜀Re^{ℎ𝑛+1/221}. (4.24) Apply (3.12) and duality estimate on 𝑋×𝑋^′ to get 𝜒𝜼_𝑛+1/2−𝜑^ℎ𝜼_𝑛+1/2,𝜙^ℎ_𝑛+1/2≤𝐶𝜒²Re𝜼_𝑛+1/2²₋₁+1𝜙𝜀Re^ℎ𝑛+1/21. (4.25) We bound the convective terms next. First, 𝐮∈𝐻²(Ω) and estimate (2.11)(b) give 𝑏^ℎ𝜑^ℎ𝜙^ℎ_𝑛+1/2,𝐮_𝑛+1/2,𝜙^ℎ_𝑛+1/2≤𝐶Re𝜑^ℎ𝜙^ℎ_𝑛+1/2²𝐮_𝑛+1/2²₂+1𝜙𝜀Re^{ℎ𝑛+1/221} (4.26) and 𝐮∈𝐿^∞(𝑋) with (2.11)(a) give 𝑏^ℎ𝜑^ℎ𝜼_𝑛+1/2,𝐮_𝑛+1/2,𝜙^ℎ_𝑛+1/2≤𝐶Re‖∇𝐮‖²_{𝐿∞(𝐿2(Ω))}𝜑^ℎ𝜼_𝑛+1/2²₁+1𝜙𝜀Re^{ℎ𝑛+1/221}. (4.27) Next, rewrite the remaining nonlinear term 𝑏^ℎ𝜑^ℎ𝐰^ℎ_𝑛+1/2,𝜼_𝑛+1/2,𝜙^ℎ_𝑛+1/2=𝑏^ℎ𝜑^ℎ𝐮_𝑛+1/2,𝜼_𝑛+1/2,𝜙^ℎ_𝑛+1/2−𝑏^ℎ𝜑^ℎ𝜼_𝑛+1/2,𝜼_𝑛+1/2,𝜙^ℎ_𝑛+1/2+𝑏^ℎ𝜑^ℎ𝜙^ℎ_𝑛+1/2,𝜼_𝑛+1/2,𝜙^ℎ_𝑛+1/2. (4.28) Once again, 𝐮∈𝐿^∞(𝑋) and (2.11)(a) give 𝑏^ℎ𝜑^ℎ𝐮_𝑛+1/2,𝜼_𝑛+1/2,𝜙^ℎ_𝑛+1/2+𝑏^ℎ𝜑^ℎ𝜼_𝑛+1/2,𝜼_𝑛+1/2,𝜙^ℎ_𝑛+1/2≤𝐶Re∇𝜑^ℎ(𝐮)_{𝐿∞(𝐿2(Ω))}𝜑^ℎ(𝐮)_{𝐿∞(𝐿2(Ω))}𝜼_𝑛+1/2²₁+𝐶Re𝜑^ℎ𝜼_𝑛+1/21𝜑^ℎ𝜼_𝑛+1/21𝜼_𝑛+1/2²₁+1𝜙𝜀Re^{ℎ𝑛+1/221}. (4.29) Lastly, estimate (2.11)(a) and inverse inequality (2.17) give 𝑏^ℎ𝜑^ℎ𝜙^ℎ_𝑛+1/2,𝜼_𝑛+1/2,𝜙^ℎ_𝑛+1/2≤𝐶Reℎ⁻¹𝜑^ℎ𝜙^ℎ_𝑛+1/2²𝜼_𝑛+1/2²₁+1𝜙𝜀Re^{ℎ𝑛+1/221}. (4.30) Apply estimates ||𝜑^ℎ(𝐯)||≤||𝐯|| and |𝜑^ℎ(𝐯)|₁≤𝑑₁|𝐯|₁ from Assumption 3.5 along with estimates (4.26)–(4.30) to get 𝑏^ℎ𝜑^ℎ𝜙^ℎ_𝑛+1/2,𝐮_𝑛+1/2,𝜙^ℎ_𝑛+1/2+𝑏^ℎ𝜑^ℎ𝜼_𝑛+1/2,𝐮_𝑛+1/2,𝜙^ℎ_𝑛+1/2+𝑏^ℎ𝜑^ℎ𝐰^ℎ_𝑛+1/2,𝜼_𝑛+1/2,𝜙^ℎ_𝑛+1/2≤4𝜙𝜀Re^{ℎ𝑛+1/221}+𝐶Re𝐮_𝑛+1/2²₂+ℎ⁻¹𝜼_𝑛+1/2²₁𝜙^ℎ_𝑛+1/2²+𝐶𝑑₁Re‖∇𝐮‖²_{𝐿∞(𝐿2(Ω))}+𝜼_𝑛+1/2²₁𝜼_𝑛+1/2²₁. (4.31) (Time discretization error): First, apply duality estimate on 𝑋×𝑋^′ to get 𝐮_𝑛+1−𝐮_𝑛−𝐮Δ𝑡^𝑡𝑛+1/2,𝜙^ℎ_𝑛+1/2𝐮≤𝐶Re^𝑛+1−𝐮_𝑛−𝐮Δ𝑡^{𝑡𝑛+1/22}₋₁+1𝜙𝜀Re^{ℎ𝑛+1/221}. (4.32) Taylor-expansion about 𝑡_𝑛+1/2 with integral remainder gives 12𝐮_𝑛+1⋅∇𝐮_𝑛+1+1,𝐯2𝐮_𝑛⋅∇𝐮_𝑛=𝐮,𝐯⋅,𝑡_𝑛+1/2⋅∇𝐮⋅,𝑡_𝑛+1/2+1,𝐯2^𝑡𝑛+1_𝑡𝑛+1/2𝑡_𝑛+1𝑑−𝑡²𝑑𝑡²(𝐮()())+1⋅,𝑡⋅∇𝐮⋅,𝑡,𝐯𝑑𝑡2^𝑡𝑛+1/2_𝑡𝑛𝑡−𝑡_𝑛𝑑²𝑑𝑡²(𝐮(⋅,𝑡)⋅∇𝐮(⋅,𝑡),𝐯)𝑑𝑡. (4.33) Add/subtract (𝐮_𝑛+1/2⋅∇𝐮(⋅,𝑡_𝑛+1/2),𝐯) and apply (4.33) to get 𝐮_𝑛+1/2⋅∇𝐮_𝑛+1/2−1,𝐯2𝐮_𝑛+1⋅∇𝐮_𝑛+1−1,𝐯2𝐮_𝑛⋅∇𝐮_𝑛=𝐮,𝐯_𝑛+1/2⋅∇𝐮_𝑛+1/2−𝐮⋅,𝑡_𝑛+1/2,𝐯+𝐮_𝑛+1/2−𝐮⋅,𝑡_𝑛+1/2⋅∇𝐮⋅,𝑡_𝑛+1/2−1,𝐯2^𝑡𝑛+1_𝑡𝑛+1/2𝑡_𝑛+1−𝑡𝐮_𝑡𝑡⋅∇𝐮+𝐮⋅∇𝐮_𝑡𝑡+2𝐮_𝑡⋅∇𝐮_𝑡−1⋅𝐯𝑑𝑡2^𝑡𝑛+1/2_𝑡𝑛𝑡−𝑡_𝑛𝐮_𝑡𝑡⋅∇𝐮+𝐮⋅∇𝐮_𝑡𝑡+2𝐮_𝑡⋅∇𝐮_𝑡⋅𝐯𝑑𝑡. (4.34) Majorize either directly or with (2.7)(a) to get 𝐮_𝑛+1/2⋅∇𝐮_𝑛+1/2−1,𝐯2𝐮_𝑛+1⋅∇𝐮_𝑛+1−1,𝐯2𝐮_𝑛⋅∇𝐮_𝑛‖‖,𝐯≤𝐶∇𝐮_𝑙^∞_(𝐿²_(Ω))𝐮_𝑛+11+𝐮_𝑛1𝐯^ℎ₁, (4.35) or with (2.7)(b), (2.7)(c) and Hölder’s inequality (in time) applied to (4.34) to get 𝐮_𝑛+1/2⋅∇𝐮_𝑛+1/2−1,𝐯2𝐮_𝑛+1⋅∇𝐮_𝑛+1−1,𝐯2𝐮_𝑛⋅∇𝐮_𝑛‖𝐮‖,𝐯≤𝐶_𝑙^∞_(𝐻²_(Ω))𝐮_𝑛+1/2−𝐮⋅,𝑡_𝑛+1/2𝐯^ℎ₁+𝐶Δ𝑡^3/2𝑡^𝑛+1/2‖𝐮‖_𝑙^∞_{(𝑛,𝑛+1;𝐻}²_(Ω))^𝑡𝑛+1_𝑡𝑛𝑡𝐮_𝑡𝑡(⋅,𝑡)²𝑑𝑡^1/2𝐯^ℎ₁+𝐶Δ𝑡^3/2𝑡^𝑛+1/2𝐮_𝑡𝑙^∞_{(𝑛,𝑛+1;𝐿2)}^𝑡𝑛+1_𝑡𝑛𝑡𝐮_𝑡(⋅,𝑡)²₂𝑑𝑡^1/2𝐯^ℎ₁. (4.36) Then, to prove Corollary 4.7, apply (4.32), (4.36) with Young’s inequality give 𝜏_𝑛⁽¹⁾𝐮,𝑝;𝜙^ℎ_𝑛+1/2≤2𝜙𝜀Re^{ℎ𝑛+1/221}𝐮+𝐶Re^𝑛+1−𝐮_𝑛−𝐮Δ𝑡^{𝑡𝑛+1/22}₋₁+𝐶ReΔ𝑡³‖𝐮‖²_𝑙^∞_(𝐻2)𝐮_𝑡𝑡²_{𝐿2(𝑡𝑛,𝑡𝑛+1;𝐿2(Ω))}+𝐶ReΔ𝑡³𝐮_𝑡²_𝑙^∞_(𝐿2(Ω))𝐮_𝑡²_{𝐿2(𝑡𝑛,𝑡𝑛+1;𝐻2(Ω))}. (4.37) We apply (4.35) instead of (4.36) to prove Theorem 4.6: 𝜏_𝑛⁽¹⁾𝐮,𝑝;𝜙^ℎ_𝑛+1/2≤2𝜙𝜀Re^{ℎ𝑛+1/221}𝐮+𝐶Re^𝑛+1−𝐮_𝑛−𝐮Δ𝑡^{𝑡𝑛+1/22}₋₁+𝐶Re‖∇𝐮‖²_𝑙^∞_(𝐿2(Ω))𝐮_𝑛+1²₁+𝐮_𝑛²₁. (4.38) (Deconvolution error): Next, add/subtract 𝜑(𝐮_𝑛+1/2) we write 𝜏_𝑛⁽²⁾𝐮,𝑝;𝐯^ℎ=−𝑏^ℎ𝐮_𝑛+1/2−𝜑𝐮_𝑛+1/2,𝐮_𝑛+1/2,𝜙^ℎ_𝑛+1/2−𝑏^ℎ𝜑𝐮_𝑛+1/2−𝜑^ℎ𝐮_𝑛+1/2,𝐮_𝑛+1/2,𝜙^ℎ_𝑛+1/2+𝜒𝐮_𝑛+1/2−𝜑𝐮_𝑛+1/2,𝜙^ℎ_𝑛+1/2+𝜒𝜑𝐮_𝑛+1/2−𝜑^ℎ𝐮_𝑛+1/2,𝜙^ℎ_𝑛+1/2. (4.39) Then duality on 𝑋×𝑋^′ and Young’s inequalities give 𝜒𝐮_𝑛+1/2−𝜑𝐮_𝑛+1/2,𝜙^ℎ_𝑛+1/2+𝜒𝜑𝐮_𝑛+1/2−𝜑^ℎ𝐮_𝑛+1/2,𝜙^ℎ_𝑛+1/2≤𝐶𝜒²Re𝐮_𝑛+1/2−𝜑𝐮_𝑛+1/2²₋₁+𝐶𝜒²Re𝜑𝐮_𝑛+1/2−𝜑^ℎ𝐮_𝑛+1/2²₋₁+1𝜙𝜀Re^{ℎ𝑛+1/221}, (4.40) and 𝐮∈𝐻²(Ω) along with (2.11)(b) give 𝑏^ℎ𝐮_𝑛+1/2−𝜑𝐮_𝑛+1/2,𝐮_𝑛+1/2,𝜙^ℎ_𝑛+1/2+𝑏^ℎ𝜑𝐮_𝑛+1/2−𝜑^ℎ𝐮_𝑛+1/2,𝐮_𝑛+1/2,𝜙^ℎ_𝑛+1/2≤𝐶Re𝐮_𝑛+1/2²₂𝐮_𝑛+1/2−𝜑𝐮_𝑛+1/2²+𝜑𝐮_𝑛+1/2−𝜑^ℎ𝐮_𝑛+1/2²+1𝜙𝜀Re^{ℎ𝑛+1/221}. (4.41) The estimates (4.37), (4.40), (4.41) with identity (4.39) give 𝜏_𝑛⁽¹⁾𝐮,𝑝;𝜙^ℎ_𝑛+1/2+𝜏_𝑛⁽²⁾𝐮,𝑝;𝜙^ℎ_𝑛+1/2≤4𝜙𝜀Re^{ℎ𝑛+1/221}+𝐶Re𝐸_𝑛. (4.42) Then estimates (3.17) and (3.23) give 𝐸_𝑛∶=𝜒²+‖𝐮‖²_𝑙^∞_(𝐻2(Ω))𝑐_1,𝑗(𝛼,𝛿)²Δ^𝑗+1𝐮_𝑛+1/2²+𝑐_2,𝑗(ℎ,𝛼,𝛿)²𝐮_𝑛+1/2²_𝑘+1+𝐮_𝑛+1−𝐮_𝑛−𝐮Δ𝑡^{𝑡𝑛+1/22}₋₁+‖𝐮‖²_𝑙^∞_(𝐻2)𝐮_𝑛+1/2−𝐮⋅,𝑡_𝑛+1/2²+Δ𝑡³‖𝐮‖²_𝑙^∞_(𝐻2)𝐮_𝑡𝑡²_{𝐿2(𝑡𝑛,𝑡𝑛+1;𝐿2(Ω))}+Δ𝑡³𝐮_𝑡²_𝑙^∞_(𝐿2(Ω))𝐮_𝑡²_{𝐿2(𝑡𝑛,𝑡𝑛+1;𝐻2(Ω))}. (4.43)
Apply estimates from (4.24), (4.25), (4.31), (4.37), (4.42) to (4.23). Set 𝜀=20 and absorb all terms including |𝜙^ℎ_𝑛+1/2|₁ from the right into left-hand-side of (4.23). Sum the resulting inequality on both sides from 𝑛=0 to 𝑛=𝑀−1 to get 𝜙^ℎ_𝑀²+Re⁻¹Δ𝑡^𝑀−1_𝑛=0𝜙^ℎ_𝑛+1/2²₁+𝜒Δ𝑡^𝑀−1_𝑛=0𝜙^ℎ_𝑛+1/2²_⋆ℎ≤𝜙^ℎ₀²+𝐶ReΔ𝑡^𝑀−1_𝑛=0𝐸_𝑛+𝐶ReΔ𝑡^𝑀−1_𝑛=0𝑝_𝑛+1/2−𝑞^ℎ_𝑛+1/2²+𝐶Re𝑑²₁Δ𝑡^𝑀−1_𝑛=0‖∇𝐮‖²_{𝐿∞(𝐿2(Ω))}+𝜼_𝑛+1/2²₁𝜼_𝑛+1/2²₁+𝐶ReΔ𝑡^𝑀−1_𝑛=0𝐮_𝑛+1/2²₂+ℎ⁻¹𝜼_𝑛+1/2²₁𝜙^ℎ_𝑛+1/2²+𝐶𝜒²ReΔ𝑡^𝑀−1_𝑛=0𝜼_𝑛+1/2²₋₁+𝐶Re⁻¹Δ𝑡^𝑀−1_𝑛=0𝜼_𝑛+1/2²₁. (4.44) Estimates (2.23), (2.16)(a) imply sup_𝑛𝜼_𝑛1≤𝐶‖∇𝐮‖_𝑙^∞_(𝐿²_(Ω)),𝜼_𝑛1≤𝐶ℎ𝐮_𝑛2. (4.45) These estimates applied to (4.44) give 𝜙^ℎ_𝑀²+Re⁻¹Δ𝑡^𝑀−1_𝑛=0𝜙^ℎ_𝑛+1/2²₁+𝜒Δ𝑡^𝑀−1_𝑛=0𝜙^ℎ_𝑛+1/2²_⋆ℎ≤𝜙^ℎ₀²+𝐶ReΔ𝑡^𝑀−1_𝑛=0𝐸_𝑛+𝐶ReΔ𝑡^𝑀−1_𝑛=0𝑝_𝑛+1/2−𝑞^ℎ_𝑛+1/2²+𝐶ReΔ𝑡^𝑀−1_𝑛=0𝑑²₁‖∇𝐮‖²_{𝐿∞(𝐿2(Ω))}+Re⁻²𝜼_𝑛+1/2²₁+𝐮_𝑛+1/2²₂𝜙^ℎ_𝑛+1/2²+𝐶𝜒²ReΔ𝑡^𝑀−1_𝑛=0𝜼_𝑛+1/2²₋₁. (4.46) Suppose that the Δ𝑡-restriction (4.16) is satisfied. Then the discrete Gronwall Lemma 2.2 applies to (4.46) and gives 𝜙^ℎ_𝑀²+Re⁻¹Δ𝑡^𝑀−1_𝑛=0𝜙^ℎ_𝑛+1/2²₁+𝜒Δ𝑡^𝑀−1_𝑛=0𝜙^ℎ_𝑛+1/2²_⋆ℎ≤𝐺_𝑀𝜙^ℎ₀²+𝐺_𝑀ReΔ𝑡^𝑀−1_𝑛=0𝐸_𝑛+𝐺_𝑀ReΔ𝑡^𝑀−1_𝑛=0𝑝_𝑛+1/2−𝑞^ℎ_𝑛+1/2²+𝑑²₁‖∇𝐮‖²_{𝐿∞(𝐿2(Ω))}+Re⁻²𝜼_𝑛+1/2²₁+𝐺_𝑀𝜒²ReΔ𝑡^𝑀−1_𝑛=0𝜼_𝑛+1/2²₋₁, (4.47) where 𝐺_𝑀=𝐶expReΔ𝑡^𝑀_𝑛=0𝑔_𝑛𝐮_𝑛²₂,𝑔_𝑛=1−𝐶ReΔ𝑡𝐮_𝑛²₂⁻¹. (4.48) Lastly, the triangle inequality and approximation theory estimates (2.23), (2.16) along with (3.23) applied to (4.47) give 𝐮_𝑀−𝐰^ℎ_𝑀²+Re⁻¹Δ𝑡^𝑀−1_𝑛=0𝐮_𝑛+1/2−𝐰^ℎ_𝑛+1/2²₁≤𝐺_𝑀𝐞^ℎ₀²+𝐶ℎ^2𝑘𝐺_𝑀𝐮_0𝑘+‖𝐮‖²_{𝑙∞(𝐻𝑘(Ω))}+Re⁻¹‖𝐮‖²_{𝑙2(𝐻𝑘+1(Ω))}+𝐺_𝑀ReΔ𝑡^𝑀−1_𝑛=0𝐸_𝑛+𝐺_𝑀Reℎ^2𝑠+2𝑝²_𝑙²_{(𝐻𝑠+1(Ω))}+𝐺_𝑀Re𝑑²₁‖∇𝐮‖²_{𝐿∞(𝐿2(Ω))}+Re⁻²+𝜒²ℎ⁴ℎ^2𝑘‖𝐮‖²_{𝑙2(𝐻𝑘+1(Ω))}. (4.49)
It remains to bound Δ𝑡∑_𝑛𝐸_𝑛.

(Theorem 4.6): Suppose that 𝜕_𝑡𝐮∈𝐿²(𝑋^′)∩𝑙²(𝑋^′). The triangle inequality and (2.18) gives Δ𝑡^𝑀−1_𝑛=𝑛0𝐮_𝑛+1−𝐮_𝑛−𝐮Δ𝑡^{𝑡𝑛+1/22}₋₁≤𝐶𝐮_𝑡²_{𝐿2(𝑋′)}+𝐮_𝑡²_{𝑙2(𝑋′)}. (4.50) Apply (4.38), instead of (4.37), to derive 𝐸_𝑛 in (4.43). Then ReΔ𝑡^𝑀−1_𝑛=0𝐸_𝑛≤𝐸, (4.51) where 𝐸∶=𝐶Re𝜒²+‖𝐮‖²_𝑙^∞_(𝐻2(Ω))𝑐_1,𝑗(𝛼,𝛿)²Δ^𝑗+1𝐮_𝑛+1/2²+𝑐_2,𝑗(ℎ,𝛼,𝛿)²𝐮_𝑛+1/2²_𝑘+1+𝐶Re𝐮_𝑡²_{𝐿2(𝑋′)}+𝐮_𝑡²_{𝑙2(𝑋′)}+‖∇𝐮‖²_𝑙^∞_(𝐿2(Ω))‖∇𝐮‖²_{𝑙2(𝐿2(Ω))}. (4.52) (Corollary 4.7): Suppose that 𝐮∈𝑙^∞(𝐻²(Ω)), 𝐮_𝑡∈𝑙^∞(𝐻¹(Ω))∩𝐿²(𝐻¹(Ω)), 𝐮_𝑡𝑡∈𝐿²(𝐿²(Ω)), and 𝐮_𝑡𝑡𝑡∈𝐿²(𝑋^′). Write 𝐸∶=𝐶Re𝜒²+‖𝐮‖²_𝑙^∞_(𝐻2(Ω))𝑐_1,𝑗(𝛼,𝛿)²Δ^𝑗+1𝐮_𝑛+1/2²+𝑐_2,𝑗(ℎ,𝛼,𝛿)²𝐮_𝑛+1/2²_𝑘+1∶=𝐶Re𝐮_𝑡𝑡𝑡²_𝐿2𝑋′+‖𝐮‖²_𝑙^∞_𝐻2(Ω)𝐮_𝑡𝑡²_{𝐿2𝐿2(Ω)}+⋯+𝐮_𝑡²_𝑙^∞_𝐿2(Ω)𝐮_𝑡²_{𝐿2𝐻2(Ω)}. (4.53) Then apply (2.19), (2.20) to bound 𝐸_𝑛 given in (4.43): ReΔ𝑡^𝑀−1_𝑛=0𝐸_𝑛≤𝐸Δ𝑡⁴. (4.54)

5. Applications

We show that the iterated (modified) Tikhonov regularization operator satisfied Assumption 3.4, 3.5 in Section 5.1 and verify the theoretical convergence rate predicted by Theorem 4.6, Corollary 4.7 in Section 5.2.

5.1. Iterated (Modified) Tikhonov-Lavrentiev Regularization

We will prove that 𝐷_𝛼,𝐽, 𝐷^ℎ_𝛼,𝐽 (Definitions 3.14, 3.15) with the differential filter 𝐺=𝐴⁻¹ satisfies Assumptions 3.4, 3.5. Proposition 3.10 implies that it is enough to show that 𝐷_𝛼,0 satisfies Assumption 3.4. Additionally, we provide sharpened estimates for 𝑑_1,𝑗, 𝑐_1,𝑗, 𝑐_2,𝑗. The key is that 𝐴⁻¹>0 is a continuous function of the Laplace operator −Δ≥0 and hence they commute (on 𝑋). Moreover, 𝐷_𝛼,0>0 is a continuous function of 𝐴⁻¹ so that 𝐷_𝛼,0 commutes with 𝐴⁻¹ and −Δ (on 𝑋).

We first characterize the spectrum of 𝐷_𝛼,0, 𝐷^ℎ_𝛼,0.

Lemma 5.1. Fix 0<𝛼≤1. Define 𝑓∶(0,1]→ℝ and 𝑔∶(0,1]→ℝ by 𝑓(𝑥)1∶=()(𝑥)𝑥1−𝛼𝑥+𝛼,𝑔∶=().1−𝛼𝑥+𝛼 (5.1) The maps 𝑓 and 𝑔 are continuous and 𝑓((0,1])=[1,𝛼⁻¹) and 𝑔((0,1])=(0,1].

Proof. The functions 𝑓, 𝑔 are clearly continuous with 𝑓  decreasing and 𝑔 increasing on (0,1]. Hence, the range of 𝑓 is [1,𝛼⁻¹) and range of 𝑔 is (0,1].

The next result shows that 𝐷_𝛼,0 and 𝐷^ℎ_𝛼,0 satisfy part of Assumptions 3.4, 3.5.

Proposition 5.2. 𝐷_𝛼,0 and 𝐷^ℎ_𝛼,0 (on 𝑋 and 𝑋^ℎ, resp.) are linear, bounded, spd, and commute with 𝐴⁻¹, (𝐴^ℎ)⁻¹ (resp.). Moreover, 𝐷_𝛼,0𝐴⁻¹≤1,𝐷_𝛼,0𝐮₁≤|𝐮|₁𝐷∀𝐮∈𝑋,_ℎ𝛼,0𝐴^ℎ−1≤1,𝐷^ℎ_𝛼,0𝐮^ℎℎ₁≤𝐮^ℎ₁∀𝐮^ℎ∈𝑋^ℎ. (5.2) Hence 𝑑₁=1 in Assumptions 3.4, 3.5.

Proof. It is immediately clear that 𝐷_𝛼,0, 𝐷^ℎ_𝛼,0 are linear. As a consequence, since 𝐴⁻¹>0 with spectrum in (0,1], then 𝐷_𝛼,0=𝑓(𝐴⁻¹) with spectrum contained in [1,𝛼⁻¹) so that 𝐷>0. Therefore, 𝐷_𝛼,0𝐴⁻¹=𝑔(𝐴⁻¹) with spectrum contained in (0,1]. A similar argument shows that (𝐴^ℎ)⁻¹ has spectrum in (0,1], 𝐷^ℎ_𝛼,0 has spectrum in [1,𝛼⁻¹), and 𝐷^ℎ_𝛼,0(𝐴^ℎ)⁻¹ has spectrum in (0,1]. Thus 𝐷_𝛼,0>0, 𝐷^ℎ_𝛼,0>0 and ||𝐷_𝛼,0𝐴⁻¹||≤1 and ||𝐷^ℎ_𝛼,0(𝐴^ℎ)⁻¹||≤1. Therefore, 𝐷_𝛼,0 and 𝐷^ℎ_𝛼,0 are bounded and commute with 𝐴⁻¹ and (𝐴^ℎ)⁻¹, respectively, as discussed above.
The second set of inequalities on each line can be proved with an appropriate choice of 𝐯 and 𝐯^ℎ in Definitions 3.2 and 3.3. Starting with Definition 3.2, take 𝜙=𝐮 and choose 𝐯=Δ𝐷_𝛼,0𝐮. Then integration by parts and the Cauchy-Schwartz inequality give the result. The discrete form is proved using Definition 3.3 and choosing 𝜙=𝐮^ℎ and 𝐯=Δ𝐷^ℎ_𝛼,0𝐮^ℎℎ.

It remains to provide estimates for 𝑐₁ and 𝑐₂, and sharpened estimates for 𝑐_1,𝑗 and 𝑐_2,𝑗. Indeed, as a direct consequence of Propositions 5.3, 5.4, we have, for each 𝑗=0,1,…,𝐽, 𝑐_1,𝑗=𝛼𝛿^2𝑗+1Δ^𝑗+1𝐯,∀𝐯∈𝐻^𝑗+1(Ω),𝑐_2,𝑗=𝐶ℎ+2^𝑗𝛼𝛿^21/2ℎ^𝑘max_{0≤𝑛≤𝑗}𝐷_𝛼,𝑛𝐯_𝑘+1,∀𝐯∈𝐻^𝑘+1(Ω). (5.3)

Proposition 5.3. Let 𝑗=0,1,…,𝐽. Then, for some 0≤𝛽≤𝑗+1, 𝐯−𝐷_𝛼,𝑗𝐯≤𝛼𝛿^2𝛽Δ^𝛽𝐯,∀𝐯∈𝐻^𝛽(Ω). (5.4)

Proof. Using (1.9), we have 𝐷⁻¹_𝛼,0𝐷_𝛼,𝐽𝐯−𝐷_{𝛼,𝐽−1}𝐯=𝐯−𝐴⁻¹𝐷_{𝛼,𝐽−1}𝐯. (5.5) Subtracting (5.5) from the identity 𝐷⁻¹_𝛼,0(𝐯−𝐯)=𝐯−𝐯, (5.6) gives us 𝐷⁻¹_𝛼,0𝐼−𝐷_𝛼,𝐽𝐴⁻¹𝐯−𝐼−𝐷_{𝛼,𝐽−1}𝐴⁻¹𝐯=−𝐴⁻¹𝐼−𝐷_{𝛼,𝐽−1}𝐴⁻¹𝐯. (5.7) Multiplying by 𝐷_𝛼,0, rearranging, simplifying, and using 𝐴−𝐼=−𝛿²Δ (Definition 3.2) gives 𝐼−𝐷_𝛼,𝐽𝐴⁻¹𝐯=−𝐷_𝛼,0𝐴⁻¹𝐼−𝐷_{𝛼,𝐽−1}𝐴⁻¹+𝐼−𝐷_{𝛼,𝐽−1}𝐴⁻¹𝐯=𝐼−𝐷_{𝛼,𝐽−1}𝐴⁻¹𝐼−𝐷_𝛼,0𝐴⁻¹𝐯=𝐼−𝐷_{𝛼,𝐽−1}𝐴⁻¹𝐷_𝛼,0𝐴⁻¹𝐷⁻¹_𝛼,0𝐯=𝐴−𝐼𝐼−𝐷_{𝛼,𝐽−1}𝐴⁻¹𝐷_𝛼,0𝐴⁻¹(1−𝛼)𝐴⁻¹𝐯=+𝛼𝐼𝐴−𝐼𝐼−𝐷_{𝛼,𝐽−1}𝐴⁻¹𝐷_𝛼,0𝐴⁻¹𝛼(𝐴−𝐼)𝐯=−𝛼𝛿²Δ𝐷_𝛼,0𝐴⁻¹𝐼−𝐷_{𝛼,𝐽−1}𝐴⁻¹𝐯. (5.8) Applying recursion, we obtain, for any 0≤𝛽<𝐽,𝐼−𝐷_𝛼,𝐽𝐴⁻¹𝐯=−𝛼𝛿^2𝛽𝐷_{𝛼,𝐽−𝛽}𝐴^−1𝛽Δ^𝛽𝐯. (5.9) Thus, taking norms and applying ||𝐷_{𝛼,𝐽−𝛽}𝐴⁻¹||≤1, we get (5.4).

Proposition 5.4. Let 𝑗=0,1,…,𝐽. Then 𝐷_𝛼,𝑗𝐰−𝐷^ℎ_𝛼,𝑗𝐰^ℎ2≤𝐶ℎ²+2^𝑗+1𝛼𝛿²ℎ^2𝑘max_{0≤𝑛≤𝑗}𝐷_𝛼,𝑛𝐰²_𝑘+1,∀𝐰∈𝐻^𝑘+1(Ω). (5.10)

Proof. Let 𝐯^ℎ_𝑗∈𝑋^ℎ be the 𝐿²-projection of 𝐷_𝛼,𝑗𝐰. Take 𝐯=𝐯^ℎ in (3.34). For 𝑗=1,…,𝐽, let e_𝑗=𝐷_𝛼,𝑗𝐰−𝐷^ℎ_𝛼,𝑗𝐰^ℎ∶=𝜼_𝑗−𝜙^ℎ_𝑗, where 𝜼_𝑗∶=𝐷_𝛼,𝑗𝐰−𝐯^ℎ_𝑗, and 𝜙^ℎ_𝑗∶=𝐷^ℎ_𝛼,𝑗𝐰^ℎ−𝐯^ℎ_𝑗. Subtract (3.34) and (3.35) to get 𝛼𝛿²∇𝜙^ℎ_𝑗,∇𝐯^ℎ+𝜙^ℎ_𝑗,𝐯^ℎ=𝛼𝛿²∇𝜼_𝑗,∇𝐯^ℎ+𝛼𝛿²∇𝐞_𝑗−1,∇𝐯^ℎ. (5.11) Take 𝐯^ℎ=𝜙^ℎ_𝑗 in (5.11) to get 𝛼𝛿²𝜙^ℎ_𝑗²₁+𝜙^ℎ_𝑗²=𝛼𝛿²∇𝜼_𝑗,∇𝜙^ℎ_𝑗+𝛼𝛿²∇𝐞_𝑗−1,∇𝜙^ℎ_𝑗. (5.12) Fix 𝜀>0. Apply Hölder’s and Young’s inequalities to (5.12) to get 𝛼𝛿²𝜙^ℎ_𝑗²₁+𝜙^ℎ_𝑗²≤𝛼𝛿²𝜼_𝑗²₁+𝜀𝛼𝛿²𝐞_𝑗−1²₁+1𝜀𝛼𝛿²𝜙^ℎ_𝑗²₁. (5.13) Taking 𝜀=1 and 𝜀=2 in (5.13) gives 𝜙^ℎ_𝑗²≤𝛼𝛿²𝜼_𝑗²₁+𝛼𝛿²𝐞_𝑗−1²₁.𝛼𝛿²𝜙^ℎ_𝑗²₁+2𝜙^ℎ_𝑗²≤2𝛼𝛿²𝜼_𝑗²₁+4𝛼𝛿²𝐞_𝑗−1²₁. (5.14)(5.15) The triangle inequality and estimate (5.14) give 𝐞_𝑗²≤𝜼_𝑗²+𝛼𝛿²𝜼_𝑗²₁+𝛼𝛿²𝐞_𝑗−1²₁. (5.16) Backward induction, estimate (5.15), and (2.23) give 𝐞_𝑗²≤𝐞₀²₁+ℎ²+𝛼𝛿²1+^𝑗_𝑖=02^𝑖max_{0≤𝑛≤𝑗}inf_{𝐯ℎ∈𝑋ℎ}𝐷_𝛼,𝑛𝐰−𝐯^ℎ2₁. (5.17) It has been shown (Estimate (2.36) in the proof of Lemma 2.7 [25]) that 𝐞₀²₁≤𝐶ℎ²+𝛼𝛿²ℎ^2𝑘𝐷₀𝐰²_𝑘+1. (5.18) Note that ∑^𝑗_𝑖=02^𝑖=2^𝑗+1−1. Then, along with application of (2.16), we prove (5.10).

Corollary 5.5 (convergence estimate). Under the assumptions of Corollary 4.7, suppose further that, for some 𝐽=0,1,…, that 𝐺=𝐴⁻¹, 𝐺^ℎ=(𝐴^ℎ)⁻¹, 𝐷=𝐷_𝛼,𝐽, 𝐷^ℎ=𝐷^ℎ_𝛼,𝐽. If Δ^𝛽𝐮∈𝐿²(Ω) for some 0≤𝛽≤𝐽+1, then 𝐮_𝑀−𝐰^ℎ_𝑀²+Re⁻¹Δ𝑡^𝑀−1_𝑛=0𝐮_𝑛+1/2−𝐰^ℎ_𝑛+1/2²₁≤𝐶1+ℎ²+2^𝐽𝛼𝛿²ℎ^2𝑘+ℎ^2𝑠+2+Δ𝑡⁴+𝛼𝛿^22𝛽. (5.19)

Proof. Apply estimates for 𝑐_1,𝑗, 𝑐_2,𝑗 from (5.3), resulting from Propositions 5.3, 5.4.

5.2. Numerical Testing

This section presents the calculation of a flow with an exact solution to verify the convergence rates of the algorithm. FreeFEM++ [32] was used to run the simulations. The convergence rates are tested against the Taylor-Green vortex problem [13, 33–35]. We use a domain of Ω=(0,1)×(0,1) and take 𝐮=(𝑢₁,𝑢₂), where 𝑢₁𝑥,𝑦,𝑡=−cos(𝑛𝜋𝑥)sin𝑛𝜋𝑦𝑒^{−2𝑛2𝜋2𝑡/𝜏},𝑢₂𝑥,𝑦,𝑡=sin(𝑛𝜋𝑥)cos𝑛𝜋𝑦𝑒^{−2𝑛2𝜋2𝑡/𝜏},𝑝1𝑥,𝑦,𝑡=−4cos(𝑛𝜋𝑥)+cos𝑛𝜋𝑦𝑒^{−2𝑛2𝜋2𝑡/𝜏}. (5.20) The pair (𝐮,𝑝) is a solution the two-dimensional NSE when 𝜏=Re and 𝐟=0.

We used CN discretization in time and P2-P1 elements in space according to Problem 1. That is, we used continuous piecewise quadratic elements for the velocity and continuous piecewise linear elements for the pressure. We chose the spatial discretization elements and parameters 𝑛=1, 𝑇=0.5, 𝜒=0.1 and Re=10,000 as a illustrative example. We chose ℎ=1/𝑚, 𝑑𝑡=(1/4)ℎ, 𝛿=√⁴ℎ and 𝛼=√ℎ, where 𝑚 is the number of mesh divisions per side of [0,1]. These were chosen so that the result of Corollary 5.5 reduces to 𝐮_𝑀−𝐰^ℎ_𝑀+Re⁻¹Δ𝑡^𝑀−1_𝑛=0𝐮_𝑛+1/2−𝐰^ℎ_𝑛+1/2²₁^1/2≤𝐶ℎ²+ℎ^𝐽+1. (5.21)

We summarize the results in Tables 1 and 2. Table 1 displays error estimates corresponding to no iterations; that is, 𝐽=0 in Definition 3.14. For the particular choice of 𝛼 and 𝛿, the computed errors ‖𝑢−𝑤^ℎ‖_∞,0 and ‖∇(𝑢−𝑤^ℎ)‖_2,0 tend to the predicted convergence rate 𝒪(ℎ). Table 2 displays error estimates corresponding to one update; that is, when 𝐽=1 in Definition 3.14. Again, for the particular choice of 𝛼 and 𝛿,the computed errors ‖𝑢−𝑤^ℎ‖_∞,0 and ‖∇(𝑢−𝑤^ℎ)‖_2,0 tend to the predicted convergence rate 𝒪(ℎ²).

Table 1 

Error and convergence rates for Leray-deconvolution with 𝐽=0 for the Taylor-Green vortex with Re=10,000, 𝛼=√ℎ, and 𝛿=√⁴ℎ. Note the convergence rate is approaching 1 as predicted by (5.21).

Table 2 

Error and convergence rates for Leray-deconvolution with 𝐽=1 for the Taylor-Green vortex with Re=10,000, 𝛼=√ℎ, and 𝛿=√⁴ℎ. Note the convergence rate is approaching 2 as predicted by (5.21).

6. Conclusion

It is infeasible to resolve all persistent and energetically significant scales down to the Kolmogorov microscale of 𝒪(Re^−3/4) for turbulent flows in complex domains using direct numerical simulations in a given time constraint. Regularization methods are used to find approximations to the solution. The modification of iterated Tikhonov-Lavrentiev to the modified iterated Tikhonov-Lavrentiev deconvolution in Definition 3.14 is a highly accurate method of solving the deconvolution problem in the Leray-deconvolution model, with errors 𝐮−𝐷_𝛼,0𝐮=𝒪((𝛼𝛿²)^𝐽+1) when applied to the differential filter. We use this result to show that under a regularity assumption, the error between the solutions to the NSE and to the Leray deconvolution model with time relaxation using the modified iterated Tikhonov-Lavrentiev deconvolution and discretized with CN in time and FE’s in space are 𝒪(ℎ^𝑘(ℎ+𝛼𝛿²)+ℎ^𝑠+1+Δ𝑡²+(𝛼𝛿²)^𝐽+1).

We also examined the Taylor-Green vortex problem using Problem 1 with the deconvolution in Definition 3.14. We use this problem because it has an exact analytic solution to the NSE. The regularization parameters 𝛼 and 𝛿 were chosen so that the convergence of the approximate solution to the error would be 𝒪(ℎ^𝐽+1) for 𝐽=0 and 𝐽=1. The convergence rates calculated correspond to those predicted, that is 𝒪(ℎ¹) for 𝐽=0 and 𝒪(ℎ²) for 𝐽=1.