A reduced-order extrapolated Crank–Nicolson finite spectral element method based on POD for the 2D non-stationary Boussinesq equations

Accepted Manuscript A reduced-order extrapolated Crank–Nicolson finite spectral element method based on POD for the 2D non-stationary Boussinesq equations

Zhendong Luo, Fei Teng, Hong Xia

PII: DOI: Reference:

S0022-247X(18)30921-1 https://doi.org/10.1016/j.jmaa.2018.10.092 YJMAA 22680

To appear in:

Journal of Mathematical Analysis and Applications

Received date:

6 August 2018

Please cite this article in press as: Z. Luo et al., A reduced-order extrapolated Crank–Nicolson finite spectral element method based on POD for the 2D non-stationary Boussinesq equations, J. Math. Anal. Appl. (2018), https://doi.org/10.1016/j.jmaa.2018.10.092

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

A reduced-order extrapolated Crank–Nicolson finite spectral element method based on POD for the 2D non-stationary Boussinesq equations$ Zhendong Luoa,∗ , Fei Tengb , Hong Xiab a School b School

of Mathematics and Physics, North China Electric Power University, Beijing 102206, China

of Control and Computer Engineering, North China Electric Power University, Beijing 102206, China

Abstract In this paper, we mainly utilize proper orthogonal decomposition (POD) to reduce the order for the coefficient vector of the classical Crank–Nicolson finite spectral element (CCNFSE) method of the two-dimensional (2D) non-stationary Boussinesq equations about vorticity-stream functions so that the reduced-order method maintains all the advantages of the CCNFSE method. Toward this end, we first establish a CCNFSE format with the second-order time accuracy for the two-dimensional (2D) non-stationary Boussinesq equations about vorticity-stream functions and analyze the existence, stability, and convergence of the CCNFSE solutions. And then, by POD, we establish a reduced-order extrapolated Crank–Nicolson finite spectral element (ROECNFSE) method and analyze the existence, stability, and convergence of the ROECNFSE solutions as well as offer the flowchart for seeking ROECNFSE solutions. Finally, we use two sets of numerical experiments to validate that the numerical computational consequences are accorded with the theoretical ones such that the effectiveness and feasibility of the ROECNFSE method are further verified. Both theory and method in this paper is new and completely different from the existing reduced-order methods.

Keywords: Proper orthogonal decomposition; Crank–Nicolson finite spectral element method; Two-dimensional non-stationary Boussinesq equations; Reduced-order extrapolated Crank–Nicolson finite spectral element method; Existence and stability as well as convergence 2000 MSC: MSC 5N12; 65M15; 65N35

1. Introduction Let Θ ⊂ R2 be a connected and bounded domain. We consider the following two-dimensional (2D) nonlinear and non-stationary Boussinesq equations (see [19, 36]):

$ This

work is jointly supported by the ”Double Tops” Postgraduate Talent Cultivation Project of North China Electric Power University and the National Science Foundation of China (11671106) ∗ corresponding author. Tel.: +86-10-61772167; Fax:+86-10-61772167 Email addresses: [email protected] (Z.D. Luo), [email protected] (F. Teng), [email protected] (H. Xia)

Preprint submitted to Journal of Mathematical Analysis and Applications

October 29, 2018

Problem 1. Find u, v, Q, and p such that  2  ⎧ ∂u ∂ u ∂ 2u ∂u ∂u ∂ p ⎪ ⎪ μ + − = 0, +u +v + ⎪ 2 2 ⎪ ∂t ∂ x ∂ y ∂ x ∂y ∂x ⎪ ⎪  2  ⎪ ⎪ ⎪ ∂v ∂ v ∂ 2v ∂v ∂v ∂ p ⎪ ⎪ ⎪ ⎪ ∂ t − μ ∂ x2 + ∂ y2 + u ∂ x + v ∂ y + ∂ y = Q, ⎪ ⎪ ⎪ ⎪ ⎪ ∂u ∂v ⎪ ⎪ + = 0, ⎪ ⎨ ∂x ∂y  2  ∂Q ∂ Q ∂ 2Q ∂Q ∂Q ⎪ ⎪ +u γ + − +v = 0, 0 ⎪ ∂t 2 2 ⎪ ∂x ∂y ∂x ∂y ⎪ ⎪ ⎪ ⎪ ⎪ u(x, y,t) = ϕu (x, y,t), v(x, y,t) = ϕv (x, y,t), ⎪ ⎪ ⎪ ⎪ ⎪ Q(x, y,t) = Q0 (x, y,t), ⎪ ⎪ ⎪ ⎪ ⎪ u(x, y, 0) = u0 (x, y), v(x, y, 0) = v0 (x, y), ⎪ ⎪ ⎩ Q(x, y, 0) = Q0 (x, y),

(x, y,t)∈Θ × (0, T ), (x, y,t)∈Θ × (0, T ), (x, y,t)∈Θ × (0, T ), (x, y,t)∈Θ × (0, T ),

(1)

(x, y,t)∈∂ Θ × (0, T ], (x, y,t)∈∂ Θ × (0, T ], (x, y)∈Θ, (x, y)∈Θ,

where (u, v)T represents the fluid velocity vector, Q denotes the temperature or heat energy, p is the pressure, √ T is the total time, μ = Pr/Re, Re is the Reynolds number, Pr is the Prandtl number, γ0 = RePr, (ϕu (x, y,t), ϕv (x, y,t))T and Q0 (x, y,t) are, respectively, the given boundary values of the fluid velocity and the temperature, and (u0 (x, y), v0 (x, y))T and Q0 (x, y) are, respectively, the given initial values of the fluid velocity and the temperature. For the sake of discussion, but without losing generality, we will assume that ϕu (x, y,t) = ϕv (x, y,t) = Q0 (x, y,t) = 0 in the following analysis procedure. The 2D nonlinear and non-stationary Boussinesq equations constitute an important mathematical model in fluid dynamics and have been successfully used to simulate the real-world engineering problems (see, e.g., [7, 19, 24, 36]). However, due to their nonlinearity, especially, when their computational domains are some irregular geometrical shapes, we can not usually find out their analytical solutions, so that we have to depend on numerical solutions. At present, finite difference (FD) scheme, finite element (FE) method, finite volume element (FVE) method, and spectral method are regarded as to be four welcome numerical approaches (see [4, 5, 11, 12, 16, 32]). However, the spectral method holds highest accuracy among four numerical methods because it adopts the smooth functions (such as, trigonometric functions, Chebyshev’s polynomials, Jacobi’s polynomials, and Legendre’s polynomials) to approximate unknown function, whereas the FE and FVE methods usually adopt standard polynomials to approximate unknown function and the FD scheme adopts difference quotient to approximate derivative. Especially, the classical finite spectral element method, just as the FE and FVE methods, can suitable for the complicated problems with irregular geometrical shape computational domains so that it is more popular (see [10]). However, so far, there has not been any report about the classical Crank–Nicolson finite spectral element (CCNFSE) method for the 2D non-stationary Boussinesq equations about vorticity-stream functions. Therefore, the first task in this paper is to build the CCNFSE method for the 2D non-stationary Boussinesq equations about vorticity-stream functions and analyze the existence, stability, and convergence of the CCNFSE solutions. Though the CCNFSE method for the 2D non-stationary Boussinesq equations about vorticity-stream functions can achieve higher accuracy than the FD scheme, the FE method, and the FVE method, it also includes great many degrees of freedom (i.e., unknowns), so that it would bring many difficulties to practical applications, e.g., its computing load is very heavy and the round-off errors in the computing process 2

would be accumulated very quickly. Therefore, a crucial issue is how to reduce the freedom degrees in the CCNFSE method so as to ease the computational burden, save CPU consuming time, and alleviate the round-off error accumulation in a way that attains the desired numerical solutions. Proper orthogonal decomposition (POD) is considered as to be an effective and feasible approach of the order reduction for numerical models (see [6, 14, 29, 33]). It can greatly reduce the unknowns in the numerical models. It has been widely applied in signal analysis and pattern recognition (see [9]), statistical computation (see [15]), and computational fluid dynamics (see [31]). In recent years, it has been successfully applied to the order reduction for the Galerkin methods (see, e.g., [17, 18]), the FE methods (see, e.g., [22, 26]), the FD schemes (see, e.g., [28, 34]), the FVE methods (see, e.g., [25, 27]), and the reduced-basis methods (see, e.g., [3, 13, 30]) for partial differential equations. However, the most existing POD-based reduced-order methods (see, e.g., [3, 6, 9, 13, 14, 15, 17, 18, 22, 25, 26, 27, 28, 29, 30, 31, 33, 34]) are built by the POD bases formed from the classical solutions at the all time nodes on [0, T ], before repetitively computing the reduced-order solutions at the same time nodes. These are actually some worthless and repetitive computations. In order to eliminate the repetitive computations, two POD-based reduced-order extrapolated Galerkin spectral methods have been presented (see, e.g., [2, 23]). However, for all we know, at the moment, there has not been any report that the reduced-order extrapolated Crank–Nicolson finite spectral element (ROECNFSE) method for the 2D non-stationary Boussinesq equations about vorticity-stream functions has been established by means of POD to reduce the order for the coefficient vector of the CCNFSE model. Hence, the second task in this paper is to develop the ROECNFSE format based on POD for the 2D non-stationary Boussinesq equations about vorticity-stream functions, which only contains very few unknowns but holds sufficiently high accuracy. Particularly, we only extract the few CCNFSE solutions at the initial several time nodes to form the snapshots. And then, we use them to produce the POD bases and build the ROECNFSE format for finding out the ROECNFSE solutions at all time nodes. Thus, we eliminate the repetitive calculations. This is equivalent to utilizing the existing information (on the quite short time interval [0, T0 ], T0  T ) to forecast the future physical law (on the time interval [T0 , T ]). More importantly, we here adopt some new techniques to analyze the existence, stability, and convergence of the ROECNFSE solutions and we use the error estimates as judgement of the selection of POD bases. Moreover, because the ROECNFSE method holds simultaneously both advantages that the POD technique can reduce the unknowns and the CCNFSE method has the higher accuracy, especially, because the CCNFSE method is completely different from the Galerkin spectral method and the non-stationary Boussinesq equations are more complicated than the parabolic and hyperbolic equations in [2, 23], the ROECNFSE method is completely different from those methods in [2, 23] and the establishment of the ROECNFSE format and the theoretical analysis of the existence, stability, and convergence of the ROECNFSE solutions require more skills and have more difficult than those in [2, 23], but the non-stationary Boussinesq equations have more important applications. Most mainly, the ROECNFSE method is established by means of POD to reduce the order of the coefficient vector of CCNFSE method so that it holds the same bases functions as the CCNFSE method and maintains all the advantages of the CCNFSE method. In addition, we adopt matrix and vector analysis method to discuss the existence, stability, and convergence of the ROECNFSE solutions so that the theoretical analysis becomes more convenient, which is completely different all the existing reduced-order methods. Therefore, the ROECNFSE method is a fire-new development and improvement over the existing reduced-order methods as mentioned before. The rest of this paper is arranged as follows. In Section 2, we build the CCNFSE method for the 2D nonstationary Boussinesq equations about the vorticity-stream functions and analyze the existence, stability, and convergence of the CCNFSE solutions. In Section 3, we first extract two sets of snapshots from the 3

initial few the coefficient vectors of CCNFSE solutions, produce two sets of POD bases from the snapshots, and establish the ROECNFSE formulation in the vector form; then, we analyze the existence, stability, and convergence of the ROECNFSE solutions by the matrix and vector analysis method; finally, we give the flowchart for finding the ROECNFSE solutions. In Section 4, we utilize two sets of numerical experiments to show that the ROECNFSE format is superior to the CCNFSE format and to verify that the numerical computing results are accorded with the theoretical analysis and that the ROECNFSE formulation is efficient for solving the 2D non-stationary Boussinesq equations about the vorticity-stream functions because it can greatly reduce the degrees of freedom, alleviate the calculation load, and save the CPU consuming time and the storage requirements in the computations. Section 5 provides the main conclusions and discussion. 2. The CCNFSE Method for the 2D Non-stationary Boussinesq Equations about the Vorticity-stream Functions 2.1. The Semi-discretized Format about Time When Θ is connected and bounded, with ∂ u/∂ x + ∂ v/∂ y = 0, there exists a unique stream function ψ such that u=

∂ψ ∂ψ , v=− . ∂y ∂x

(2)



Further, there is a vorticity function ω such that ω = ∂ v/∂ x − ∂ u/∂ y = − ∂ ψ 2 /∂ x2 +∂ ψ 2 /∂ y2 . Thus, Problem 1 can be transformed into the following form about the vorticity-stream functions (ψ , ω ) and the temperature or heat energy Q: ⎧ 2 2 ⎨− ∂ ψ − ∂ ψ = ω , (x, y,t) ∈ Θ × (0, T ), (3) ∂ x2 ∂ y2 ⎩ ψ = 0, (x, y,t) ∈ ∂ Θ × [0, T ], ⎧  2  ∂ω ∂ ω ∂ 2ω ∂ψ ∂ω ∂ψ ∂ω ∂Q ⎪ ⎪ μ + − − = , (x, y,t) ∈ Θ × (0, T ), + ⎪ ⎨ ∂t ∂ x2 ∂ y2 ∂y ∂x ∂x ∂y ∂x ⎪ω = 0, (x, y,t) ∈ ∂ Θ × [0, T ], ⎪ ⎪ ⎩ω = ω 0 , (x, y) ∈ Θ,

(4)

⎧  2  ∂Q ∂ Q ∂ 2Q ∂ψ ∂Q ∂ψ ∂Q ⎪ ⎪ ⎪ ⎨ ∂ t − γ0 ∂ x2 + ∂ y2 + ∂ y ∂ x − ∂ x ∂ y = 0, (x, y,t) ∈ Θ × (0, T ), Q = 0, (x, y,t) ∈ ∂ Θ × [0, T ], ⎪ ⎪ ⎪ ⎩Q = Q0 , (x, y) ∈ Θ,

(5)

where ω 0 = ∂ v0 /∂ x − ∂ u0 /∂ y. A step of most key for establishing the CCNFSE formulation is to discretize the first derivatives of time in (4) and (5). Toward this end, let M be a positive integer, Δt = T /M be the time step, and let ω n (x, y), ψ n (x, y), and Qn (x, y) be the approximation of ω (x, y,t), ψ (x, y,t), and Q(x, y,t) at tn = nΔt (n = 0, 1, 2, ..., M), respectively. First, from the first equation in (4), we attain    2 ∂ 2ω ∂ ∂ ω ∂ 2ω ∂ψ ∂ω ∂ψ ∂ω ∂Q = μ + 2 − + + . (6) ∂ t2 ∂t ∂ x2 ∂y ∂y ∂x ∂x ∂y ∂x 4

Thus, by Taylor’s formula and (6), we have

∂ ω n ω n+1 − ω n Δt ∂ 2 ω n + o(Δt 2 ) = − ∂t Δt 2 ∂ t2   ω n+1 − ω n μ ∂ 2 (ω n+1 − ω n ) ∂ 2 (ω n+1 − ω n ) ≈ + − Δt 2 ∂ x2 ∂ y2 1 ∂ ψ n+1 ∂ ω n+1 1 ∂ ψ n ∂ ω n 1 ∂ ψ n+1 ∂ ω n+1 1 ∂ ψ n ∂ ω n 1 ∂ (Qn+1 − Qn ) − − + − . 2 ∂y ∂x 2 ∂y ∂x 2 ∂x ∂y 2 ∂x ∂y 2 ∂x Further, by inputting (7) into the first equation in (4), we have  ω n+1 − ω n μ ∂ 2 (ω n+1 + ω n ) ∂ 2 (ω n+1 + ω n ) + − Δt 2 ∂ x2 ∂ y2 +

+

(7)

1 ∂ ψ n+1 ∂ ω n+1 1 ∂ ψ n ∂ ω n 1 ∂ ψ n+1 ∂ ω n+1 1 ∂ ψ n ∂ ω n + − − 2 ∂y ∂x 2 ∂y ∂x 2 ∂x ∂y 2 ∂x ∂y

1 ∂ (Qn+1 + Qn ) . (8) 2 ∂x Because, when Δt is sufficiently small, ω n+1 = ω n + O(Δt), ψ n−1 = ψ n + O(Δt), and Qn+1 = Qn + O(Δt), from (8) we attain the following semi-discretized formulation with the second-order accuracy in time.  ω n+1 − ω n μ ∂ 2 (ω n+1 + ω n ) ∂ 2 (ω n+1 + ω n ) ∂ ψ n ∂ ω n ∂ ψ n ∂ ω n ∂ Qn + − − = . (9) + Δt 2 ∂ x2 ∂ y2 ∂y ∂x ∂x ∂y ∂x =

By using the same technique as the deriving equation (9), we can attain the following semi-discretized formulation with the second-order accuracy in time.  ∂ ψ n ∂ Qn ∂ ψ n ∂ Qn Qn+1 − Qn γ0 ∂ 2 (Qn+1 + Qn ) ∂ 2 (Qn+1 + Qn ) + − − = 0. (10) + Δt 2 ∂ x2 ∂ y2 ∂y ∂x ∂x ∂y The Sobolev spaces and norms in this paper used are normal (see [1, 19]). Let V = H01 (Θ). Then the variational format for semi-discretized formulation with the second-order accuracy in time is the following. Problem 2. Find (ω n , ψ n , Qn ) ∈ V ×V ×V (n = 1, 2, ..., M) that satisfy 



∂ ψ n−1 ∂ w ∂ ψ n−1 ∂ w (11) + dxdy = ω n−1 wdxdy, ∀w ∈ V, n = 1, 2, ..., M + 1; ∂x ∂x ∂y ∂y Θ Θ   

μ Δt ∂ ω n ∂ w ∂ ω n ∂ w ω nw + + dxdy 2 ∂x ∂x ∂y ∂y Θ   



 μ Δt ∂ ω n−1 ∂ w ∂ ω n−1 ∂ w ∂ ω n−1 ∂ w ∂ ω n−1 ∂ w ω n−1 w + + dxdy − μ Δt + dxdy = 2 ∂x ∂x ∂y ∂y ∂x ∂x ∂y ∂y Θ Θ 

 ∂ ψ n−1 ∂ ω n−1 ∂ ψ n−1 ∂ ω n−1 ∂ Qn−1 + Δt − + wdxdy, ∀w ∈ V, n = 1, 2, ..., M, (12) ∂x ∂y ∂y ∂x ∂x Θ   

γ0 Δt ∂ Qn ∂ w ∂ Qn ∂ w Qn w + + dxdy 2 ∂x ∂x ∂y ∂y Θ    

 γ0 Δt ∂ Qn−1 ∂ w ∂ Qn−1 ∂ w ∂ Qn−1 ∂ w ∂ Qn−1 ∂ w n−1 + dxdy − γ0 Δt + dxdy Q w+ = 2 ∂x ∂x ∂y ∂y ∂x ∂x ∂y ∂y Θ Θ 

 ∂ ψ n−1 ∂ Qn−1 ∂ ψ n−1 ∂ Qn−1 + Δt − wdxdy, ∀w ∈ V, n = 1, 2, ..., M. (13) ∂x ∂y ∂y ∂x Θ 5

Let a1 (ψ , ω , φ ) = e.g., [7, 19]):

 Θ

∂ψ ∂ω ∂ψ ∂ω − ∂x ∂y ∂y ∂x



φ dxdy. Then a1 (ψ , ω , φ ) has the following properties (see,

a1 (ψ , ω , φ ) = a1 (ψ , φ , ω ), a1 (ψ , ω , ω ) = 0, ∀ψ , φ , ω ∈ H01 (Θ);

(14)

|a1 (ψ , ω , φ )|  C ∇(curlψ ) 0 ∇ω 0 ∇φ 0 , ∀ψ , φ , ω ∈

(15)

H01 (Θ).

¯ is given, it is known that the iterative equations (11), (12), and (13) are three relatively When ω 0 ∈ C0 (Θ) independent variational problems of normal elliptic equations so that they have a unique series of solutions (ω n , ψ n , Qn ) ∈ V ×V ×V (n = 1, 2, ..., M) meeting the following stability: ψ n (x, y) 0,∞ + ω n (x, y) 0,∞ + Qn (x, y) 0,∞  σ ,

(16)

where σ and follow-on is a general positive constant, only dependent on ω 0 , Q0 , Re, Pr, and T . Moreover, from the above semi-discretized procedure in time, we can easily attain that the series of solutions (ω n , ψ n , Qn ) ∈ V ×V ×V (n = 1, 2, ..., M) have the second-order accuracy, i.e., ψ (x, y,tn ) − ψ n (x, y) 0,∞ + ω (x, y,tn ) − ω n (x, y) 0,∞ + Q(x, y,tn ) − Qn (x, y) 0,∞  σ Δt 2 . (17) 2.2. The Fully Discretized CCNFSE Format Let ℑN be the quasi-uniform quadrilateral subdivision for Θ (see [19]) and the spectral element subspace be chosen as the following: ¯ : wN |K j ∈ P1 (K j ), K j ∈ ℑN , j = 1, 2, ..., N}, VN = {wN ∈ H01 (Θ) ∩C0 (Θ)

(18)

where N is the number of elements and P1 (K j )’s are formed by the quadrilateral spectral element, i.e.,   P1 (K j ) = span Ni j : 1  i  4 , 4

where Ni j = Nˆ i ◦ Fj−1 (x, y), Nˆ i (ξ , η ) = 14 [1 + cos π (ξ − ξi )][1 + cos π (η − ηi )], (x, y) = Fj (ξ , η ) = ( ∑ Nˆ i (ξ , i=1

4

η )xi j , ∑ Nˆ i (ξ , η )yi j ) is a reversible transformation from K j ∈ ℑN to the referencing quadrilateral Kˆ = i=1

ˆ respectively. [−1, 1] × [−1, 1], and (xi j , yi j ) and (ξi , ηi ) are the vertexes of K j and K, Let RN : H01 (Θ) → VN be the H 1 -orthogonal projection, i.e., for any ϕ ∈ H01 (Θ), there holds

Θ

∇(RN ϕ − ϕ ) · ∇vN dxdy = 0, ∀vN ∈ VN .

Moreover, because ℑN is the quasi-uniform quadrilateral subdivision for Θ, the number of nodes is equal to the number of elements (see [19]), RN has the following important property (see, e.g., [11, Chapter III]). Theorem 3. For any ϕ ∈ H q (Θ) with q  2, we have ∇RN ϕ 0,r  σr ∇ϕ 0,r , ∂ k (RN ϕ − ϕ ) 0  σ N k−q , 0  k  q  N + 1, where σr (r = 2 or ∞, and when r = 2, σr = 1) is also a general positive constant and N is the number of nodes in ℑN . By using the spectral element subspace, we can establish the CCNFSE formulation as follows. 6

Problem 4. Find (ωNn , ψNn , TNn ) ∈ VN ×VN ×VN (n = 1, 2, ..., M) that satisfy  

∂ ψNn−1 ∂ wN ∂ ψNn−1 ∂ wN dxdy = ωNn−1 wN dxdy, ∀wN ∈ VN , n = 1, 2, ..., M + 1; (19) + ∂x ∂x ∂y ∂y Θ Θ  n 

 μ Δt ∂ ωN ∂ wN ∂ ωNn ∂ wN ωNn wN + + dxdy 2 ∂x ∂x ∂y ∂y Θ   

μ Δt ∂ ωNn−1 ∂ wN ∂ ωNn−1 ∂ wN n−1 dxdy = ωN wN + + 2 ∂x ∂x ∂y ∂y Θ  

∂ ψNn−1 ∂ ωNn−1 ∂ ψNn−1 ∂ ωNn−1 ∂ Qn−1 N + Δt wN dxdy − + ∂x ∂y ∂y ∂x ∂x Θ  

∂ ωNn−1 ∂ wN ∂ ωNn−1 ∂ wN dxdy, ∀wN ∈ VN , n = 1, 2, ..., M, − μ Δt (20) + ∂x ∂x ∂y ∂y Θ  

 γ0 Δt ∂ QnN ∂ wN ∂ QnN ∂ wN n + dxdy QN wN + 2 ∂x ∂x ∂y ∂y Θ   

∂ Qn−1 γ0 Δt ∂ Qn−1 n−1 N ∂ wN N ∂ wN dxdy QN w + = + 2 ∂x ∂x ∂y ∂y Θ  

∂ Qn−1 ∂ Qn−1 N ∂ wN N ∂ wN − γ0 Δt dxdy + ∂x ∂x ∂y ∂y Θ  

∂ ψNn−1 ∂ Qn−1 ∂ ψNn−1 ∂ Qn−1 N N + Δt wN dxdy, ∀wN ∈ VN , n = 1, 2, ..., M, (21) − ∂x ∂y ∂y ∂x Θ where ωN0 = RN ω 0 and Q0N = RN Q0 . Because the left hand sides of (19), (20), and (21) are, respectively, three bounded and positive definite bilinear functionals in VN and, for the known ωNn−1 , ψNn−1 , and Qn−1 , the right hand sides of (19), (20), and (21) are, respectively, three bounded linear functionals about wN , by the Lax–Milgram theorem (see, e.g., [11, 19]), we conclude that Problem 2 exists a unique series of solutions (ωNn , ψNn , QnN ) ∈ VN ×VN ×VN (n = 1, 2, ..., M) satisfying the following stability: ψNn 0,∞ + ωNn 0,∞ + QnN 0,∞  σ .

(22)

by using the Moreover, when (ψ , ω , Q) ∈ finite element and CCNFSE methods (see, e.g., [11, 19])), we can attain the following error estimates [W q,∞ (Θ) ∩ H01 (Θ)] × [W q,∞ (Θ) ∩ H01 (Θ)] × [W q,∞ (Θ) ∩ H01 (Θ)],

ω (x, y,tn ) − ωNn (x, y) 0,∞  σ (Δt 2 + N −q ), 2  q  N + 1, n = 1, 2, ..., M,

(23)

 σ (Δt 2 + N −q−1 ), 2  q  N + 1, n = 1, 2, ..., M,

(24)

ψ (x, y,tn ) − ψNn (x, y) 0,∞ + N −1 ψ (x, y,tn ) − ψNn (x, y) 1,∞

Q(x, y,tn ) − QnN (x, y) 0,∞ + N −1 Q(x, y,tn ) − QnN (x, y) 1,∞  σ (Δt 2 + N −q−1 ), 2  q  N + 1, n = 1, 2, ..., M. Set n n n n n n n n n n n n T , ω21 , ω31 , ω41 , ω12 , ω22 , ω32 , ω42 , ..., ω1N , ω2N , ω3N , ω4N ) , ωn = (ω11

7

(25)

n n n n n n n n n n n n T ψn = (ψ11 , ψ21 , ψ31 , ψ41 , ψ12 , ψ22 , ψ32 , ψ42 , ..., ψ1N , ψ2N , ψ3N , ψ4N ) ,

Qn = (Qn11 , Qn21 , Dn31 , Qn41 , Qn12 , Qn22 , Qn32 , Qn42 , ..., Qn1N , Qn2N , Qn3N , Qn4N )T . N˜ = (N11 , N21 , N31 , N41 , N12 , N22 , N32 , N42 , ..., N1N , N2N , N3N , N4N )T . Then we have

ωNn =

N

4

∑ ∑ ωi j Ni j = N˜ · ωn ,

j=1 i=1

ψNn =

N

4

∑ ∑ ψi j Ni j = N˜ · ψn ,

QnN =

j=1 i=1

N

4

∑ ∑ Qi j Ni j = N˜ · Qn .

j=1 i=1

Thus, Problem 4 can rewritten as the following matrix form. Problem 5. Find (ωn , ψn , Qn ) ∈ R4N × R4N × R4N (n = 1, 2, ..., M) that satisfy ˜ n−1 = Cωn−1 , n = 1, 2, ..., M + 1; Aψ

(26)

ˆ n−1 , n = 1, 2, ..., M, Aωn = Aωn−1 + Δt B(ψNn−1 )ωn−1 + Δt CQ ˆ n−1 + Δt B( ˆ ψ n−1 )Qn−1 , n = 1, 2, ..., M, ˆ n = AQ AQ N

(27) (28)

where A˜ = diag{ A˜ 11 , A˜ 22 , ..., A˜ NN }, C = diag{C11 , C22 , ..., CNN }, A = diag {A11 , A22 , ..., ANN }, Aˆ = diag ˆ ψ n−1 ) = diag { Aˆ 11 , Aˆ 22 , ..., Aˆ NN }, Cˆ = diag{Cˆ 11 , Cˆ 22 , ..., Cˆ NN }, B(ψNn−1 ) = diag{B11 , B22 , ..., BNN }, B( N ˆ ˆ ˆ { B11 , B22 , ..., BNN }, and     ∂ NiI CIJ = NiI N jJ dxdy , Cˆ IJ = , N jJ dxdy Θ Θ ∂x 4×4 4×4    ∂ NiI ∂ N jJ ∂ NiI ∂ N jJ , dxdy A˜ IJ = + ∂x ∂x ∂y ∂y Θ 4×4

μ Δt ˜ γ0 Δt ˜ AIJ , Aˆ IJ = CIJ + AIJ , AIJ = CIJ + 2 2     ∂ ψNn−1 ∂ NiI ∂ ψNn−1 ∂ NiI N jJ dxdy − μ A˜ IJ , − BIJ = ∂x ∂y ∂y ∂x Θ 4×4     n−1 n−1 ∂ ψN ∂ NiI ∂ ψN ∂ NiI N jJ dxdy Bˆ IJ = − γ0 A˜ IJ . − ∂x ∂y ∂y ∂x Θ 4×4

Remark 6. As long as of nodes N, the Reynolds number Re and the Prandtl  the time step Δt,√the number 0 number Pr, i.e., μ = Pr/Re and γ0 = RePr, Q (x, y), and ω 0 = ∂ v0 /∂ x − ∂ u0 /∂ y are given, by solving Problem 5, we can attain three series of the coefficient vectors ψn , ωn , and Qn (n = 1, 2, ..., M). Finally, we can attain the CCNFSE solutions QnN and (unN , vnN ) (n = 1, 2, ..., M) for the heat energy and the velocity by unN = ∂ ψNn /∂ y and vnN = −∂ ψNn /∂ x. However, the CCNFSE format, i.e., Problem 4 or Problem 5, includes many unknowns, so that we have to lessen the unknowns of Problem 5 by POD.

8

3. The ROECNFSE Method Based on POD for the 2D Non-stationary Boussinesq Equations 3.1. Production of POD Bases First, we extract, respectively, the initial L coefficient vectors ωn and Qn (n = 1, 2, ..., L) from the series of the coefficient vectors ωn and Qn (n = 1, 2, ..., M) and use them to form two snapshot matrices P1 = (ω1 , ω2 , ..., ωL ) and P2 = (Q1 , Q2 , ..., QL ) with volume 4N ×L. Let λi, j > 0 (i = 1, 2; j = 1, 2, · · · , ri = rank(Pi )) be the positive eigenvalues of Pi PiT arranged degressively and let Ui = (φi1 , φi2 , · · · , φiri ) ∈ R4N×ri be the homologous orthonormal eigenvectors of Pi PiT . Then, two sets of POD basis Φi = (φi1 , φi2 , · · · , φidi ) (di  ri ; i = 1, 2) are attained from the initial di vectors in Ui that satisfies the following formulas (see [28]):  Pi − Φi ΦTi Pi 2,2 = λi,di +1 , i = 1, 2, (29) where Pi 2,2 = supχ0 Pi χ 2 / χ 2 and χ 2 is the norm of vector χ. Furthermore, we have ωn − Φ1 ΦT1 ωn 2 = (P1 − Φ1 ΦT1 P1 )en 2   P1 − Φ1 ΦT1 P1 2,2 en 2  λ1,d1 +1 , n = 1, 2, ..., L,

(30)

Qn − Φ2 ΦT2 Qn 2 = (P2 − Φ2 ΦT2 P2 )en 2   P2 − Φ2 ΦT2 P2 2,2 en 2  λ2,d2 +1 , n = 1, 2, ..., L,

(31)

where en (n = 1, 2, · · · , L) are the unit vectors with n-th component being 1. Therefore, Φi = (φi1 , φi2 , · · · , φid ) (i = 1, 2) are two sets of optimal POD basis. Remark 7. Because the order 4N of the matrix Pi Pi T is greatly larger than the order L of the matrix Pi T Pi , i.e., the number of the nodes 4N in ℑN is greatly larger than that of extracted snapshots L, nevertheless both positive eigenvalues λi, j ( j = 1, 2, ..., ri ; i = 1, 2) are identical, we may first seek out the eigenvalues λi, j ( j = 1, 2, ..., ri ; i = 1, 2) of Pi T Pi and the relative eigenvectors ϕij ( j = 1, 2, ..., ri ; i = 1, 2), then, with the formula  φij = Pi ϕij / λi, j ( j = 1, 2, ..., ri ; i = 1, 2), we can attain the eigenvectors φij ( j = 1, 2, ..., ri ; i = 1, 2) corresponding to the positive eigenvalues λi, j ( j = 1, 2, ..., ri ; i = 1, 2) of Pi Pi T . Thus, we can expediently attain the POD basis. 3.2. Establishment of the ROECNFSE Format From (30) and (31) in Section 3.1, we have, respectively, attained the first L (L  N) coefficient vectors of ROECNFSE solutions ωnd = Φ1 ΦT1 ωn =: Φ1 βnω and Qnd = Φ2 ΦT2 Qn =: Φ2 βnQ (n = 1, 2, . . . , L), where n , ω n , ω n , ω n , ω n , ω n , ω n , ω n , ..., ω n , ω n , ω n , ω n )T , Qn = (Qn , Qn , Qn , ωnd = (ωd11 d21 d31 d41 d12 d22 d32 d42 d1N d2N d3N d4N d d11 d21 d31 n n n , Qd41 , Qd12 , Qnd22 , Qnd32 , Qnd42 , ..., Qnd1N , Qnd2N , Qnd3N , Qnd4N )T , βnω = (βωn 1 , βωn 2 , . . . , βωn d1 )T , and βnQ = (βQ1 n , . . . , β n )T . As long as we replace the coefficient vectors ωn and Qn in Problem 5 with ωn = Φ βn βQ2 1 ω Qd2 d and Qnd = Φ2 βnQ (n = L + 1, L + 2, . . . , N), respectively, we can obtain the following ROECNFSE format. ⎧ ⎪ Φ1 βnω = Φ1 ΦT1 ωn , Φ2 βnQ = Φ2 ΦT2 Qn , ⎪ ⎪ ⎪ ⎪ ˜ 1  n  L; ⎪ ψnd = A˜ −1 CΦ1 βnω , ψdn = ψnd · N, ⎪ ⎪ ⎪ ⎨Φ βn = Φ βn−1 + Δt A−1 B(ψ n−1 )Φ βn−1 + Δt A−1 CΦ ˆ 2 βn−1 , L + 1  n  M; 1 ω 1 ω 1 ω Q d (32) n n−1 ˆ ψ n−1 )Φ2 βn−1 , L + 1  n  M; ˆ −1 B( ⎪ β = Φ β + Δt A Φ ⎪ 2 2 Q Q Q d ⎪ ⎪ ⎪ ˜ n = L + 1, L + 2, ..., M; ⎪ ψnd = A˜ −1 CΦ1 βnω , ψdn = ψnd · N, ⎪ ⎪ ⎪ ⎩ωn = Φ βn , Qn = Φ βn , n = 1, 2, ..., M, 1 ω 2 Q d d 9

where ωn and Qn (n = 1, 2, . . . , L) are the initial L coefficient vectors in Problem 5, respectively, and ψnd = n , ψ n , ψ n , ψ n , ψ n , ψ n , ψ n , ψ n , ..., ψ n , ψ n , ψ n , ψ n )T , and the matrices A, ˜ A, A, ˆ B, (ψd11 d21 d31 d41 d12 d22 d32 d42 d1N d2N d3N d4N ˆ ˆ B, C, and C are provided in Problem 5. Farther, the format (32) is simplified into ⎧ ⎪ βnω = ΦT1 ωn , βnQ = ΦT2 Qn , ⎪ ⎪ ⎪ ⎪ ˜ 1  n  L; ⎪ ψnd = A˜ −1 CΦ1 βnω , ψdn = ψnd · N, ⎪ ⎪ ⎪ ⎪ ⎪ T −1 ˆ n−1 ⎪ ⎨βnω = βωn−1 + ΔtΦT1 A−1 B(ψdn−1 )Φ1 βn−1 ω + ΔtΦ1 A CΦ2 βQ , L + 1  n  M; (33) ˆ ψ n−1 )Φ2 βn−1 , L + 1  n  M; ⎪ βnQ = βQn−1 + ΔtΦT2 Aˆ −1 B( ⎪ Q d ⎪ ⎪ ⎪ ⎪ n −1 n n n ˜ ⎪ ⎪ n = L + 1, L + 2, ..., M; ⎪ψd = A˜ CΦ1 βω , ψd = ψd · N, ⎪ ⎪ ⎪ ⎩ n ωd = Φ1 βnω , Qnd = Φ2 βnQ , n = 1, 2, ..., M. Remark 8. After we have obtained ψnd and Qnd (n = 1, 2, ..., M) from (33), we can attain the ROECNFSE ˜ ∂ y and vn = −ψn · ∂ N/ ˜ ∂x solutions (und , vnd ) and Qnd of the velocity and the heat energy by und = ψnd · ∂ N/ d d ˜ respectively. Because the CCNFSE format, Problem 5, has 2 × 4N unknowns at as well as Qnd = Qnd · N, each time level, whereas the ROECNFSE format (33) at the same time level only has d1 + d2 unknowns (di  4N; i = 1, 2), the ROECNFSE format, (32) or (33), is clearly superior to the CCNFSE format, Problem 5. 3.3. Existence, Stability, and Convergence of the ROECNFSE Solutions To analyze the existence, stability, and convergence of the ROECNFSE solutions, we consider the maxnorms of matrix and vector (for the more details, see [38]), which are, respectively, defined dy D ∞ = max

l

∑ |di j |,

1im j=1

∀ D = (di j )m×l ∈ Rm × Rl ,

χ ∞ = max |χ j |, ∀χ = (χ1 , χ2 , ..., χm )T ∈ Rm . 1 jm

For the ROECNFSE solutions, we have the following result on the existence, stability, and convergence. Theorem 9. If ω 0 = ∂ v0 /∂ x − ∂ u0 /∂ y ∈ W 0,∞ (Θ), equivalently, (u0 , v0 ) ∈ W 1,∞ (Θ) ×W 1,∞ (Θ), and Q0 ∈ W 1,∞ (Θ), then the ROECNFSE solutions (ωdn , ψdn , Qnd ) are existing and unique and meet the following stability:

(34) ωdn 0,∞ + ψdn 0,∞ + Qnd 0,∞  σ ω 0 ∞ + Q0 ∞ , n = 1, 2, ..., M. Farther, when (ω , ψ , Q) ∈ [W q,∞ (Θ) ∩ H01 (Θ)] × [W q,∞ (Θ) ∩ H01 (Θ)] × [W q,∞ (Θ) ∩ H01 (Θ)] (2  q  N + 1), we have the following error estimates ω (x, y,tn ) − ωdn 0,∞ + Q(x, y,tn ) − Qnd 0,∞     2 −q  σ Δt + N + λ1,d1 +1 + λ2,d2 +1 , n = 1, 2, ..., M;

(35)

ψ (x, y,tn ) − ψdn 0,∞ + N −1 ψ (x, y,tn ) − ψdn 1,∞      σ Δt 2 + N −q−1 + N −1 λ1,d1 +1 + λ2,d2 +1 , n = 1, 2, ..., M.

(36)

10

Proof. First, we can easily know from the format (33) and Remark 8 that the ROECNFSE solutions (ωdn , ψdn , Qnd ) or (und , vnd , Qnd ) (n = 1, 2, ..., M) are existing and unique. Next, we analyze the stability of the ROECNFSE solutions. With ωnd = Φ1 βnω and Qnd = Φ2 βnQ (n = 1, 2, ..., M), the system of equations (32) can be reverted into the following format: ⎧ ⎪ ωnd = Φ1 ωn , Qnd = Φ2 Qn , ⎪ ⎪ ⎪ ⎪ ⎪ ˜ 1  n  L; ⎪ ψn = A˜ −1 Cωnd , ψdn = ψnd · N, ⎪ ⎨ d ˆ n−1 , L + 1  n  M; (37) ωnd = ωn−1 + Δt A−1 B(ψdn−1 )ωn−1 + Δt A−1 CQ d d d ⎪ ⎪ ⎪ n n−1 −1 n−1 n−1 ⎪ ˆ ψ )Q , L + 1  n  M; Qd = Qd + Δt Aˆ B( ⎪ d d ⎪ ⎪ ⎪ ⎩ψn = A˜ −1 Cωn , L + 1  n  M. d d Moreover, from the FE method (see, e.g., [19, Lemma 1.22 and Theorems 1.31 and 1.33]) and the spectral element method (see, e.g., [11, Chapter III]), we can attain the following inequalities ˆ ∞  σ N; A˜ −1 ∞ + A−1 ∞ + Aˆ −1 ∞  σ N −1 ; B ∞ + B −1 ˆ C ∞  σ , C ∞  σ ; C ∞  σ N.

(38)

Thus, by (37) and (38), we obtain ψnd ∞  σ N −1 ωnd ∞ , n = 1, 2, ..., M;

(39)

ωnd ∞  σ ω0 ∞ , n = 1, 2, ..., L,

(40)

 σ Q ∞ , n = 1, 2, ..., L,

(41)

n−1 n−1  ωn−1 ∞ , n = L + 1, L + 2, ..., M, d ∞ + σ Δt ωd ∞ + σ Δt Q n−1 −1 n−1  Qn−1 ∞ , n = L + 1, L + 2, ..., M. d ∞ + σ Δt ωd ∞ + σ ΔtN Q

(42)

Qnd ∞ ωnd ∞ Qnd ∞

0

(43)

Summing (42)+(43) from L + 1 to n and using Theorem 3, (40), and (41), we attain ωnd ∞ + Qnd ∞  σ ( ω0 ∞ + Q0 ∞ ) + σ Δt

n−1

∑ ( ωid ∞ + Qi ∞ ),

i=L

n = L + 1, L + 2, ..., M.

(44)

Using the Gronwall inequality (see [19, 38]) for (44), we gain

ωnd ∞ + Qnd ∞  σ ω0 ∞ + Q0 ∞ exp[Δt(n − L)]  σ ( ω0 ∞ + Q0 ∞ ), n = L + 1, L + 2, ..., M. Uniting (45) with (39), (40), and (41), we attain

ωnd ∞ + ψnd ∞ + Qnd ∞  σ ω0 ∞ + Q0 ∞ , n = 1, 2, ..., M.

(45)

(46)

˜ ∞  1, from (46), we immediately attain (34). Because ωdn = N˜ · ωnd , ψdn = N˜ · ψnd , Qnd = N˜ · Qnd , and N Finally, we analyze the convergence of the ROECNFSE solutions. Because ωNn = N˜ · ωn , ψNn = N˜ · ψn , n ˜ ∞  1, and χ ∞  χ 2 (∀χ ∈ R4N ), when n = 1, 2, ..., L, with (30), (31), (37), and (38), QN = N˜ · Qn , N we have  ˜ ∞ ωn − ωnd ∞  ωn − Φ1 ΦT1 ωn 2  λ1,d1 +1 , n = 1, 2, ..., L, ωNn − ωdn 0,∞  N (47)  ˜ ∞ Qn − Qnd ∞  Qn − Φ2 ΦT2 Qn 2  λ2,d2 +1 , n = 1, 2, ..., L, QnN − Qnd 0,∞  N (48) ˜ ∞ ψn − ψnd ∞  ( A˜ −1 C)(ωn − ωnd ) ∞ ψNn − ψdn 0,∞  N 11

 A˜ −1 ∞ C ∞ ωn − ωnd ∞  σ N −1

 λ1,d+1 , n = 1, 2, ..., L.

Furthermore, by the inverse estimate theorem (see, e.g., [11, Chapter III]), we have  ψNn − ψdn 1,∞  σ N ψNn − ψdn ∞  σ λd+1 , n = 1, 2, ..., L.

(49)

(50)

From (26) and the fourth equation in (37), with (38), we have ˜ ∞ ψn − ψnd ∞  ( A˜ −1 C)(ωn − ωnd ) ∞ ψNn − ψdn 0,∞  N  A˜ −1 ∞ C(ωn − ωnd ) ∞  σ N −1 ωNn − ωdn 0,∞ , n = L + 1, L + 2, ..., M.

(51)

Farther, with the inverse estimate theorem (see, e.g., [11, Chapter III]), we attain ψNn − ψdn 1,∞  σ ωNn − ωdn 0,∞ , n = L + 1, L + 2, ..., M.

(52)

From (27) and the second equation in (37), with (34), (38), (52), we have ωn − ωnd ∞ −1 ˆ n−1  ωn−1 − ωn−1 + Δt A−1 [B(ψNn−1 )ωn−1 − B(ψdn−1 )ωn−1 − Qn−1 d d ] ∞ + Δt A C(Q d ) ∞ −1 n−1 n−1 n−1  ωn−1 − ωn−1 + B(ψdn−1 )C−1 C(ωn−1 − ωn−1 d ∞ + Δt A [B(ψN − ψd )ωd d )] ∞

ˆ n−1 − Qn−1 ) ∞ + Δt A−1 C(Q d n−1 n−1 n−1 n−1 n−1  ωn−1 − ωn−1 − Qn−1 d ∞ + σ Δt( ψN − ψd 1,∞ + ωN − ωd 0,∞ + Q d ∞ ) n−1 n−1 n−1  ωn−1 − ωn−1 − Qn−1 d ∞ + σ Δt( ωN − ωd 0,∞ + Q d ∞ ),

n = L + 1, L + 2, ..., M.

(53)

From (28) and the third equation in (37), with (34), (38), (52), we have ˆ ψ n−1 )Qn−1 − B( ˆ ψ n−1 )Qn−1 ] ∞ Qn − Qnd ∞  Qn−1 − Qn−1 + Δt Aˆ −1 [ B( N d d d n−1 n−1 ˆ −1 ˆ n−1 ˆ ψ n−1 )(Qn−1 − Qn−1 )] ∞  Qn−1 − Qn−1 + B( d ∞ + Δt A [ B(ψN − ψd )Q d d n−1 n−1 n−1  Qn−1 − Qn−1 − Qn−1 d ∞ + σ Δt( ψN − ψd 1,∞ + Q d ∞ ) n−1 n−1 n−1 n−1  Qn−1 − Qn−1 d ∞ + σ Δt( ωN − ωd 0,∞ + QN − Qd ∞ ),

n = L + 1, L + 2, ..., M.

(54)

Summing (53)+(54) from L + 1 to n and using (47) and (48), we attain   ωn − ωnd ∞ + Qn − Qnd ∞  λ1,d1 +1 + λ2,d2 +1 +σ Δt

n−1

∑ ( ωNi − ωdi 0,∞ + QiN − Qid ∞ ), n = L + 1, L + 2, ..., M.

(55)

i=L

Further, we obtain ˜ n − ωnd ) 0,∞ + N(Q ˜ n − Qnd ) 0,∞ ωNn − ωdn 0,∞ + QnN − Qnd 0,∞ = N(ω  ωn − ωnd ∞ + Qn − Qnd ∞   n−1  λ1,d1 +1 + λ2,d2 +1 + σ Δt ∑ ( ωNi − ωdi 0,∞ + QiN − Qid ∞ ), n = L + 1, L + 2, ..., M. (56) i=L

12

Applying the Gronwall inequality (see [19, 38]) to (56), we attain    n n n n λ1,d1 +1 + λ2,d2 +1 exp[σ Δt(n − L)] ωN − ωd 0,∞ + QN − Qd 0,∞  σ



λ1,d1 +1 +





λ2,d2 +1 , n = L + 1, L + 2, ..., M.

(57)

Uniting (57) with (51) and (52), we attain ψNn − ψdn 0,∞ + N −1 ψNn − ψdn 1,∞     σ N −1 λ1,d1 +1 + λ2,d2 +1 , n = L + 1, L + 2, ..., M.

(58)

Uniting (23) and (24) with (57) and (58), respectively, we attain (35) and (36). This finishes the proof of Theorem 9.  Because u = ∂ ψ /∂ y and v = −∂ ψ /∂ x, unN = ∂ ψNn /∂ y and vnN = −∂ ψNn /∂ x, and und = ∂ ψdn /∂ y and vnd = −∂ ψdn /∂ x, we immediately attain the following result. Theorem 10. Under the same conditions as Theorem 9, the 2D non-stationary Boussinesq equations, Problem 1, has a unique set of the ROECNFSE solutions (und , vnd , Qnd ) that meet the following stability: und 0,∞ + vnd 0,∞ + Qnd 0,∞  σ ( u0 0,∞ + v0 0,∞ + Q0 0,∞ ), n = 1, 2, ..., M,

(59)

and the following convergence: u(x, y,tn ) − und 0,∞ + v(x, y,tn ) − vnd 0,∞ + Q(x, y,tn ) − Qnd 0,∞     2 −q  σ Δt + N + λ1,d1 +1 + λ2,d2 +1 , n = 1, 2, ..., M,

(60)

where 2  q  N + 1.   Remark 11. The error factor λ1,d1 +1 + λ2,d2 +1 in the error estimates in Theorem 10 is caused for the order reduction of the CCNFSE format, which could be served as the suggestions of the POD basis choice. The more concrete explanations are as follows. (1) In order to ensure solutions with the desired accuracy, it is necessary to choose d1  the ROECNFSE  and d2 that satisfies λ1,d1 +1 = λ2,d2 +1  max{Δt 2 , N −q }. (2) In order to ensure that the error estimates in Theorems 9 and 10 attain an optimal order in the numerical computations, we usually choose Δt and N satisfying Δt 2 = O(N −q ). In addition, a lot of numerical experiments have shown that the eigenvalues λi, j ( j  = 1, 2, ..., L;  i = 1, 2) of the symmetrical matrix T Pi Pi are usually decreasing quickly to be close to zero, so λ1,d1 +1 + λ2,d2 +1  max{Δt 2 , N −q } can be satisfied when di  6 (i = 1, 2). 3.4. The Flowchart for Finding the ROECNFSE Solutions Now, we provide the flowchart of finding the ROECNFSE solutions for the 2D non-stationary Boussinesq equations, which consists of the following five steps. Step 1. For given the number Re and the Prandtl  time step Δt, the √ number 0of nodes N, 0the Reynolds 0 number Pr, i.e., μ = Pr/Re and γ0 = RePr, Q (x, y), and ω = ∂ v /∂ x − ∂ u0 /∂ y, solving the classical 13

n n Problem coefficient 15 at2 the initial the

L steps obtains

vectors (ω , Q ) (1  n  L) and form the matrices L 1 2 L P1 = ω , ω , ..., ω 4N×L and P2 = Q , Q , ..., Q 4N×L .

Step 2. Find the positive eigenvalues λi,1  λi,2  . . .  λi,ri > 0 (ri = rank(Pi )) and the corresponding eigenvectors ϕij ( j = 1, 2, . . . , r) of the matrix PiT Pi (i = 1, 2).  2 −q Step 3. Determine the numbers di of POD basis by λi,di +1  max{Δt  , N } (2  q  N + 1) and produce i i i i i the POD basis Φi = (φ1 , φ2 , ..., φdi ) with the formulas φ j = Pi ϕ j / λi, j (1  j  di , i = 1, 2). Step 4. Attain (ψnd , Qnd ) from the coefficient vector equation (33), which only includes d1 + d2 degrees of freedom at each time step, and then, achieve the ROECNFSE solutions (und , vnd )T and Qnd of the velocity and ˜ ∂ y, vn = −ψn · ∂ N/ ˜ ∂ x, and Qn = Qn · N˜ (n = 1, 2, ..., M). heat energy by the formulas und = ψnd · ∂ N/ d d d d n+1 n+1 n n 2 −q Step 5. If und − un+1 d 0,∞ + vd − vd 0,∞ + Qd − Qd 0,∞  max{Δt , N } (n = L, L + 1, . . . , M − 1), n−L−i n−L−i i i and Q = Qd (i = 1, 2, . . . , L) and return to Step 2. end. Else, let ω = ωd

4. Some Numerical Experiments In this section, we employ two set of numerical experiments to explain the superiorities of the ROECNFSE format for the 2D non-stationary Boussinesq equations about the vorticity-stream functions. 4.1. The numerical experiments of square cavity flow In this experiment, we take the computational domain as Θ = (0, 1) × (0, 1), Re = 103 , Pr = 0.7, the side length Δx = Δy = 0.01 of quadrilateral elements in ℑN , i.e., N = 104 , the temporal step Δt = 0.0001, the initial values and boundary values of the velocity vector and heat energy all as zero except for Q0 (1, 1) = Q0 (1, 1,t) = 1. Thus, from (23), (24), and (25), we can deduce that the theoretical errors for the CCNFSE solutions are O(10−8 ) when q = 2. We first compute out 20 CCNFSE solutions (ωn , Qn ) (n = 1, 2, ..., 20) at the initial 20 steps by the CCN

FSE format and use them to form two snapshot matrices P1 = ω1 , ω2 , ..., ω20 and P2 = Q1 , Q2 , ..., Q20 . And then, we compute out two sets of eigenvalues λi, j (their orders are degressively arranged) and eigenT vectors ϕi, j ( j = 1, 2, ..., 20; i = 1, 2) of the matrices P i Pi , further attain the POD bases according on Steps 2 and 3 in Section 3.4. By reckoning, we achieve that λi,7  2 × 10−8 (i = 1, 2). Therefore, we only need to choose the initial six eigenvectors ϕij ( j = 1, 2, ..., 6) to product the POD basis Φi = (φi1 , φi2 , ..., φi6 ) with  the formula φij = Pi ϕij / λi, j (1  j  6, i = 1, 2). Finally, we compute out the ROECNFSE solutions (und , vnd ) and Qnd of the velocity and heat energy (n = 30000, i.e., at t = 3) according on Step 4 in Section 3.4 and depict them in (b)’s of Figs. 1 and 2, respectively. In order to show the advantage of the ROECNFSE format, we also compute out the CCNFSE solutions (unN , vnN ) and QnN of the velocity and heat energy (when n = 30000, i.e., at t = 3) by the CCNFSE format, Problem 5, and Remark 6 and depict them in (a)’s of Figs. 1 and 2, respectively. Both photos (a)’s and (b)’s of Figs. 1 and 2 are approximately identical, but the results of the ROECNFSE solutions are better than those of the CCNFSE solutions. Fig. 3 shows the errors between between the ROECNFSE solutions of the velocity and heat energy with the variational number of POD basis and the corresponding CCNFSE solutions when t = 3, Re = 1000, and Pr = 0.7. It also shows that the numerically computing results are accorded with the theoretical ones, because theoretical and numerical errors are (10−8 ). Moreover, by comparing the unknowns and the CPU elapsed time of the CCNFSE format with the ROECNFSE format, we can easily realize that the ROECNFSE format is superior to the CCNFSE format. First, in this experiments, the ROECNFSE format at each time node only has twelve unknowns (i.e., twelve degrees of freedom), whereas the CCNFSE format at the same time nodes includes 2 × 4 × 104 unknowns. 14

Fig. 1. (a) The CCNFSE velocity solution when t = 3. (b) The ROECNFSE velocity solution when t = 3.

(a)

(b)

Fig. 2. (a) The CCNFSE heat energy solution when t = 3. (b) The ROECNFSE heat energy solution when t = 3.

Fig. 3. The errors between the ROECNFSE solutions with different number of POD basis and the CCNFSE solutions when t = 3.

Next, when t = 3, the CPU elapsed time on Think-Pad E530 laptop for the CCNFSE format is about 3654 seconds, but the CPU elapsed time on the same laptop for the ROECNFSE format only needs about 18 seconds. These sufficiently explain that the ROECNFSE format is far superior to the CCNFSE format. It is 15

shown that the ROECNFSE format is very effective for seeking the numerical solutions of the non-stationary Boussinesq equations about the vorticity-stream functions. 4.2. The numerical experiments of channel flow with two rectangular protrusions The computational domain Θ consists of a channel with a width of 6 and a total length of 20, with two identical rectangular protrusions at the bottom and at the top of the channel. The two rectangular protrusions both have a width of 2 and a length of 4 (see Fig. 4). A structured mesh with side length x = y = 0.01 is used. Except for the inflow from the left boundary with a velocity of (u, v) = (0.1(y − 2)(8 − y), 0) (x = 0, 2  y  8) and the outflow on the right boundary with velocity of u(x, y,t) = u(19, y,t) (19  x  20, 2  y  8, 0  t  T ), all of the initial and other boundary value conditions are taken as 0. The time-step increment is also taken as Δt = 0.0001. We yet take Re = 1000 and Pr = 0.7. Under the circumstances, the theoretical errors for the CCNFSE solutions also are O(10−8 ).

Fig. 4. The computational domain and boundary conditions of flow.

Fig. 5. (a) The CCNFSE solution of the velocity when t = 4. (b) The ROECNFSE solution of the velocity when t = 4.

First, we also compute out 20 CCNFSE solutions (ωn , Qn ) (n = 1, 2, ..., 20) at the initial

20 steps 1 , ω2 , ..., ω20 and P = = ω by the CCNFSE format and use them to form two snapshot matrices P 1 2

1 2 Q , Q , ..., Q20 . And then, we compute out the eigenvalues λi, j (their orders are degressively arranged) 16

Fig. 6. (a) The CCNFSE solution of the heat energy when t = 4. (b) The ROECNFSE solution of the heat energy when t = 4.

and eigenvectors ϕij ( j = 1, 2, ..., 20) of the matrix PiT Pi , further attain the POD bases according on Steps 2  and 3 in Section 3.4. By reckoning, we also achieve that λi,7  2 × 10−8 so that we only need to choose the initial six eigenvectors ϕij ( j = 1, 2, ..., 6) to formulate the POD basis Φi = (φi1 , φi2 , ..., φi6 ) by means of  the formula φ j = Pi ϕij / λi, j (1  j  6; i = 1, 2). Finally, we find out the ROECNFSE solution (und , vnd ) and Qnd of the velocity and heat energy (n = 40000, i.e., at t = 4) with Step 4 in Section 3.4 and depict them in photos (b)’s of Figs. 5 and 6, respectively. In order to show the preponderance of the ROECNFSE format, we also find out the CCNFSE solution (unN , vnN ) and QnN of the velocity and heat energy (n = 40000, i.e., at t = 4) by the CCNFSE format, Problem 5, and Remark 6 and depict them in photo (a)’s of Figs. 5 and 6, respectively. Both photos (a)’s and (b)’s of Figs. 5 and 6 are approximately identical, but the results of the ROECNFSE solutions are far better than those of the CCNFSE solutions. However, their errors all attain O(10−8 ) and are accorded with the theoretical results. Especially, in this experiment, the ROECNFSE format at each time-level has also only 2 × 6 unknowns, whereas the CCNFSE format at the same time-level has 6 × 136 × 104 unknowns, i.e., the unknowns for the ROECNFSE format are far fewer than those for the CCNFSE format, so that the ROECNFSE format can significantly reduce the round-off error accumulating in the computing, alleviate the calculating load, shorten the CPU time required for the calculations, and improve the computational accuracy. For instance, when the durative time T = 4, the required CPU time on Think-Pad E530 laptop for the CCNFSE format is about 6372 seconds, the ROECNFSE format spends only 36 seconds in CPU on the same laptop. Therefore, it is further confirmed that the ROECNFSE format is far superior to the CCNFSE format. Remark 12. The accuracy of the ROECNFSE solutions is far higher than that of other reduced-order numerical solutions, for example, the reduced-order FE, FVE, and the FD solutions. For instance, for the reduced-order FE method in [20], the reduced-order FVE method in [21], and the reduced-order FD scheme in [37], though their space-step is taken as Δx = Δy = 0.01, the accuracy of the the reduced-order FE, FVE, or FD solutions only attains O(10−2 ), whereas our space-step is also taken as Δx = Δy = 0.01, but the accuracy of the ROECNFSE solutions can attain O(10−8 ). 17

5. Conclusions and Discussion In this paper, by using the POD technique to reduce the order of the coefficient vector for the CCNFSE method, we have studied the order reduction for the CCNFSE format of the 2D non-stationary Boussinesq equations about the vorticity-stream functions. We have established the ROECNFSE format for the 2D non-stationary Boussinesq equations, analyzed the existence, uniqueness, stability, and convergence of the ROECNFSE solutions by the matrix and vector analysis method, and also given the flowchart for solving the ROECNFSE format. Moreover, we have used two sets of numerical experiments to verify the correctness of the theoretical analysis, to explain that the ROECNFSE format is far superior to the CCNFSE format because the unknowns of the ROECNFSE format are far fewer than those of the CCNFSE format so that, compared to the CCNFSE format, the ROECNFSE format can greatly reduce the computational load, retard the round-off error accumulation, and save the CPU consuming time in the computation. Especially, the establishment of the ROECNFSE format about the vorticity-stream functions is only to use the initial few coefficient vectors of the CCNFSE vorticity function and heat energy so that the constitution for the POD basis here is far simpler and more convenient than the existing other POD-based reduced-order methods (see, e.g., [6, 13, 18, 22, 28, 29, 30]). In addition, the establishment of the ROECNFSE format and the theoretical analysis here are different from the existing POD-based reduced-order methods as mentioned above and others. Especially, the ROECNFSE model for the 2D non-stationary Boussinesq equations is first presented in this paper and is a development and improvement over the existing reduced-order methods because the ROECNFSE model has higher accuracy than other reduced-order methods, such as, the reduced-order FE and FVE methods, and the reduced-order FD scheme. Although we only study ROECNFSE method for the 2D non-stationary Boussinesq equations, the technique here can be extend to the three-dimensional problems or more complex fluid dynamics equations, even be applied in the more complex real-world engineering problems. Therefore, the technique here has important applied prospect. References [1] R.A. Adams, Sobolev Spaces, Academic Press, New York, 1975. [2] J.An, Z.D. Luo, H. Li, P. Sun, Reduced-order extrapolation spectral-finite difference scheme based on POD method and error estimation for three-dimensional parabolic equation, Front. Math. China 10(5) (2015) 1025–1040. [3] P. Benner, A. Cohen, M. Ohlberger, A.K. Willcox, Model Reduction and Approximation: Theory and Algorithm, Computational Science & Engineering, SIAM, 2017. [4] J.R. Cannon, Y. Lin, A priori L2 error estimates for finite-element methods for nonlinear diffusion equations with memory, SIAM Journal on Numerical Analysis 27(3) (1990) 595–607. [5] Z. Cai, S. McCormick, On the accuracy of the finite volume element method for diffusion equations on composite grid, SIAM J. Numer. Anal. 27(3) (1990) 636–655. [6] W. Cazemier, R.W.C.P. Verstappen, A.E.P. Veldman, Proper orthogonal decomposition and low-dimensinal models for driven cavity flows, Physics of Fluids 10(7) (1998) 1685–1699. [7] A. Cohen, Numerical Analysis of Wavelet Methods, Elsevier Science B.V. Amsterdam, Netherlands, 2003. [8] V. Girault, P.A. Raviart, Finite Element Methods for Navier-Stokes Equations: Theory and Algorithms, Springer-Verlag, Berlin Heidelberg, 1986. [9] K. Fukunaga, Introduction to Statistical Recognition, Academic Press, New York, 1990. [10] O. Guba, M. Taylor, A. St-Cyr, Optimization-based limiters for the spectral element method, J. Comput. Phys. 267(267) (2014) 176–195. [11] B.Y. Guo, Spectral Methods and Their Applications, Singapore, World Scientific, 1998. [12] B.Y. Guo, Some progress in spectral methods, Sci. China-Math. 56(12) (2013) 2411–2438. [13] J.S. Hesthaven, G. Rozza, B. Stamm, Certified Reduced Basis Methods for Parametrized Partial Differential Equations, Springer International Publishing, 2016. [14] P. Holmes, J.L. Lumley, G. Berkooz, Turbulence, Coherent Structures, Dynamical Systems and Symmetry, Cambridge University Press, Cambridge, 1996. [15] I.T. Jolliffe, Principal Component Analysis, Springer-Verlag, Berlin, 2002.

18

[16] W.P. Jones, K.R. Menziest, Analysis of the cell-centred finite volume method for the diffusion equation, J. Comput. Phys. 165(1) (2000) 45–68. [17] K. Kunisch, S. Volkwein, Galerkin proper orthogonal decomposition methods for parabolic problems, Numer. Math. 90(1) (2001) 117–148. [18] K. Kunisch, S. Volkwein, Galerkin proper orthogonal decomposition methods for a general equation in fluid dynamischs, SIAM J. Numer. Anal. 40(2) (2002) 492–515. [19] Z.D. Luo, The Foundations and Applications of Mixed Finite Element Methods, Chinese Science Press, Beijing, 2006 ( in Chinese). [20] Z.D. Luo, A POD-based reduced-order finite difference extrapolating model for the non-stationary incompressible Boussinesq equations, Adv. Differ. Equ. 2014(272) (2014) 1–18. [21] Z.D. Luo, Proper orthogonal decomposition-based reduced-order stabilized mixed finite volume element extrapolating model for the nonstationary incompressible Boussinesq equations, J. Math. Anal. Appl. 425(1) (2015) 259–280. [22] Z.D. Luo, J. Chen, I.M. Navon, X.Y. Yang, Mixed finite element formulation and error estimates based on proper orthogonal decomposition for the non-stationary Navier–Stokes equations, SIAM J. Numer. Anal. 47(1) (2008) 1–19. [23] Z.D. Luo, S.J. Jin, A reduced-order extrapolation spectral-finite difference scheme based on the POD method for 2D secondorder hyperbolic equations, Math. Model. Anal. 22(5) (2017) 569–586. [24] Z.D. Luo, H. Li, P. Sun, A fully discrete stabilized mixed finite volume element formulation for the non-stationary conduction–convection problem, J. Math. Anal. Appl. 404(1) (2013) 71–85. [25] Z.D. Luo, H. Li, Y.J. Zhou, X.M. Huang, A reduced FVE formulation based on POD method and error analysis for twodimensional viscoelastic problem, J. Math. Anal. Appl. 385(1) (2012) 310–321. [26] Z.D. Luo, H. Li, Y.J. Zhou, Z.H. Xie, A reduced finite element formulation based on POD method for two-dimensional solute transport problems, J. Math. Anal. Appl. 385(1) (2012) 371–383. [27] Z.D. Luo, Z.H. Xie, Y.Q. Shang, J. Chen, A reduced finite volume element formulation and numerical simulations based on POD for parabolic problems, J. Comput. Appl. Math. 235(8) (2011) 2098–2111. [28] Z.D. Luo, X.Z. Yang, Y.J. Zhou, A reduced finite difference scheme based on singular value decomposition and proper orthogonal decomposition for Burgers equation, J. Comput. Appl. Math. 229(1) (2009) 97–107. [29] H.V. Ly, H.T. Tran, Proper orthogonal decomposition for flow calculations and optimal control in a horizontal CVD reactor, Quart. Appl. Math. 60(4) (1989) 631–656. [30] A. Quarteroni, A. Manzoni, F. Negri, Reduced Basis Methods for Partial Differential Equations, Springer International Publishing, 2016. [31] F.M. Selten, Baroclinic empirical orthogonal functions as basis functions in an atmospheric model, J. Atmos. Sci. 54(16) (1997) 2099–2114. [32] J. Shen, T. Tang, Spectral and High-Order Methods with Applications, Science Press, Beijing, 2006. [33] L. Sirovich, Turbulence and the dynamics of coherent structures: part I-III, Quart. Appl. Math. 45(3) (1987) 561–590. [34] P. Sun, Z.D. Luo, Y.J. Zhou, Some reduced finite difference schemes based on a proper orthogonal decomposition technique for parabolic equations, Appl. Numer. Math. 60(1-2) (2010) 154–164. [35] J.P. Wang, A finite spectral finite element method for incompressible Navier-Stokes equations, Int. J. Mech. Research 3(3) (2014) 33–42. [36] J.H. Wu, The 2D Incompressible Boussinesq Equations, Peking University Summer School Lecture Notes, Beijing, 2012. [37] H. Xia, Z.D. Luo, A stabilized MFE reduced-order extrapolation model based on POD for the 2D unsteady conductionconvection problem. J. Inequal. Appl. 2017(124) (2017) 1–17. [38] W.S. Zhang, Finite Difference Methods for Patial Differential Equations in Science Computation, Higher Education Press, Beijing, 2006.

19