Reduced-Order Extrapolation Finite Volume Element Methods Based on Proper Orthogonal Decomposition

Chapter 3

Reduced-Order Extrapolation Finite Volume Element Methods Based on Proper Orthogonal Decomposition To begin this new chapter, we briefly review what we have done in the preceding two chapters: 1. In Chapter 1, we have introduced the POD-based reduced-order extrapolation finite difference (FD) (PODROEFD) schemes for the following PDEs: the two-dimensional (2D) parabolic equation; the 2D unsteady Stokes equation; the 2D shallow water equation with sediment concentration. 2. In Chapter 2, we have studied the POD-based reduced-order extrapolation finite element (FE) (PODROEFE) method or the POD-based reducedorder stabilized Crank–Nicolson (CN) extrapolation mixed FE (MFE) (PODROSCNEMFE) method for the following PDEs: the 2D viscoelastic wave equation (by PODROEFE); the 2D unsteady Burgers equation (by PODROEFE); the 2D parabolized Navier–Stokes equation (by PODROSCNEMFE). Here, in contrast to the FD and FE methods, we will develop the finite volume element (FVE) method. The FVE method is sometimes referred to as a box method (see [12]) or a generalized difference method (see [65,67]). It is regarded as a highly effective discretization tool for PDEs because it is easy to implement and offers flexibility in handling complicated domains (see [23, 24,61,153]). Most importantly, it can maintain local mass or other conservation laws, which is highly desirable in many applications. It has been widely applied to find numerical solutions of various types of PDEs (see, e.g., [4,8,12,23,24,32, 49,50,61,64,65,67,86,73,146,153,187,188]). However, the classical FVE method contains many degrees of freedom (i.e., unknowns) related to the nodes of grids. To reduce computational complexity, therefore, it is necessary to reduce the degrees of freedom of the classical FVE method. In this Chapter 3, we study the POD-based reduced-order extrapolation FVE (PODROEFVE) methods. We first introduce the construction, the theoretical analysis, and the implementation of algorithms for the Proper Orthogonal Decomposition Methods for Partial Differential Equations https://doi.org/10.1016/B978-0-12-816798-4.00008-5 Copyright © 2019 Elsevier Inc. All rights reserved.

157

158 Proper Orthogonal Decomposition Methods for Partial Differential Equations

PODROEFVE method for the 2D hyperbolic equation and the Sobolev equation and introduce the POD-based reduced-order extrapolation stabilized CN mixed FVE (PODROSCNEMFVE) method for the 2D nonstationary incompressible Boussinesq equation including the velocity vector field and the pressure field as well as the temperature field, respectively. Then, we provide numerical examples to show that the PODROEFVE methods and PODROSCNEMFVE method are more advantageous than the standard FVE methods. Furthermore, it is shown that the PODROEFVE methods and PODROSCNMEFVE method are workable and efficient for finding numerical solutions of PDEs.

3.1 POD-BASED REDUCED-ORDER EXTRAPOLATING FINITE VOLUME ELEMENT ALGORITHM FOR TWO-DIMENSIONAL HYPERBOLIC EQUATION In this section, we first establish a semidiscretized scheme with respect to time for the 2D hyperbolic equation. Then, we build a fully discretized FVE formulation directly from the semidiscretized scheme about time for the 2D hyperbolic equation. Thus, we can avoid the semidiscretized FVE formulation with respect to spatial variables for the 2D hyperbolic equation such that our method becomes simpler and more convenient than the existing methods (see, e.g., [28, 177]). Next, we establish the PODROEFVE algorithm for the 2D hyperbolic equation and provide the error estimates of solutions of the PODROEFVE algorithm, and the implementation for solving the PODROEFVE algorithm. Finally, we provide some numerical experiments to verify that the results of numerical computation are consistent with theoretical analysis. These contents are based on reference [115].

3.1.1 Classical Finite Volume Element Method for the 2D Hyperbolic Equation and Generation of Snapshots In this subsection, we first give a semidiscretized formulation with respect to time for the 2D second-order hyperbolic equation and the existence and error estimates of its solutions, and then we start directly from the semidiscretized formulation with respect to time to establish a classical fully discretized FVE formulation and to analyze errors.

3.1.1.1 The Second-Order Hyperbolic Equation Let  ⊂ R 2 be a bounded domain with piecewise smooth boundary ∂ and consider the following initial boundary value problems for the 2D hyperbolic equation in  × [0, T ].

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

Problem 3.1.1. Find u such that ⎧ ⎪ ⎨ utt − εu = f (x, y, t), (x, y, t) ∈  × (0, T ], u(x, y, t) = ϕ(x, y, t), (x, y, t) ∈ ∂ × (0, T ], ⎪ ⎩ u(x, y, 0) = ϕ (x, y), u (x, y, 0) = ϕ (x, y), (x, y) ∈ , 0 t 1

159

(3.1.1)

where utt = ∂ 2 u/∂t 2 ,  is the 2D Laplacian, ε is a positive constant, T is the total time, and source term f (x, y, t), initial value functions ϕ0 (x, y) and ϕ1 (x, y), and boundary value function ϕ(x, y, t) all are given. The hyperbolic equation is commonly used to describe acoustics and wave phenomena in nature, such as hydrodynamics, displacement problems in porous media and vibrations of a membrane, and electromagnetic wave propagation. For the sake of convenience and without loss of generality, we assume that ϕ0 (x, y), ϕ1 (x, y), and ϕ(x, y, t) are identically zero functions in the following theoretical analysis.

3.1.1.2 Semidiscretized Formulation with Respect to Time for the 2D Hyperbolic Equation The Sobolev spaces and their norms for functional analysis used in here are standard as in Chapter 1 (see [1]). Let U = H01 () be the Hilbert space on domain . The variational formulation for Problem 3.1.1 can be written as follows. Problem 3.1.2. Find u(t) : [0, T ] → U such that  (utt , w) + εa(u, w) = (f, w), ∀w ∈ U, u(x, y, 0) = 0, ut (x, y, 0) = 0, (x, y) ∈ ,

(3.1.2)

where a(u, w) = (∇u, ∇w), and (·, ·) denotes the inner product in L2 (). By using the same approach as in the proof of Theorem 2.2.1, we can deduce the existence and uniqueness of the solution for Problem 3.1.2 and obtain the following estimates (see [41,133,149]): | u(m+1) (t) |s + | u(m) (t) |s+1  Ct f H m (H s ) , 0  t  T , m = 0, 1, 2, 3, s = 0, 1,

(3.1.3)

s where | · |s is the · H m (H s ) is the norm in H m (0, T ;  seminorm in H (),  (m) s H ()), Ct = 2 exp(t)/ min{1, ε}, and u (t) is the mth order derivative of u with respect to time t. Note that  · 1 is equivalent to ∇ · 0 in H01 () by the Poincaré inequality, namely, there is a constant C0  0 such that

∇v0  v1  C0 ∇v0 , ∀v ∈ H01 ().

(3.1.4)

160 Proper Orthogonal Decomposition Methods for Partial Differential Equations

Let N be the positive integer and k = T /N denote the time step size. For any function g(x, y, t), let g n be the semidiscretized approximation with respect to time of g(x, y, t) at time tn = nk (0  n  N ). Put g n,1/2 =

g n+1 + g n−1 ¯ n g n+1 − g n ¯ ¯ n g n+1 − 2g n + g n−1 . , ∂t g = , ∂t ∂t g = 2 k k2

Then, the semidiscretized formulation with respect to time t for Problem 3.1.2 can be written as follows. Problem 3.1.3. Find un+1 ∈ U (n = 1, 2, · · · , N − 1) such that (∂¯t ∂¯t un , v) + εa(un,1/2 , v) = (f n , v) ∀v ∈ U,

(3.1.5)

u0 (x, y) = u1 (x, y) = 0, (x, y) ∈ ,

(3.1.6)

subject to

where f n = f (x, y, tn ). If ϕ0 (x, y) and ϕ1 (x, y) are nonzero functions, it is only necessary to define u0 = ϕ0 (x, y) and u1 = ϕ1 (x, y). Let ˜ n+1 , v) = 2(un+1 , v) + k 2 εa(un+1 , v), A(u F˜ (v) = (4un − 2un−1 , v) − k 2 εa(un−1 , v) + 2k 2 (f n , v). Then, Problem 3.1.3 can be rewritten as follows. Problem 3.1.4. Find un+1 ∈ U (n = 1, 2, · · · , N − 1) such that ˜ n+1 , v) = F˜ (v), ∀v ∈ U, A(u

(3.1.7)

u0 (x, y) = u1 (x, y) = 0, (x, y) ∈ .

(3.1.8)

subject to

From Problem 3.1.3 or 3.1.4, we have the following result. Theorem 3.1.1. If f ∈ L2 (0, T ; L2 ()), Problem 3.1.4, or equivalently, Problem 3.1.3 has a unique sequence of solutions un ∈ U (n = 1, · · · , N). If f ∈ H 3 (0, T ; L2 ()), then the following error estimates hold: ˜ 2 f H 3 (L2 ) , u(tn ) − un 0 + k∇(u(tn ) − un )0  Ck where C˜ = Ct T (1 + 3ε)(T ε + 2C0 )/(24ε) and n = 1, 2, · · · , N .

(3.1.9)

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

161

˜ ·) is a symmetric positive definite and bounded Proof. Since, for fixed k, A(·, bilinear function on U × U and F˜ (·) is also a bounded linear function on U , Problem 3.1.4, or equivalently, Problem 3.1.3, has a unique sequence of solutions un ∈ U (n = 1, · · · , N) from the Lax–Milgram theorem, Theorem 2.1.12 (see [31,65,80]). Taking t = tn in Problem 3.1.2 and using Taylor’s formula, we obtain (u(tn+1 ) − 2u(tn ) + u(tn−1 ), v) + = k 2 (f n , v) + (R1 , v) +

εk 2 a(u(tn+1 ) + u(tn−1 ), v) 2

εk 2 a(R2 , v), 2

∀v ∈ U,

(3.1.10)

where R1 = k 4 [u(4) (ξ1n ) + u(4) (ξ2n )]/24, R2 = k 2 [u(2) (ξ3n ) + u(2) (ξ4n )]/4, and ξin ∈ [tn−1 , tn ] (i = 1, 2, 3, 4). Set en = u(tn ) − un (n = 0, 1, 2, · · · , N ). Subtracting (3.1.7) from (3.1.10) yields (en+1 − 2en + en−1 , v) + = (R1 , v) +

εk 2 a(en+1 + en−1 , v) 2

εk 2 a(R2 , v), ∀v ∈ U. 2

(3.1.11)

Taking v = en+1 − en−1 in (3.1.11) and using the Hölder inequality gives en+1 − en 20 − en − en−1 20 +  (R1 0 +

εk 2 (∇en+1 20 − ∇en−1 20 ) 2

εk 2 R2 2 )en+1 − en−1 0 . 2

(3.1.12)

Note that en+1 − en−1 0  en+1 − en 0 + en − en−1 0 and en+1 − en−1 0  en+1 0 + en−1 0  C0 (∇en+1 0 + ∇en−1 0 ). Thus, from (3.1.12), we have en+1 − en 0 − en − en−1 0 +  R1 0 + 

εk 2 (∇en+1 0 − ∇en−1 0 ) 2C0

εk 2 R2 2 2

k 4  (4) u (ξ1n )0 + u(4) (ξ2n )0 24 εk 4  (2) + u (ξ3n )2 + u(2) (ξ4n )2 , n = 1, 2, · · · , N − 1. 8

(3.1.13)

162 Proper Orthogonal Decomposition Methods for Partial Differential Equations

Note that e1 = e0 = 0. Summing (3.1.13) from 1 to n, we obtain εk 2 (∇en+1 0 + ∇en 0 ) 2C0 n k 4  (4)  u (ξ1i )0 + u(4) (ξ2i )0 24

en+1 − en 0 +

i=1

n εk 4  (2) + u (ξ3i )2 + u(2) (ξ4i )2 . 8

(3.1.14)

i=1

Then, from (3.1.14), we have k 4  (4) u (ξ1i )0 + u(4) (ξ2i )0 24 n

en+1 0 − en 0 

i=1

+

εk 4 8

n 

u(2) (ξ3i )2 + u(2) (ξ4i )2 , n = 1, 2, · · · , N − 1.

(3.1.15)

i=1

Summing (3.1.15) from 1 to n, we obtain en+1 0 

n nk 4  (4) u (ξ1i )0 + u(4) (ξ2i )0 24 i=1

+

n εnk 4  (2) u (ξ3i )2 + u(2) (ξ4i )2 . 8

(3.1.16)

i=1

Thus, combining (3.1.14) and (3.1.16) with (3.1.3), we obtain the error estimates (3.1.9) of Theorem 3.1.1. The proof of Theorem 3.1.1 is complete.

3.1.1.3 Classical Fully Discretized Finite Volume Element Method for the 2D Hyperbolic Equation In the following, we start directly from the semidiscretized formulation with respect to time to establish the classical fully discretized FVE scheme and to analyze errors. To this end, we need to rehash the FVE theory for an FVE approximation of Problem 3.1.4 according to [65], where one can also find more detailed FVE theory. First, let h = {K} be a quasiuniform triangulation on  with h = max hK , where hK is the diameter of the triangle K ∈ h (see Definition 2.1.10 or references [21,65,80]). In order to describe the FVE method, we introduce a dual partition ∗h based on h whose elements are called the control volumes. We construct the control volume in the same way as in [65]. Let zK be the barycenter of K ∈ h . We connect zK with line segments to the midpoints of the edges of K, thus partitioning K into three quadrilaterals Kz (z = (xz , yz ) ∈ Zh (K), where

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

163

FIGURE 3.1.1 The left figure is a triangle K partitioned into three subdomains Kz . The right figure is a sample domain with dotted lines indicating the corresponding control volume Vz .

 Zh (K) are the vertices of K). Then with each vertex z ∈ Zh = K∈ h Zh (K) we associate a control volume Vz , which consists of the union of the subdomains Kz , sharing the vertex z. Finally, we obtain a set of control volumes covering the domain , which is called the dual partition ∗h of the triangulation h (see Fig. 3.1.1). We denote the set of interior vertices of Zh by Zh◦ . We call the partition ∗h regular or quasiuniform, if there exist two positive constants C1 and C2 , being independent of the spatial mesh size h and temporal mesh size k, such that C1 h2  mes(Vz )  C2 h2 , ∀Vz ∈ ∗h ,

(3.1.17)

where mes(Vz ) denotes the measure of the dual element Vz . The barycenter-type dual partition can be introduced for any FE triangulation h and leads to relatively simple calculations. Besides, if the FE triangulation h is quasiuniform, then the dual partition ∗h is also quasiuniform (see [65,80]). The trial function space Uh and the test function space U˜ h are, respectively, chosen as follows: Uh = {wh ∈ C() ∩ H01 (); wh |K ∈ P1 (K), ∀K ∈ h },

U˜ h = wh ∈ L2 (); wh |Vz ∈ P0 (Vz )(∀Vz ∈ ∗h ), wh |Vz = 0 (if Vz ∩ ∂ = ∅) , where Pm (e) (m = 0, 1) are the spaces of polynomials of degree  m on closed subset e (e = K or Vz ). It is obvious that Uh ⊂ U . For w ∈ U , let h w be the interpolation projection of w onto the trial function space Uh . By the interpolation theory of Sobolev spaces (see Theorem 2.1.10 or references [21,31,65,80]), we have | w − h w |m  Ch2−m | w |2 , m = 0, 1, if w ∈ H 2 (),

(3.1.18)

where C is a generic positive constant independent of the spatial mesh size h and temporal mesh size k.

164 Proper Orthogonal Decomposition Methods for Partial Differential Equations

U˜ h is spanned by the following basis func-

In fact, the test function space tions:  1, φz (x, y) = 0,

(x, y) ∈ Vz , z ∈ Zh◦ . elsewhere,

For any wh ∈ U˜ h , we have

wh =

wh (z)φz .

z∈Zh◦

For w ∈ U , let ∗h w be the L2 interpolation projection of w onto the test space U˜ h , namely,

w(z)φz . (3.1.19) ∗h w = z∈Zh◦

By the interpolation theory again (see Theorem 2.1.10 or [65,80]), we have w − ∗h w0  Ch | w |1 .

(3.1.20)

Moreover, the interpolation projection ∗h satisfies the following. Lemma 3.1.2. If vh ∈ Uh , then  (vh − ∗h vh )dxdy = 0, K ∈ h , K

vh − ∗h vh Lr ()  ChK vh W 1 ,r() , 1  r  ∞. Though the trial function space Uh satisfies Uh ⊂ U just like the FE methods, the test function space, however, has the property that U˜ h ⊂ Uh . As in the case of nonconforming FE methods, this is due to the loss of continuity of the functions in U˜ h on the boundary of two neighboring elements. So the bilinear form a(u, w) must be revised accordingly. Using the idea of nonconforming FE methods, the integral on the whole domain is written as a sum of integrals on the dual elements Vz as follows: 

 wudxdy = wudxdy 

Vz ∈ ∗h

=− +

Vz



∇u∇wdxdy

Vz ∈ ∗h

Vz

Vz ∈ ∗h

∂Vz





 ∂u ∂u wdy − wdx , ∂x ∂y

(3.1.21)

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

165

 where ∂Vz denotes a line integral, in the counterclockwise direction, on the boundary ∂Vz of the dual element. So we have a(u, w) =

 Vz ∈ ∗h

∇u∇wdxdy

Vz





−

∂Vz

Vz ∈ ∗h

 ∂u ∂u wdy − wdx , ∂x ∂y

∀u, w ∈ U.

(3.1.22)

Since U˜ h is the space of piecewise constant functions with the characteristic functions of the dual elements Vz as the basis functions, we have a(uh , wh ) = − =−





wh

Vz ∈ ∗h

∂Vz

Vz ∈ ∗h

z∈Zh◦



∂uh ∂uh dy − dx ∂x ∂y



wh (z)a(u ˜ h , φz ), ∀uh ∈ Uh , ∀wh ∈ U˜ h ,

(3.1.23)

 ∂uh ∂uh where a(u ˜ h , φz ) = dy − dx , z = (xz , yz ). ∂x ∂y ∂Vz Now, the classical fully discretized FVE formulation for Problem 3.1.2 can be stated as follows. 



∈ Uh (n = 1, 2, · · · , N − 1) such that Problem 3.1.5. Find un+1 h (∂¯t ∂¯t unh , ∗h vh ) + εa(uh

n,1/2

, ∗h vh ) = (f n , ∗h vh ), ∀vh ∈ Uh ,

(3.1.24)

subject to u0h (x, y) = u1h (x, y) = 0, (x, y) ∈ .

(3.1.25)

In order to discuss the existence, uniqueness, and error estimates of solutions for Problem 3.1.5, we need the following two lemmas of properties about a(uh , ∗h vh ), cf. [64,65]. Lemma 3.1.3. The bilinear form a(uh , ∗h vh ) is symmetric, bounded, and positive definite on Uh , i.e., a(uh , ∗h vh ) = a(vh , ∗h uh ) = a(uh , vh ), ∀uh , vh ∈ Uh , satisfying a(uh , vh )  ∇uh 0 ∇vh 0 ,

a(uh , uh )  ∇uh 20 ,

∀uh , vh ∈ Uh .

166 Proper Orthogonal Decomposition Methods for Partial Differential Equations

Proof. Because ∗h is the L2 interpolation projection operator, it follows from the definition of a(uh , ∗h vh ), Lemma 3.1.2, and Green’s formula that   



∂uh ∂uh ∗ ∗ a(uh , h vh ) = − h vh (z) dy − dx ∂x ∂y ∂Vz Vz ∈ ∗h z∈Zh◦   



∂uh ∂uh =− vh (z) dy − dx ∂x ∂y ∂Vz ◦ ∗ Vz ∈ h z∈Zh

= a(uh , vh ). For the rest, the properties of boundedness and positive definiteness about a(uh , vh ) are obvious. Lemma 3.1.4. Following the same notation as before, we have the following: i. (uh , ∗h u¯ h ) = (u¯ h , ∗h uh ), ∀uh , u¯ h ∈ Uh ; ii. for any u ∈ H m () (m = 0, 1) and vh ∈ Uh , | (u, vh ) − (u, ∗h vh ) | Chm+n um vh n , n = 0, 1; iii. define ||| uh |||0 = (uh , ∗h uh )1/2 ; then ||| · |||0 is equivalent to  · 0 on Uh , i.e., there exist two positive constants C3 and C4 such that C3 uh 0 ||| uh |||0  C4 uh 0 , ∀uh ∈ Uh . Proof. Because ∗h is the L2 interpolation projection operator, it follows from (2.1.1) and Green’s formula that 



∗ ∗ (uh , h u¯ h ) = h u¯ h (z) uh dxdy Vz ∈ ∗h z∈Zh◦

=



Vz ∈ ∗h

=

u¯ h (z)

z∈Zh◦









uh dxdy K

 u¯ h dxdy

uh (z) K

K∈ h z∈Zh◦

=

uh dxdy Vz

u¯ h (z)

K∈ h z∈Zh◦

=

Vz



 u¯ h dxdy

uh (z)

Vz ∈ ∗h z∈Zh◦

Vz

= (∗h u¯ h , uh ). For any u ∈ H m ()2 (m = 0, 1) and vh ∈ Uh , by using the Hölder inequality, Lemma 3.1.2, and the interpolation theory (see Theorem 2.1.10 or [65,80]), we

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

167

have | (u, vh ) − (u, ∗h vh ) | um vh − ∗h vh −m  Chm+n um vh n , n = 0, 1. For uh ∈ Uh , from (2.1.1), we have ||| uh |||20 = (uh , ∗h uh ) = =



z∈Zh◦

Vz ∈ ∗h



∗h uh (z) 

=



Vz

=

⎛ ⎝

K∈ h

 uh dxdy

uh (z) K

K∈ h z∈Zh◦ ∩∂K



uh dxdy Vz

uh dxdy

uh (z)

Vz ∈ ∗h z∈Zh◦





⎞2



uh (z)⎠ mes(K).

z∈Zh◦ ∩∂K

On the other hand, from (2.1.1), we have

 || uh ||20 = u2h dxdy K∈ h K



=

K∈ h

⎡

⎢1 mes(K) ⎣ 4

z∈Zh◦ ∩∂K

⎛ u2h (z) +

1⎝ 4

z∈Zh◦ ∩∂K

⎞2 ⎤ ⎥ uh (z)⎠ ⎦ .

Thus, by using the Cauchy–Schwarz inequality, we obtain that ||| · |||0 is equivalent to  · 0 on Uh . In order to derive error estimates for the solutions of Problem 3.1.5, we need to introduce a generalized elliptic projection. For given solutions un ∈ U (n = 0, 1, 2, · · · , N ) of Problem 3.1.3, let Rh : U → Uh satisfy Rh u1 = h u1 , Rh u0 = h u0 , and (Rh un+1 − un+1 − 2(Rh un − un ) + Rh un−1 − un−1 , vh ) εk 2 (∇(Rh un+1 − un+1 + Rh un−1 − un−1 ), ∇vh ) 2 = 0, ∀vh ∈ Uh , n = 1, 2, · · · , N. +

(3.1.26)

Then, by using standard FE arguments (see, e.g., [21,31,80]), we have ∇Rh un 0  C∇ui 0 , n = 1, 2, · · · , N,

(3.1.27)

168 Proper Orthogonal Decomposition Methods for Partial Differential Equations

Rh un − un 0 + k

n

∇(Rh un − un )0  Ckh2

i=1

n

ui 2 ,

i=1

n = 1, 2, · · · , N.

(3.1.28)

The following theorem holds for Problem 3.1.5. Theorem 3.1.5. Under the hypotheses of Theorem 3.1.1, Problem 3.1.5 has a unique set of solutions {unh }N n=1 ⊂ Uh satisfying the following stability result: unh 0 + k 2 ∇unh 0  C˜ 1 k 2

n

f (ti )0 ,

(3.1.29)

i=1

where C˜ 1 = C4 nk/ min{C3 , ε/(2C0 C4 )}  C4 T / min{C3 , ε/(2C0 C4 )}. Further, if k 2 = O(h), the following error bounds between the solution u(t) of Problem 3.1.2 and the solutions unh of Problem 3.1.5 hold:  u(tn ) − unh 0 + k∇(u(tn ) − unh )0  C h2 + k 2 , n = 1, 2, · · · , N, (3.1.30) where C is a constant independent of h but dependent on the data ε, T , and f of Problem 3.1.2. Proof. Using the same technique as Theorem 3.1.1, by Lemmas 3.1.3 and 3.1.4, we easily prove that Problem 3.1.5 has a unique sequence of solutions {unh }N n=1 ⊂ n+1 n−1 Uh . Taking vh = uh + uh in Problem 3.1.5 and using the Hölder inequality and Lemmas 3.1.3 and 3.1.4, we have εk 2 2 + un−1 ∇(un+1 h h )0 2  k 2 ||| f n |||0 ||| un+1 + un−1 |||0 . h h

|||20 − ||| un−1 |||20 + ||| un+1 h h

(3.1.31)

n+1 + un−1 |||0  C4 un+1 + un−1 + un−1 Note that ||| un+1 h h h h 1  C4 C0 ∇(uh h )0 . From (3.1.31), we obtain

εk 2 n−1 (∇un+1 h 0 − ∇uh 0 ) 2C4 C0 εk 2 ||| un+1 |||0 − ||| un−1 |||0 + ∇(un+1 + un−1 h h h h )0 2C4 C0  k 2 ||| f n |||0 . (3.1.32)

||| un+1 |||0 − ||| un−1 |||0 + h h

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

169

Note that u0h = u1h = 0. Summing (3.1.32) from 1 to n yields |||0 + ||| unh |||0 + ||| un+1 h  k2

n

εk 2 n (∇un+1 h 0 + ∇uh 0 ) 2C4 C0

||| f i |||0 ,

(3.1.33)

i=1

which gives (3.1.29). Subtracting Problem 3.1.5 from Problem 3.1.3 by taking v = vh and using Lemma 3.1.3, we have − 2(un − unh ) + (un−1 − un−1 (un+1 − un+1 h h ), vh ) εk 2 + un−1 − un−1 (∇(un+1 − un+1 h h ), ∇vh ) 2 ∗ = k 2 (f n , vh − ∗h vh ) − (un+1 − 2unh + un−1 h h , vh − h vh ), +

∀vh ∈ Uh , n = 1, 2, · · · , N − 1.

(3.1.34)

Set un − unh = (un − Rh un ) + (Rh un − unh ) = ρ n + en . Then, by (3.1.26) and (3.1.34), we obtain εk 2 a(en+1 + en−1 , vh ) 2 = (Rh un+1 − un+1 − 2(Rh un − un ) + Rh un−1 − un−1 , vh )

(en+1 − 2en + en−1 , ∗h vh ) +

+ (Rh un+1 − un+1 − 2(Rh un − un ) + Rh un−1 − un−1 , ∗h vh − vh ) − 2(un − unh ) + un−1 − un−1 + (un+1 − un+1 h h , vh ) ∗ − 2(un − unh ) + un−1 − un−1 + (un+1 − un+1 h h , h vh − vh )

εk 2 (∇(Rh un+1 − un+1 + Rh un−1 − un−1 ), ∇vh ) 2 εk 2 + + un−1 − un−1 (∇(un+1 − un+1 h h ), ∇vh ) 2 = k 2 (f n , vh − ∗h vh ) − (ρ n+1 − 2ρ n + ρ n−1 , vh − ∗h vh ) +

+ (un+1 − 2un + un−1 , ∗h vh − vh ), ∀vh ∈ Uh .

(3.1.35)

Taking vh = en+1 − en−1 in (3.1.35) and using Lemma 3.1.2, we have εk 2 (∇en+1 20 − ∇en−1 20 ) 2 = k 2 (f n − ∗h f n , en+1 − en−1 − ∗h (en+1 − en−1 ))

||| en+1 − en |||20 − ||| en − en−1 |||20 +

− (ρ n+1 − 2ρ n + ρ n−1 , en+1 − en−1 − ∗h (en+1 − en−1 ))

170 Proper Orthogonal Decomposition Methods for Partial Differential Equations

  C h2 k 2 f n 1 + h2 ∇(ρ n+1 − 2ρ n + ρ n−1 )0  + h2 ∇(un+1 − 2un + un−1 )0 ∇(en+1 − en−1 )0 .

(3.1.36)

By Taylor’s formula and (3.1.28), we obtain ∇(ρ n+1 − 2ρ n + ρ n−1 )0  C∇(un+1 − 2un + un−1 )0  Ck 2 f H 3 (L2 ) . (3.1.37) By the inverse inequality, we have ∇(en+1 − en−1 )0  Ch−1 en+1 − en−1 0  Ch−1 (en+1 − en 0 + en − en−1 0 ). Since ∇(en+1 − en−1 )0  ∇en+1 0 + ∇en−1 0 , by simplifying (3.1.36), we obtain Ch(||| en+1 − en |||0 − ||| en − en−1 |||0 ) εk 2 (∇en+1 0 − ∇en−1 0 ) + 2    Ck 2 h2 f L∞ (H 1 ) + f H 3 (L2 ) .

(3.1.38)

If k 2 = O(h), by summing (3.1.38) from 1 to n and using Lemma 3.1.4, we obtain en+1 − en 0 + ∇en+1 0 + ∇en 0  Cnhk 2  Ck 3 .

(3.1.39)

Applying the triangle inequality to (3.1.39) yields en+1 0 − en 0 + ∇en+1 0  Ck 3 .

(3.1.40)

Summing (3.1.40) from 1 to n − 1, we obtain en 0 +

n

∇ei 0  Cnk 3  Ck 2 .

(3.1.41)

i=1

Combining (3.1.28) with (3.1.41) and Theorem 3.1.1, we have obtained (3.1.30).

Remark 3.1.1. If ϕ0 (x, y) and ϕ1 (x, y) are nonzero functions, then it is necj essary to assume ϕj (x, y) ∈ W 2,3 () and uh (x, y) = h ϕj (x, y) (j = 0, 1), in order to get the same results as Theorem 3.1.5. Thus, if f (x, y, t), ε, ϕ0 (x, y), and ϕ1 (x, y), the triangulation parameter h, the time step increment k, and the trial function space Uh are given, we can obtain a sequence of fully discretized FVE solutions unh ∈ Uh (n = 1, 2, · · · , N ) by solving Problem 3.1.5. From them, we can choose the first L (in general, L  N , for example, L = 20, N = 200) solutions unh (1  n  L) as the snapshots.

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

171

Remark 3.1.2. In the same way as we noted before, when one computes an actual problem, one may obtain a collection of snapshots from physical system trajectories by drawing samples from experiments and interpolation (or data assimilation).

3.1.2 Formulating the POD Basis and Establishing the Reduced-Order Algorithm for the 2D Hyperbolic Equation In this subsection, we recapitulate the POD method (see Sections 2.2.4 and 2.3.4 or References [90,108,120,123]) and establish the PODROEFVE algorithm for the 2D second-order hyperbolic equation. For unh (x, y) (n = 1, 2, · · · , L) in Section 3.1.3, let Wi (x, y) = uih (x, y) (1  i  L) and V = span{W1 , W2 , · · · , WL } and refer to V as the space generated by the snapshots {Wi }L i=1 , at least one of l which is assumed to be a nonzero function. Let {ψj }j =1 denote an orthonormal basis of V with l = dim{V}. Then each member Wi of the set of snapshots can be expressed as Wi =

l

(Wi , ψj )U ψj , i = 1, 2, · · · , L,

(3.1.42)

j =1

where (Wi , ψj )U = (∇uih , ∇ψj ) and (·, ·) is the L2 -inner product. Definition 3.1.1 (POD method and POD basis). The POD method consists in finding an orthonormal basis ψj (i = 1, 2, · · · , l), (1  d  l) such that the mean square error between the elements Wi (1  i  L) and the corresponding dth partial sum of (3.1.42) are minimized on average, i.e.,

1

(Wi , ψj )U ψj min Wi − d {ψj }j =1 L L

d

i=1

j =1

2

,

(3.1.43)

U

subject to (ψr , ψj )U = δrj , 1  r  d, 1  j  r,

(3.1.44)

where Wi 2U = ∇uih 20 . A solution {ψj }dj =1 satisfying (3.1.43) and (3.1.44) is called a POD basis of rank d. We formulate the correlation matrix A = (Aij )L×L ∈ R L×L corresponding to the snapshots {Wi }L i=1 by Aij = (Wi , Wj )U /L. Since the matrix A is positive semidefinite and has rank l, the solution of (3.1.43) and (3.1.44) can be found.

172 Proper Orthogonal Decomposition Methods for Partial Differential Equations

In addition, the following result holds (see Propositions 2.2.4 and 2.3.8 or [90, 108,120,123]). Proposition 3.1.6. Let λ1  λ2  · · ·  λl > 0 denote the positive eigenvalues of A and v 1 , v 2 , · · · , v l the associated orthonormal eigenvectors. Then a POD basis of rank d  l is given by 1 (W1 , W2 , · · · , WL ) · v i , ψi = √ Lλi

1  i  d  l.

(3.1.45)

Furthermore, the following identity holds: L d l



1

Wi − (Wi , ψj )U ψj 2U = λj . L j =1

i=1

(3.1.46)

j =d+1

Let U d = span {ψ1 , ψ2 , · · · , ψd }. For uh ∈ Uh , define a Ritz projection operator R d : Uh → U d as follows: (∇R d uh , ∇wd ) = (∇uh , ∇wd ), ∀wd ∈ U d . According to functional analysis (see [143]), there is an extension R h : U → Uh of R d such that R h |Uh = R d : Uh → U d defined by (∇R h u, ∇wh ) = (∇u, ∇wh ),

∀wh ∈ Uh ,

(3.1.47)

where u ∈ U . Due to (3.1.47), the operator R h is well defined and bounded (see Theorem 2.1.11 or [80]) such that ∇(R h u)0  ∇u0 ,

∀u ∈ U.

(3.1.48)

We have the following inequality (see Lemma 2.2.5 or [90,123]): u − R h u0  Ch∇(u − R h u)0 ,

∀u ∈ U.

(3.1.49)

The following result holds (see Lemmas 2.2.5 and 2.3.9, or [21,90,108,120, 123]). Lemma 3.1.7. For every d (1  d  l), the projection operator R d satisfies L l

1

∇(uih − R d uih )20  λj , L i=1

(3.1.50)

j =d+1

where uih ∈ V (i = 1, 2, · · · , L) are the solutions of Problem 3.1.5. Moreover, if un ∈ U (n = 1, 2, · · · , N ) are the solutions to Problem 3.1.3, then the following estimates hold: un − R h un s  Ch2−s un 2  Ch2−s f n 0 ,

s = 0, 1.

(3.1.51)

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

173

By using U d , we can formulate the PODROEFVE method for Problem 3.1.2 as follows. Problem 3.1.6. Find und ∈ U d (n = 1, 2, · · · , N ) such that und = R d unh =

d

(∇uih , ∇ψj )ψj , n = 1, 2, · · · , L,

(3.1.52)

j =1 n,1/2 (∂¯t ∂¯t und , ∗h vd ) + εa(ud , vd )

= (f n , ∗h vd ), ∀vd ∈ U d , n = L, L + 1, · · · , N − 1. (3.1.53) Let und = α1n ψ1 + α2n ψ2 + · · · + αdn ψd (n = 1, 2, · · · , N). By using the property of ∗h , Problem 3.1.6 can be rewritten as follows. Problem 3.1.7. Find (α1n , α2n , · · · , αdn )T ∈ R d (n = 1, 2, · · · , N ) such that αin = (∇unh , ∇ψi ), i = 1, 2, · · · , d, n = 1, 2, · · · , L, d

(3.1.54)

αjn+1 [2Bij + k 2 εCij ] = Gn,n−1 + Fin , i

j =1

i = 1, 2, · · · , d, n = L, · · · , N − 1,

(3.1.55)

, and Fin (i = 1, 2, · · · , d, n = L, L + 1, · · · , N − 1) are where Bij , Cij , Gn,n−1 i given by Bij =





ψi (xz , yz )

Cij =



ψi (xz , yz ) ∂Vz

 ∂ψj (x, y) ∂ψj (x, y) dx − dy , ∂y ∂x

d



4Bij αjn − (k 2 εCij + 2Bij )αjn−1 ,

j =1

Fin = 2k 2





Vz ∈ ∗h

= Gn,n−1 i

ψj (x, y)dxdy, Vz

Vz ∈ ∗h



Vz ∈ ∗h

 ψi (xz , yz )

f (x, y, tn )dxdy. Vz

Remark 3.1.3. Since Uh is the space of piecewise linear functions, if h is a uniformly regular triangulation, then the number of total degrees of freedom, i.e., the number of unknowns, for Problem 3.1.5 at each time level is Nh (where Nh is the number of vertices of triangles in h ; see [31] or [80]), whereas the number of total degrees of freedom for Problem 3.1.6 at each time level (when n > L) is d (d  l  L). For real-world scientific and engineering problems,

174 Proper Orthogonal Decomposition Methods for Partial Differential Equations

the number of vertices of triangles in h is of the order of tens of thousands or even a few hundred millions, while d here is only the number of a few maximal eigenvalues so d is really small (for example, in Section 3.1.5, d = 6, while Nh  4 × 106 ). In particular, it needs no repetitive computation and one only uses the first few L given solutions of Problem 3.1.5 to extrapolate all the other (n > L) solutions. This is completely different from the existing reduced-order models (see, e.g., [2,26,62,63,88,90–93,100,109,112,113,118,122–124,155]).

3.1.3 Error Estimates of the Reduced-Order Solutions for the 2D Hyperbolic Equation In this subsection, we resort to the classical FVE method to derive the error estimates of the solutions for Problem 3.1.6. We have the following main results. Theorem 3.1.8. Under the same assumptions of Theorem 3.1.5, Problem 3.1.6 has a unique sequence of solutions und ∈ U d such that ∇und 0  ∇unh 0  C˜ 1

n

f (ti )0 , 1  n  L,

(3.1.56)

i=1

und 0

+k

2

n

uid 1  Cf L∞ (L2 ) ,

L + 1  n  N. (3.1.57)

i=1

Consequently, the sequence of solutions to Problem 3.1.6 is stable. If k 2 = O(h), then we have the following error estimates between the reduced-order solutions und to Problem 3.1.6 and the classical FVE solutions unh to Problem 3.1.5: unh − und 0 + k ⎛  C ⎝k 2 L

n

∇(uih − uid )0

i=i

⎞1/2

l

λj ⎠

, n = 1, 2, · · · , L,

(3.1.58)

j =d+1

unh − und 0 + k

n

i=0

∇(uih − uid )0 ⎛

 Ckh(n − L) + C ⎝k 2 L

l

⎞1/2 λj ⎠

, n = L + 1, · · · , N.

(3.1.59)

j =d+1

Proof. Note that U d ⊂ Uh . For fixed k, let n+1 n+1 ∗ 2 Ad (un+1 d , vd ) = 2(ud , h vd ) + k εa(ud , vd ), n−1 ∗ 2 n ∗ 2 Fd (vd ) = (4und − 2un−1 d , h vd ) + 2k (f , h vd ) − k εa(ud , vd ).

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

175

Thus, Problem 3.1.6 can be restated as follows. Problem 3.1.8. Find und ∈ U d (n = 1, 2, · · · , N ) such that d Ad (un+1 d , vd ) = Fd (vd ), ∀vd ∈ U , n = L, L + 1, · · · , N − 1,

(3.1.60)

subject to und

= R d unh

=

d

(∇uih , ∇ψj )ψj , n = 1, 2, · · · , L.

(3.1.61)

j =1

By Lemmas 3.1.3 and 3.1.4, we obtain Ad (vd , vd ) = 2 ||| vd |||20 +k 2 ε∇vd 20  αund 21 , ∀vd ∈ U d ,

(3.1.62)

where α = min{2C3 , k 2 ε}, namely, the bilinear form Ad (·, ·) on the left-hand side of (3.1.60) is a symmetric positive definite bilinear functional on U d , and Ad (·, ·) and Fd (·) are obviously bounded. Thus, by the Lax–Milgram theorem (see Theorems 2.1.12 and 2.1.13 or [31,80]), Problem 3.1.8, or equivalently, Problem 3.1.6, has a unique sequence of solutions und ∈ U d (n = 1, 2, · · · , N). When n = 1, 2, · · · , L, from (3.1.48), we have ∇und 20 = ∇R d unh 20 = ∇R h unh 20  ∇unh 20 .

(3.1.63)

From (3.1.63) and Theorem 3.1.5, we immediately obtain (3.1.56). −und . Taking vd = un+1 − When n = L+1, L+2, · · · , N −1, let Qnd = un+1 d d n+1 n + un − un−1 = Qn + Qn−1 in Problem 3.1.6, we have un−1 = u − u d d d d d d d k2ε n+1 + un−1 − un−1 a(un+1 d d , ud d ) 2 = k 2 (f n , ∗h (Qnd + Qn−1 (3.1.64) d )), L  n  N − 1.

n−1 ∗ n (Qnd − Qn−1 d , h (Qd + Qd )) +

By using Lemmas 3.1.3 and 3.1.4 and the Hölder and Cauchy–Schwarz inequalities, from (3.1.64), we obtain k2ε n−1 2 2 [∇un+1 d 0 − ∇ud 0 ] 2 = k 2 (f n , ∗h (Qnd + Qn−1 d ))    k 2 f (tn )0 ||| Qnd − Qn−1 ||| . (3.1.65) 0 d

|||20 + ||| Qnd |||20 − ||| Qn−1 d

Note that ||| Qnd − Qn−1 |||0 ||| un+1 |||0 + ||| un−1 |||0 or ||| Qnd − Qn−1 |||0  d d d d ||| Qnd |||0 + ||| Qn−1 ||| . From (3.1.65), we obtain 0 d ||| Qnd |||0 − ||| Qn−1 |||0 + d  k 2 f (tn )0 .

 k2ε  n−1 ∇un+1 d 0 − ∇ud 0 2 (3.1.66)

176 Proper Orthogonal Decomposition Methods for Partial Differential Equations

Note that H 2 (L2 ) →→ L∞ (L2 ). By Taylor’s formula, Theorem 3.1.5, and Lemma 3.1.7, we have d L−1 ||| QL−1 |||0  C4 R d uL h − R uh 0 d L L  C4 [R d uL h − uh 0 + uh − u(tL )0 + u(tL ) − u(tL−1 )0 L−1 + u(tL−1 ) − uL−1 − R d uL−1 h 0 + uh h 0 ]

 C(k + h2 )f L∞ (L2 ) .

(3.1.67)

Thus, summing (3.1.66) from L to n (n  N − 1) and using (3.1.56), we have 

k2ε  n 2  + ∇u  f (ti )0 ∇un+1  k 0 0 d d 2 n

||| Qnd |||0 +

i=L

k2ε L + C(k + h2 )f L∞ (L2 ) + (∇uL−1 h 0 + ∇uh 0 ). 2

(3.1.68)

Since ||| Qnd |||0 =||| un+1 − und |||0 ||| un+1 |||0 − ||| und |||0 , summing (3.1.68) d d from L to n (n  N − 1) yields ||| un+1 |||0 +k 2 ε d

n+1

uid 1

i=L

 k2

j n



f (ti )0 + C(n − L)(k + h2 )f L∞ (L2 )

j =L i=L

(n − L)k 2 ε L (∇uL−1 h 0 + ∇uh 0 ) 2 n

f (ti )0 + C(n − L)(k + h2 )f L∞ (L2 )  nk 2 +

i=L

(n − L)k 2 ε L + (∇uL−1 h 0 + ∇uh 0 )  Cf L∞ (L2 ) , 2

(3.1.69)

which gives (3.1.57). Since U d ⊂ Uh , subtracting Problem 3.1.6 from Problem 3.1.5 by taking vh = vd ∈ U d , we obtain (∂¯t ∂¯t (unh − und ), ∗h vd ) + kεa(uh

n,1/2

n,1/2

− ud

, vd )

= 0, ∀vd ∈ U , n = L, L + 1, · · · , N − 1, d

(3.1.70)

subject to unh − und = unh − R d unh ,

n = 1, 2, · · · , L.

(3.1.71)

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

177

For n = 1, 2, · · · , L, by Lemma 3.1.7 and (3.1.49), we obtain unh − und 0 + k

n

∇(uih − uid )0

i=1

= unh − R d unh 0 + k

n

∇(uih − R d ui )0

i=1

 Ch∇(unh − R d unh )0 + k

n

∇(uih − R d ui )0

i=1

 Ck

n

∇(uih − R d ui )0

i=1

!

 Ck n ⎛

n

"1/2 ∇(uih

i=1

 C ⎝k 2 L

l

−R u

d i

)20

⎞1/2 λj ⎠

, n = 1, 2, · · · , L.

(3.1.72)

j =d+1

Let θ n = unh − und = (unh − R d unh ) + (R d unh − und ) = ρ n + n . On the one hand, by using Lemmas 3.1.3 and 3.1.4, we have ε (∂¯t ∂¯t θ n , ∗h (θ n+1 − θ n−1 )) + a(θ n+1 + θ n−1 , θ n+1 − θ n−1 ) 2 ε n n−1 ¯ n ¯ ¯ ¯ = (∂t θ − ∂t θ , ∂t θ + ∂t θ n−1 ) + a(θ n+1 + θ n−1 , θ n+1 − θ n−1 ) 2 ε =||| ∂¯t θ n |||20 − ||| ∂¯t θ n−1 |||20 + [∇θ n+1 20 − ∇θ n−1 20 ] 2 ε n 2 n−1 2 ||| ∂¯t θ |||0 − ||| ∂¯t θ |||0 + [∇θ n+1 20 − ∇θ n−1 20 ]. (3.1.73) 2 On the other hand, from the definition of R d , Theorems 3.1.8 and 3.1.9, the Hölder inequality, the system of error equations (3.1.70), and (3.1.49), if k 2 = O(h), we have ε (∂¯t ∂¯t θ n , ∗h (θ n+1 − θ n−1 )) + a(θ n−1 + θ n−1 , θ n+1 − θ n−1 ) 2 ε n ∗ n+1 n−1 ¯ ¯ = (∂t ∂t θ , h (ρ − ρ )) + a(θ n+1 + θ n−1 , ρ n+1 − ρ n−1 ) 2 = (∂¯t ∂¯t θ n , ∗h (ρ n+1 − ρ n−1 )) ε + a(ρ n+1 + ρ n−1 + n+1 + n−1 , ρ n+1 − ρ n−1 ) 2 ε = (∂¯t ∂¯t θ n , ∗h (ρ n+1 − ρ n−1 )) + a(ρ n+1 + ρ n−1 , ρ n+1 − ρ n−1 ) 2

178 Proper Orthogonal Decomposition Methods for Partial Differential Equations

1 ε = (∂¯t θ n+1 − ∂¯t θ n , ∗h (ρ n+1 − ρ n−1 )) + a(ρ n+1 + ρ n−1 , ρ n+1 − ρ n−1 ) k 2 Ch [||| ∂¯t θ n |||20 + ||| ∂¯t θ n−1 |||20 ] + Ch−1 (||| ρ n+1 |||20 + ||| ρ n−1 |||20 )  k ε + (∇ρ n+1 20 − ∇ρ n−1 20 ) 2  Ck[||| ∂¯t θ n |||20 + ||| ∂¯t θ n−1 |||20 ] + Ch−1 (||| ρ n+1 |||20 + ||| ρ n−1 |||20 ) ε (3.1.74) + (∇ρ n+1 20 − ∇ρ n−1 20 ). 2 Combining (3.1.73) with (3.1.74) yields ε ||| ∂¯t θ n |||20 − ||| ∂¯t θ n−1 |||20 + [∇θ n+1 20 − ∇θ n−1 20 ] 2  Ck[||| ∂¯t θ n |||20 + ||| ∂¯t θ n−1 |||20 ] + Ch−1 (||| ρ n+1 |||20 ε + ||| ρ n−1 |||20 ) + (∇ρ n+1 20 − ∇ρ n−1 20 ). (3.1.75) 2 Note that, by Lemma 3.1.7 and Theorem 3.1.5, we have ∇ρ i 0  ∇(uih − ui )0 + ∇(ui − R h ui )0 + ∇R h (ui − uih )0  Ch and ρ i 0  Ch2 . Thus, summing (3.1.75) from L to n (L  n  N − 1) and using Lemma 3.1.7 and Proposition 3.1.6, we have ε ||| ∂¯t θ n |||20 + [∇θ n+1 20 + ∇θ n 20 ] 2 n

 Ck ||| ∂¯t θ i |||20 +C(n − L)h3 + Ch2 .

(3.1.76)

i=L

If k is sufficiently small so that Ck  1/2, from (3.1.76), we obtain ||| ∂¯t θ n |||20 +ε[∇θ n+1 20 + ∇θ n 20 ]  Ck

n−1

||| ∂¯t θ i |||20 +C(n − L)h3 + Ch2 .

(3.1.77)

i=L

Applying Lemma 1.4.1 (the discrete Gronwall lemma) to (3.1.77), we have ||| ∂¯t θ n |||20 +[∇θ n+1 20 + ∇θ n 20 ]    C (n − L)h3 + Ch2 exp[Ck(n − L − 1)]    C (n − L)h3 + h2 ,

(3.1.78)

that is,  1/2 . ||| ∂¯t θ n |||0 +∇θ n+1 0  C (n − L)h3 + h2

(3.1.79)

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

179

Since ||| ∂¯t θ n |||0  [||| θ n+1 |||0 − ||| θ n |||0 ]/k, from (3.1.79), we obtain  ||| θ n+1 |||0 +k∇θ n+1 0 ||| θ n |||0 +Ckh + Ck 2 h (n − L).

(3.1.80)

If k 2 = O(h), summing (3.1.80) from L to n − 1 and noting that θ L 0   1/2 # C k 2 L lj =d+1 λj , we have θ n 0 + k

n

∇θ i 0

i=0

⎛

l

 Ckh(n − L) + Ck 2 h(n − L)3/2 + C ⎝k 2 L ⎛  Ckh(n − L) + C ⎝k 2 L

l

λj ⎠

j =d+1

⎞1/2 λj ⎠

⎞1/2

(3.1.81)

,

j =d+1

where n = L + 1, L + 2, · · · , N . This completes the proof of Theorem 3.1.8. Combining Theorem 3.1.8 with Theorems 3.1.5, we obtain the following. Theorem 3.1.9. Under the same assumptions of Theorem 3.1.8, we have the error estimates between the solutions u(tn ) at tn (n = 1, 2, · · · , N ) for Problem 3.1.2 and the solutions und for the PODROEFVE Problem 3.1.6 as follows: u(tn ) − und 0 + k

n

∇(u(ti ) − uid )0

i=0

⎛

 C(h2 + k 2 ) + C ⎝k 2 L

l

⎞1/2 λj ⎠

, n = 1, 2, · · · , L,

j =d+1

u(tn ) − und 0 + k

n

i=L

∇(u(ti ) − uid )0 ⎛

 C(h2 + k 2 ) + Ckh(n − L) + C ⎝k 2 L

l

⎞1/2 λj ⎠

, n = L + 1, · · · , N.

j =d+1

Remark 3.1.4. We make the following additional comments on Theorems 3.1.8 and 3.1.9: 1. The errors in Theorems 3.1.8 and 3.1.9 provide a suggestion what to choose as the number d of POD bases and the number L of snapshots, namely, for k 2 = O(h) and k  1, those d and L should be required to satisfy

180 Proper Orthogonal Decomposition Methods for Partial Differential Equations

 # 1/2 = O(k). Therefore, the number L of snapshots should L lj =d+1 λj √ not be too large. It would be best to satisfy, e.g., L < 5 (usually L = 20) so that it does not affect total errors and can reduce the computational load solving the eigenvalue problem. 2. The errors kh(n − L) (L + 1  n  N ) in Theorems 3.1.8 and 3.1.9 are caused by the extrapolation of the algorithm. They can serve as a guide for the timely renewal of POD basis generated # with the newest solutions of Problem 3.1.6. That is, if h(n − L)  (L lj =d+1 λj )1/2 (L + 1  n  N ), und (1  n  N ) is the sequence of solutions of the reduced-order Problem 3.1.6 satisfying the accuracy requirement; else, if h(n − L) > # (L lj =d+1 λj )1/2 (L + 1  n  N ), it is necessary to renew the POD basis, (i = 1, 2, · · · , L), and to update generated with new snapshots Wi = un−i d the PODROEFVE algorithm Problem 3.1.6. 3. The condition k 2 = O(h) is not strict; it is only a working condition for the theoretical error analysis and may change from case to case.

3.1.4 Implementation of the Reduced-Order Algorithm for the 2D Hyperbolic Equation In the following, we address implementation for solving the PODROEFVE algorithm Problem 3.1.6 for the 2D second-order hyperbolic equation. The flowchart consists of seven steps. Step 1. Classical FVE computation and extraction of snapshots For given the coefficient ε, the source term f (x, y, t), the initial states ϕ0 (x, y) and ϕ1 (x, y), the boundary value function ϕ(x, y, t), the total time duration T , and the desirable accuracy δ = max{k 2 , h2 } = k 2 (since condition k 2 = O(h)), solving the following classical FVE formulation for the first few L steps (L  N = T /k, usually, take L = 20), we have (∂¯t ∂¯t unh , ∗h vh ) + εah (uh

n,1/2

, ∗h vh ) = (f n , ∗h vh ),

∀vh ∈ Uh , n = 1, 2, · · · , L − 1, subject to u0h (x, y) = h ϕ0 (x, y), u1h (x, y) = h ϕ1 (x, y), (x, y) ∈ , and unh = h ϕ(x, y, tn ), (x, y) ∈ ∂, 1  n  L, generates the collection of snapshots Wi (x, y) = uih , i = 1, 2, · · · , L,

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

181

which may, indeed, also be samples from the physical system trajectories, from experiments with interpolation (or data assimilation), or previously acquired results. Step 2. Formulation of correlation matrix A The correlation matrix A = (Aij )L×L , where Aij = (Wi , Wj )U /L = j (∇uih , ∇uh )/L and (·, ·) is the L2 -inner product. Step 3. Computing eigenvalues and eigenvectors of A Find the eigenvalues λ1  λ2  · · ·  λl > 0 (l = dim{W1 , W2 , · · · , WL }) j j j and corresponding eigenvectors v j = (a1 , a2 , · · · , aL )T (j = 1, 2, · · · , l) of A. Step 4. Determination of the number of POD bases # The number d of POD bases should satisfy L lj =d+1 λj  δ 2 . Step 5. Choice of POD basis # j i  The POD bases ψj = L i=1 ai uh / Lλj (j = 1, 2, · · · , d). Step 6. Computation of PODROEFVE solutions Solve the following system of equations that only includes d degrees of freedom at every time level: αin = (∇Wn , ∇ψi ) = (∇unh , ∇ψi ), i = 1, 2, · · · , d, n = 1, 2, · · · , L, 

Bij = ψi (xz , yz ) ψj (x, y)dxdy, i, j = 1, 2, · · · , d, Vz

Vz ∈ ∗h

Cij =



ψi (xz , yz )

Vz ∈ ∗h

Fin = 2k 2







Vz ∈ ∗h

∂Vz

 ∂ψj (x, y) ∂ψj (x, y) dx − dy , i = 1, 2, · · · , d, ∂y ∂x



ψi (xz , yz )

f (x, y, tn )dxdy, Vz

i = 1, 2, · · · , d, n = L, · · · , N − 1, d

n−1 n 2 4B , Gn,n−1 = α − (k εC + 2B )α ij j ij ij j i j =1

i = 1, 2, · · · , d, n = L, · · · , N − 1, d

αjn+1 (2Bij + k 2 εCij ) = Gn,n−1 + Fin , i

j =1

i = 1, 2, · · · , d, n = L, L + 1, · · · , N − 1, to obtain (α1n , α2n , · · · , αdn )T ∈ R d (n = 1, 2, · · · , N). We then obtain the solu# tions und = di=1 αin ψi (x, y) (n = 1, 2, · · · , L, L + 1, · · · , N ).

182 Proper Orthogonal Decomposition Methods for Partial Differential Equations

Step 7. Renewal of POD basis and circulation or end # If h(n−L)  (L lj =d+1 λj )1/2 (L+1  n  N ), then und (n = 1, 2, · · · , N) form a sequence of solutions for Problem # 3.1.6 satisfying desirable accuracy. End. Else, namely, if h(n − L) > (L lj =d+1 λj )1/2 (L + 1  n  N ), let Wi = uid (i = n − L, n − L + 1, · · · , n − 1). Return to Step 2. − und 0  und − Remark 3.1.5. Indeed, Step 7 could be done so that if un−1 d n+1 n ud 0 (n = L, L + 1, · · · , N − 1), then ud (n = 1, 2, · · · , N) are the solutions − und 0 < for Problem 3.1.6 satisfying desirable accuracy. Else, namely, if un−1 d n+1 n i ud − ud 0 (n = L, L + 1, · · · , N − 1), let Wi = ud (i = n − L, n − L − 1, · · · , n − 1). Return to Step 2.

3.1.5 Numerical Experiments for the 2D Hyperbolic Equation In the following, we present some numerical experiments for the 2D hyperbolic equation with piecewise continuous initial value and boundary value functions to validate the practicability and efficiency of the PODROEFVE algorithm. Moreover, we verify that the results of numerical computations are consistent with theoretical analysis. We consider the 2D hyperbolic equation on a computational domain 2 cm × 2 cm, i.e.,  = [0, 2] × [0, 2]. We first divide the closed domain  = {(x, y); 0  x  2, 0  y  2} into 2000 × 2000 small squares with side length x = y = 0.001, then link the diagonal of the square to divide each square into two triangles along the same direction, and adopt local refining meshes such that the meshes nearby (x, 2) (0  x  2) are one-third of the meshes nearby √ (x, 0) (0  x  2), which makes up the triangularization h . Thus h = 2 × 0.001. In order to satisfy k 2 = O(h), we take the time step as k = 0.04, T = 500k = 20. The dual decomposition ∗h is taken as barycenter dual decomposition, i.e., the barycenter of the right triangle K ∈ h is taken as the node of the dual decomposition. We take ε = 1, the source term f (x, y, t) = 0, the initial value functions ⎧ ⎪ (x − 2)2 + 0.5, 1  x  2, y = 2, ⎪ ⎪ ⎨ 0.5, 0.5  x  1, y = 2, ϕ0 (x, y) = ⎪ 0, 0.0  x  0.5, y = 2, ⎪ ⎪ ⎩ 0, elsewhere,  2(x − 2), 1  x  2, y = 2, ϕ1 (x, y) = 0, elsewhere, and the following boundary value function:  ϕ0 (x, y), (x, y) ∈ ∂, t = 0, ϕ(x, y, t) = n−1 uh (x, y), (x, y) ∈ ∂, t = tn , n = 1, 2, · · · , n.

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

183

FIGURE 3.1.2 The left graphic is the classical FVE solution unh when t = 10, while the right one is the PODROEFVE solution und when d = 6 and t = 10.

FIGURE 3.1.3 The left graphic is the classical FVE solution unh when t = 20, while the right one is the PODROEFVE solution und when d = 6 and t = 20.

We first find the numerical solutions unh by the classical FVE formulation, namely, Problem 3.1.5 for t = 10 and t = 20, depicted graphically in the left columns of Figs. 3.1.2 and 3.1.3, respectively. Next, we choose only the first 20 states at time t = 0.04, 0.08, · · · , 0.80 to constitute a collection of snapshots. # After computing, we have achieved, for d = 6 and k = 0.04, the error 1/2  3 × 10−3 in Theorem 3.1.8, which shows that we need (20k 2 20 j =7 λj ) only to take six main POD bases to expand into subspaces U d . Finally, we find two PODROEFVE solutions at t = 250k = 10 and t = 500k = 20 with the PODROEFVE algorithm Problem 3.1.6 following the seven steps of Section 3.1.4, but when reaching t = 15 it is only necessary to renew the POD basis once, as depicted in the right columns of Figs. 3.1.2 and 3.1.3. Every pair of two figures in Figs. 3.1.2 and 3.1.3 exhibits close similarity, but the PODROEFVE solutions are in many ways better than the classical FVE solutions due to the use and content of more initial numerical solutions, namely, six POD bases and the reduction of the truncation error accumulation in the computational process. Fig. 3.1.4 shows the errors between 20 solutions und of the PODROEFVE algorithm Problem 3.1.6 with 20 different numbers of POD bases and the solutions unh of the classical FVE formulation, namely, Problem 3.1.5 at t = 10 and 20. It shows that the results for numerical examples are consistent with those obtained for the theoretical cases since theoretical and numerical errors do not exceed 3 × 10−3 .

184 Proper Orthogonal Decomposition Methods for Partial Differential Equations

FIGURE 3.1.4 For t = 10 and 20, the absolute errors between the classical FVE solutions and the POD-based reduced-order FVE solutions with different numbers of POD bases for a group of 20 snapshots.

Comparing the classical FVE formulation, Problem 3.1.5, with the PODROEFVE algorithm, Problem 3.1.6, using six POD bases and implemented by the numerical simulations when the total time T = 20, we find that, for the classical FVE formulation, Problem 3.1.5, which includes more than 4 × 106 degrees of freedom, the required computing time is approximately 24 minutes on a laptop, while for the PODROEFVE algorithm, Problem 3.1.6, with six POD basis elements, which has only six degrees of freedom, the corresponding time duration is less than 6 seconds on the same laptop, so the required computing time to solve the classical FVE formulation, Problem 3.1.5, is about 240 times that of the PODROEFVE algorithm, Problem 3.1.6, with six POD bases, while the errors between their respective solutions do not exceed 3 × 10−3 . Thus, the time consuming calculations and resource demands for the computational process could be greatly reduced and the truncated error (TE) accumulation in the computational process is also lowered. We have shown that finding the approximate solutions for the 2D hyperbolic equation with the PODROEFVE algorithm, Problem 3.1.6, is computationally effective. Remark 3.1.6. In this Section 3.1, we have first developed the semidiscretized formulation with respect to time for the 2D hyperbolic equation. Thereafter, we have established the classical fully discretized FVE formulation that starts directly from the semidiscretized formulation with respect to time and circumvents the semidiscretized FVE formulation with respect to the space variable such that the process of argumentation could be significantly simplified. Next, we have employed the POD method to establish the PODROEFVE algorithm for the 2D hyperbolic equation, analyzed the errors between the classical FVE solutions and the PODROEFVE solutions, provided the implementation of algorithm for solving the PODROEFVE algorithm, and discussed theoretically the relation between the number of snapshots and the number of solutions at all time instants. Finally, we checked the correctness of theoretical results by means of some numerical experiments.

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

185

The major difference between the PODROEFVE algorithm here and the existing POD-based reduced-order methods consists in that most existing PODbased reduced-order formulations (see, e.g., [2,26,62,63,88,90–93,100,109,112, 113,118,122–124,155]) employ numerical solutions obtained from the classical numerical methods on the total time span [0, T ] to formulate the POD bases and establish the reduced-order models and then recompute the solutions on the same time span [0, T ], which actually leads to repeating computations on the same time span [0, T ]; while the PODROEFVE algorithm here only uses the first few given numerical solutions on a very short time span [0, T0 ] (T0  T ) as snapshots to formulate the POD basis and establish the PODROEFVE algorithm with much fewer degrees of freedom but sufficiently high accuracy for finding the numerical solutions on the total time span [0, T ]. Thus, the PODROEFVE algorithm has sufficiently adopted the advantage of the POD method, namely, utilizing the given data (on a very short time span [0, T0 ] and T0  T ) to forecast (or infer) future physical phenomena (on the time span [T0 , T ]). Therefore, the method here has many potential applications to dynamical processes, including big data.

3.2 POD-BASED REDUCED-ORDER EXTRAPOLATION FINITE VOLUME ELEMENT ALGORITHM FOR THE TWO-DIMENSIONAL SOBOLEV EQUATION In this section, the POD technique is used to treat the classical FVE formulation for the 2D Sobolev equation. A PODROEFVE algorithm with much fewer degrees of freedom and sufficiently high accuracy is developed for the 2D Sobolev equation. The error estimates with respect to the norm in H01 () between the PODROEFVE solutions and the classical FVE solutions are provided. The implementation for solving the PODROEFVE algorithm is given. By comparing the numerical results of the PODROEFE algorithm, the classical FVE formulation, the PODROEFE formulation, the classical FE formulation, the PODROEFD scheme, and the classical FD scheme for the 2D Sobolev equation, it is shown (see, e.g., Fig. 3.2.4 in Section 3.2.5) that the PODROEFVE algorithm is a highly effective numerical method. Moreover, it is shown that the PODROEFVE algorithm is feasible and efficient for solving the 2D Sobolev equation. These results were published in [99].

3.2.1 The Classical Finite Volume Method for the 2D Sobolev Equation 3.2.1.1 Physical Background and Generalized Solution for the 2D Sobolev Equation Let  ⊂ R 2 be a bounded domain with piecewise smooth boundary ∂. Consider the following initial and boundary value problems for the 2D Sobolev equation on  × [0, T ].

186 Proper Orthogonal Decomposition Methods for Partial Differential Equations

Problem 3.2.1. Find u such that ut − εut − γ u = f, (x, y, t)∈ × (0, T ],

(3.2.1)

subject to the boundary condition u(x, y, t) = 0, (x, y, t) ∈ ∂ × (0, T ]

(3.2.2)

and the initial condition u(x, y, 0) = ϕ0 (x, y), (x, y) ∈ ,

(3.2.3)

where ut = ∂u/∂t, ε and γ are two positive constants, T is the total time duration, f (x, y, t) is the source term, ϕ0 (x, y) is the initial value function, and f (x, y, t) and ϕ0 (x, y) all are smooth enough to ensure the subsequent analysis validity. For the sake of convenience and without loss of generality, we assume that ϕ0 (x, y) is the zero function in the subsequent theoretical analysis. The Sobolev equation, i.e., Problem 3.2.1, is useful to many practical engineering fields (see [147,148]); for example, it can be used to describe the fluid flow penetrating rocks, soils, or different viscous media. It usually involves a complex computing domain, initial and boundary value functions, or source terms that are dependent on the media. Analytic solutions for the practical Sobolev equation are usually not acquirable, so one must resort to numerical solutions. The FVE method (see [23,24,61,153]) is also one of the most effective numerical methods for the 2D Sobolev equation, as it has the following advantages. First, it conserves mass and total energy. Second, it has higher accuracy and is more suitable for computations involving complicated boundary conditions than the FD scheme. Third, it has the same accuracy as the FE method but is simpler and more convenient to apply than the FE method. A semidiscretized FVE formulation (see [25]) and a fully discretized FVE formulation (see [66]) for the 2D Sobolev equation are already available. They are more effective than its FE formulations and its FD schemes. However, they also include many degrees of freedom. Thus, the truncation errors in the computational process will accumulated quickly so that the classical FVE solutions can produce more deviations after just a few computing steps. Therefore, in this section, we utilize the POD method to establish the PODROEFVE algorithm for the 2D Sobolev equation. To this end, we first establish a variational formalism. Let U = H01 () be a Hilbert space on domain . Then, the variational formulation for Problem 3.2.1 can be written as follows. Problem 3.2.2. Find u(t) : [0, T ] → U such that  (ut , w) + εa(ut , w) + γ a(u, w) = (f, w), ∀w ∈ U, u(x, y, 0) = 0, (x, y) ∈ ,

(3.2.4)

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

187

where a(u, w) = (∇u, ∇w) and (·, ·) denotes the inner product either in L2 () or in [L2 ()]2 . By using the same approach as in the proof of Theorem 2.2.1, we can deduce the existence and uniqueness of solution for Problem 3.2.2 and obtain the following estimates (see also [13,147,148]): utt L2 (H 1 ) + ∇ut 0  Cft L2 (L2 ) , where  · H m (H l ) is the norm in H m (0, T ; H l ()) and C =

√

2/ min{1, 2ε, γ }.

3.2.1.2 Semidiscretized Solutions in Time for the 2D Sobolev Equation Let N be a positive integer and let k = T /N denote the time step size. For any function g(x, y, t), we define g n = g(x, y, tn ) at time tn = nk (0  n  N ). Then Problem 3.2.2 has the following time-discretized formulation. Problem 3.2.3. Find un ∈ U (n = 1, 2, · · · , N) such that (un , v) + εa(un , v) + γ ka(un , v) = k(f (tn ), v) + (un−1 , v) + εa(un−1 , v) + (R1 , v) + εa(R1 , v), ∀v ∈ U, (3.2.5) subject to u0 (x, y) = u(x, y, 0) = 0, (x, y) ∈ ,

(3.2.6)

where R1n = ∂t un − unt = O(k 2 ∂ 2 u/∂t 2 ). If ϕ0 (x, y) is a nonzero function, it is only necessary to define u0 = ϕ0 (x, y). Then, the semidiscretized formulation with respect to time t for Problem 3.2.3 is written as follows. Problem 3.2.4. Find un ∈ U (n = 1, 2, · · · , N ) such that (un , v) + εa(un , v) + γ ka(un , v) = k(f (tn ), v) + (un−1 , v) + εa(un−1 , v), ∀v ∈ U,

(3.2.7)

subject to u0 (x, y) = u(x, y, 0) = 0, (x, y) ∈ .

(3.2.8)

From Problem 3.2.4, we have the following. Theorem 3.2.1. If f ∈ L2 (0, T ; H 1 ()), Problem 3.2.4 has a unique sequence of solutions un ∈ U (n = 1, 2, · · · , N ). Under the further assumption that u ∈ H 2 (0, T ; H 1 ()), we have the following error estimates: u(tn ) − un 1  Ckuntt 1 , n = 1, 2, · · · , N,

(3.2.9)

188 Proper Orthogonal Decomposition Methods for Partial Differential Equations

where C = β0−1 γ −1 max{1, ε} and β0 is the constant in the Poincaré inequality (Theorem 2.1.4). Proof. Since the left-hand side of Problem 3.2.4 is a symmetric positive definite and bounded bilinear functional on U , by the Lax–Milgram theorem, i.e., Theorem 2.1.12, we deduce that Problem 3.2.4 has a unique sequence of solutions un ∈ U (n = 1, · · · , N ). By subtracting (3.2.5) of Problem 3.2.4 from (3.2.7) of Problem 3.2.2 by taking t = tn and v = u(tn ) − un , using Taylor’s formula, we have γ a(u(tn ) − un , u(tn ) − un ) = (R(t), u(tn ) − un ) + εa(R(t), u(tn ) − un ), t ∈ [tn−1 , tn ],

(3.2.10)

t where R(t) = tn−1 utt (s)(t − s)ds. Thus, when u ∈ H 2 (0, T ; H01 ()), by the Hölder inequality, we have γ ∇(u(tn ) − un )20 = γ a(u(tn ) − un , u(tn ) − un ) = (R(t), u(tn ) − un ) + εa(R(t), u(tn ) − un )  kutt 0 u(tn ) − un 0 + kε∇utt 0 ∇(u(tn ) − un )0  max{1, ε}kutt 1 u(tn ) − un 1 .

(3.2.11)

By the Poincaré inequality (Theorem 2.1.4), we find that  · 1 is equivalent to ∇(·)0 in H01 (), i.e., β0 v1  ∇v0 (for any v ∈ H01 ()). Thus, from (3.2.11), we obtain (3.2.9), where C = β0−1 γ −1 max{1, ε}.

3.2.1.3 Fully Discretized Finite Volume Element Solutions for the 2D Sobolev Equation In order to derive the fully discretized FVE formulation for Problem 3.2.2, it is necessary to introduce an FVE approximation for the spatial variables of Problem 3.2.4 as we did in Section 3.1.3. Let h = {K} be a quasiuniform triangulation of  with h = max hK , where hK is the diameter of the triangle K ∈ h , and ∗h be a dual partition based on h whose elements are called control volumes as first given in Section 3.1.3. The trial function space Uh and test function space U˜ h are, respectively, defined as follows:

Uh = uh ∈ C() ∩ H01 (); uh |K ∈ P1 (K), ∀K ∈ h ,

U˜ h = wh ∈ L2 (); wh |Vz ∈ P0 (Vz )(∀Vz ∈ ∗h ), wh |Vz = 0(Vz ∩ ∂ = ∅) , where Pl (l = 0, 1) are the spaces of polynomials of degree l on K. Using the same technique as in Section 3.1.3, we can establish the classical fully discretized FVE formulation for Problem 3.2.2 as follows.

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

189

Problem 3.2.5. Find unh ∈ Uh (n = 1, 2, · · · , N) such that (unh , ∗h vh ) + εa(unh , ∗h vh ) + γ ka(unh , ∗h vh ) n−1 ∗ ∗ = k(f (tn ), ∗h vh ) + (un−1 h , h vh ) + εa(uh , h v), ∀vh ∈ Uh ,

u0h (x, y) = 0, (x, y) ∈ ,

(3.2.12) (3.2.13)

where a(uh , wh ) = −



wh (z)a(u ˜ h , φz ), ∀uh ∈ Uh , ∀wh ∈ U˜ h ,

Vz ∈ ∗h

(3.2.14)

  ∂uh ∂uh where a(u ˜ h , φz ) = dy − dx , z = (xz , yz ), and ∂Vz denotes the ∂x ∂y ∂Vz line integral, in the counterclockwise direction, on the boundary ∂Vz of the dual element. 



For Problem 3.2.5, we have the following. Theorem 3.2.2. If the solution u ∈ H 2 (0, T ; W 2,3 ()) of Problem 3.2.2, we have the following stability result of solutions for Problem 3.2.5 if ϕ0 = 0: ! unh 1

 C hϕ0 1 + ϕ0 0 + k

n

" f (ti )0 .

(3.2.15)

i=1

Further, if k = O(h), we have the following error estimates between the solution u of Problem 3.2.2 and the solutions unh of Problem 3.2.5: u(tn ) − unh 1  C (h + k) , n = 1, 2, · · · , N,

(3.2.16)

where C is a generic constant independent of h and k but dependent on the other data ε, γ , and f of Problem 3.2.2. Proof. Since the left-hand side of Problem 3.2.5 constitutes a symmetric positive definite and bounded bilinear function on Uh , it follows from the Lax– Milgram theorem, i.e., Theorem 2.1.12 (or Theorem 2.1.13) that Problem 3.2.5 has a unique set of solutions {unh }N n=1 ⊂ Uh . By taking vh = uh in Problem 3.2.5 and using Lemmas 3.1.3 and 3.1.4, we immediately obtain (3.2.15). By taking v = ∗h vh in Problem 3.2.4, we obtain ((un , ∗h vh ) + εa(un , ∗h vh ) + kγ a(un , ∗h vh ) = k(f (tn ), ∗h vh ) + (un−1 , ∗h vh ) + εa(un−1 , ∗h vh ), ∀vh ∈ Uh .

(3.2.17)

190 Proper Orthogonal Decomposition Methods for Partial Differential Equations

By subtracting Eq. (3.2.12) of Problem 3.2.5 from (3.2.17) and using Lemma 3.1.3, we obtain the following error equation: (un − unh , ∗h vh ) + εa(un − unh , ∗h vh ) + kγ a(un − unh , ∗h vh ) ∗ n−1 ∗ − un−1 = (un−1 − un−1 h , h vh ) + εa(u h , h vh ), ∀vh ∈ Uh .

(3.2.18)

Because kγ a(Rh un −unh , ∗h (Rh un −unh )) = kγ a(Rh un −unh , Rh un −unh )  0, from the properties of Ritz projection, Lemmas 3.1.3 and 3.1.4, and (3.2.18), we have ||| Rh un − unh |||20 +ε || ∇(Rh un − unh ) ||20  (Rh un − unh , ∗h (Rh un − unh )) + εa(Rh un − unh , ∗h (Rh un − unh ))  (Rh un − unh , ∗h (Rh un − unh )) + εa(Rh un − unh , ∗h (Rh un − unh )) + kγ a(Rh un − unh , ∗h (Rh un − unh )) = (un − unh , ∗h (Rh un − unh )) + εa(un − unh , ∗h (Rh un − unh )) + kγ a(un − unh , ∗h (Rh un − unh )) + (Rh un − un , ∗h (Rh un − unh )) ∗ n n n−1 ∗ n n = (un−1 − un−1 − un−1 h , h (Rh u − uh )) + εa(u h , h (Rh u − uh )) ∗ n n n−1 ∗ n n − un−1 = (Rh un−1 − un−1 h , h (Rh u − uh )) + εa(Rh u h , h (Rh u − uh ))

+ (un−1 − Rh un−1 , ∗h (Rh un − unh )) 1 2  |||20 +ε∇(Rh un−1 − un−1 ||| Rh un−1 − un−1 h h )0 2  + ||| Rh un − unh |||20 +ε∇(Rh un − unh )20

+ Ch−1 un−1 − Rh un−1 20 + Ch ||| Rh un − unh |||20 .

(3.2.19)

By simplifying (3.2.19), we get ||| Rh un − unh |||20 +ε∇(Rh un − unh )20 2 ||| Rh un−1 − un−1 |||20 +ε∇(Rh un−1 − un−1 h h )0

+ Ch3 + Ch ||| Rh un − unh |||20 .

(3.2.20)

By summing (3.2.20) from 1 to n, using Theorem 2.1.4, for k = O(h) and Ch  1/2, we have Rh un − unh 20 + ε∇(Rh un − unh )20  Cnkh2 + Ck

n−1

Rh ui − uih 20 .

i=0

(3.2.21) By applying Theorem 2.4.1 to (3.2.21), we have Rh un − unh 21  CT h exp(Cnk)  Ch.

(3.2.22)

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

191

By combining (3.2.22) and Theorem 3.2.1 with Theorem 2.1.20, we obtain (3.2.16). Remark 3.2.1. If ϕ0 (x, y) is a nonzero function, then it is necessary to assume ϕ0 (x, y) ∈ W 2,3 () and let uh (x, y, 0) = h ϕ0 (x, y); we can get the same results as in Theorem 3.2.2. Thus, if f (x, y, t), ε, γ , ϕ0 (x, y), time step size k, the triangulation parameter h, and the trial function space Uh are given, we can obtain a solution time sequence unh ∈ Uh by solving Problem 3.2.5. We then choose √ the first L (in general, L  N and L < 5, for example, L = 20, N = 200) solutions unh (1  n  L) as snapshots from N solutions unh (1  n  N ). Remark 3.2.2. When we compute actual problems, we may obtain an ensemble of snapshots from physical system trajectories by drawing samples from experiments. For example, when we solve a “real-world” Sobolev equation, we can use the samples drawing from physical systems in evolution to construct the ensemble of snapshots and then reconstruct the POD basis for the ensemble of snapshots by using the POD method to be given in the following, and finally the trial function space Uh is replaced with the subspace generated with the POD basis in order to derive a reduced-order model with lower dimensions. Thus, the future evolution of the physical system can be quickly simulated, which is a nice feature for real-life applications.

3.2.2 Formulation of the POD Basis and the Reduced-Order Algorithm for the 2D Sobolev Equation In this subsection, we will employ the POD method to formulate the POD basis (see Sections 2.2.4 and 2.3.4, or [90,108,120,123]) and to establish the PODROEFVE algorithm for the 2D Sobolev equation. For unh (x, y) (n = 1, 2, · · · , L) in Section 3.2.1, let Wi (x, y) = uih (x, y) (1  i  L) and V = span{W1 , W2 , · · · , WL }

(3.2.23)

and refer to V as the space generated by the snapshots {Wi }L i=1 , at least one of l which is a nonzero vector. Let {ψj }j =1 denote an orthonormal basis of V with l = dim{V}. Then each member of the set can be expressed as Wi =

l

(Wi , ψj )U ψj , i = 1, 2, · · · , L,

(3.2.24)

j =1

where (Wi , ψj )U = (∇uih , ∇ψj ); (·, ·) is the L2 -inner product. Definition 3.2.1 (POD Method and POD Basis). The POD method involves searching for the orthonormal basis ψj (i = 1, 2, · · · , l) such that for every d

192 Proper Orthogonal Decomposition Methods for Partial Differential Equations

(1  d  l) the mean square error between the elements Wi (1  i  L) and the corresponding dth partial sum of (3.2.24) is minimized on average, i.e.,

1

(Wi , ψj )U ψj min Wi − {ψj }dj =1 L L

d

i=1

j =1

2

,

(3.2.25)

U

subject to (ψr , ψj )U = δrj , 1  r  d, 1  j  r,

(3.2.26)

where Wi 2U = ∇uih 20 . A set of solutions {ψj }dj =1 of (3.2.25) and (3.2.26) is called the POD basis of rank d. We formulate the Gramian matrix A = (Aij )L×L ∈ R L×L corresponding to the snapshots {Wi }L i=1 by  1 Aij = ∇Wi (x, y) · ∇Wj (x, y)dxdy. (3.2.27) L  Since the matrix A is positive semidefinite and has rank l, the set of solutions {ψj }dj =1 of (3.2.25) and (3.2.26) can be found and we have the following results (see Proposition 2.2.4 or [90,108,120,123]). Proposition 3.2.3. Let λ1  λ2  · · ·  λl > 0 denote the positive eigenvalues of A and v 1 , v 2 , · · · , v l the associated orthonormal eigenvectors. Then a POD basis of rank d  l is given by 1 (W1 , W2 , · · · , WL ) · v i , ψi = √ Lλi

1  i  d  l.

(3.2.28)

Furthermore, the following error formula holds: L d l



1

2 Wi − (Wi , ψj )U ψj U = λj . L i=1

j =1

(3.2.29)

j =d+1

Let U d = span {ψ1 , ψ2 , · · · , ψd }. For uh ∈ Uh , define a Ritz projection P d : Uh → U d as follows: (∇P d uh , ∇wd ) = (∇uh , ∇wd ), ∀wd ∈ U d .

(3.2.30)

Then, by functional analysis (see, e.g., [143]), there is an extension P h : U → Uh of P d such that P h |Uh = P d : Uh → U d , denoted by (∇P h u, ∇wh ) = (∇u, ∇wh ), ∀wh ∈ Uh ,

(3.2.31)

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

193

where u ∈ U . Due to (3.2.31), the operator P h is well defined and bounded, i.e., ∇(P h u)0  ∇u0 , ∀u ∈ U.

(3.2.32)

Thus, we have the following inequality (see Lemma 2.2.5 or [90,123]): u − P h u0  Ch∇(u − P h u)0 , ∀u ∈ U.

(3.2.33)

The following lemma holds (see [21,90,108,120,123]). Lemma 3.2.4. For every d (1  d  l), the projection operator P d satisfies L l

1

∇(uih − P d uih )20  λj , L i=1

(3.2.34)

j =d+1

where uih ∈ V (i = 1, 2, · · · , L) are the solutions of Problem 3.2.5. Further, if u(t) ∈ H 2 () (0  t  T ) is the solution to Problem 3.2.2, we have the following results (see [31,80]) for the Ritz projection P h defined by (3.2.31): | u(tn ) − P h u(tn ) |s  Ch2−s , n = 1, 2, · · · , N, s = 0, 1.

(3.2.35)

Thus, using U d , we can proceed to formulate the PODROEFVE algorithm for the 2D Sobolev equation as follows. Problem 3.2.6. Find und ∈ U d (n = 1, 2, · · · , N ) such that und = P d unh =

d

(∇unh , ∇ψj )ψj , n = 1, 2, · · · , L,

(3.2.36)

j =1

(und , ∗h vd ) + εa(und , ∗h vd ) + γ ka(und , ∗h vd ) n−1 ∗ ∗ d = k(f (tn ), ∗h vd ) + (un−1 d , h vd ) + εa(ud , h vd ), ∀vd ∈ U , n = L + 1, L + 2, · · · , N, (3.2.37)

where unh (n = 1, 2, · · · , L) are the first L solutions for Problem 3.2.5. Write und = α1n ψ1 + α2n ψ2 + · · · + αdn ψd (n = 1, 2, · · · , N). By using ∗h , Problem 3.2.6 can be rewritten as follows. Problem 3.2.7. Find (α1n , α2n , · · · , αdn )T ∈ R d (n = 1, 2, · · · , N ) such that αjn = (∇unh , ∇ψj ), j = 1, 2, · · · , d, n = 1, 2, · · · , L, d

j =1

αjn [Bij + (ε + γ k)Cij ] = Fin +

d

(3.2.38)

αjn−1 (Bij + εCij ),

j =1

i = 1, 2, · · · , d, n = L + 1, L + 2, · · · , N,

(3.2.39)

194 Proper Orthogonal Decomposition Methods for Partial Differential Equations

where

Cij =

ψi (xz , yz )

Vz ∈ ∗h

Fin =







∂Vz

ψi (xz , yz )

 ∂ψj (x, y) ∂ψj (x, y) dx − dy , ∂y ∂x

f (x, y, tn )dxdy,

(3.2.40)

Vz

Vz ∈ ∗h

Bij =







ψj (x, y)dxdy, i, j = 1, 2, · · · , d.

ψi (xz , yz ) Vz

Vz ∈ ∗h

Remark 3.2.3. Since Uh is the space of piecewise linear functions and h is a uniformly regular triangulation, the total degrees of freedom, i.e., the number of unknowns, for Problem 3.2.5 at each time level is Nh (where Nh is the number of vertices of triangles in h ; see [31,80]), whereas the total number of degrees of freedom for Problem 3.2.6 at each time level is d (d  l  L). For practical scientific and engineering problems, the number Nh of vertices of triangles in h is often more than tens of thousands or even a few hundred millions, but our d here is only the number of a few main eigenvalues, so d is very small (e.g., in Section 3.2.5, d = 6, but Nh = 2000 × 2000 = 4 × 106 ). Thus, Problem 3.2.6 constitutes the PODROEFVE algorithm for the 2D Sobolev equation, i.e., Problem 3.2.1. In particular, Problem 3.2.6 uses only the first few given L solutions of Problem 3.2.5 to extrapolate all other (n > L) solutions, so it does not involve repetitive computations. Therefore, it is completely different from the existing reduced-order formulations (see, e.g., [4,8,12,23,24,32,49,61,64,65, 67,86,73,146,153,187,188]) based on the POD method.

3.2.3 Error Estimations of the Reduced-Order Solutions for the 2D Sobolev Equation In the following, we resort to the classical FVE method to deduce the error estimates of solutions for Problem 3.2.6. We have the following main results for Problem 3.2.6. Theorem 3.2.5. Under the same assumptions of Theorem 3.2.2, Problem 3.2.6 has a unique set of solution und ∈ U d such that  und 21

C

ϕ0 21

+k

n

$ f (ti )20

, 1  n  N.

(3.2.41)

i=1

√ Further, if k = O(h), L  N (for example, L2 = O(N)), and L < 5, we have the following error estimates with respect to the norm in H01 ():

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

% & l &

& n n uh − ud 1  M(n) + C 'L λj , n = 1, 2, · · · , N,

195

(3.2.42)

j =d+1

√ where M(n) = 0 (n = 1, 2, · · · , L), while M(n) = C(k + h) n − L (n = L + 1, L + 2, · · · , N). Proof. Since the left-hand side of (3.2.37) for Problem 3.2.6 is also a symmetric positive definite and bounded bilinear functional on U d , it has a unique set of d solutions {und }N n=1 ⊂ U by the Lax–Milgram theorem, i.e., Theorem 2.1.12. For n = 1, 2, · · · , L, by using (3.2.32) and (3.2.15), we have ∇und 20 = ∇P d uh 20  ∇unh 20 ! " n

2 2 2 2 2  C h ϕ0 1 + ϕ0 0 + k f (ti )0 .

(3.2.43)

i=1

in (3.2.37) of Problem 3.2.6 For n = L + 1, L + 2, · · · , N , taking vd = un+1 d and using Lemmas 3.1.3 and 3.1.4 and the Hölder and Cauchy–Schwarz inequalities, we have ||| und |||20 +ε∇und 20 + γ k∇und 20 n ||| un−1 |||0 ||| und |||0 +kf (tn )0 ||| und |||0 +ε∇un−1 d d 0 ∇ud 0 1 2 n 2  [||| un−1 |||20 + ||| und |||20 +ε∇un−1 d d 0 + ε∇ud 0 2 + kf (tn )20 + k ||| und |||20 ], (3.2.44)

i.e., ||| und |||20 +ε∇und 20 + 2γ k∇und 20 2  kf (tn )20 + k ||| und |||20 + ||| un−1 |||20 +ε∇un−1 d d 0 . (3.2.45)

Summing (3.2.45) from L + 1 to n, we have ||| und |||20 +ε∇und 20 n

2 L 2 k [f (ti )20 + ||| uid |||20 ]+ ||| uL d |||0 +ε∇ud 0 .

(3.2.46)

i=L+1

If k is sufficiently small such that k  1/2, from (3.2.46), we obtain ||| und |||20 +ε∇und 20 ||| und |||20 +2ε∇und 20 2 L 2  2 ||| uL d |||0 +2ε∇ud 0 + 2k

n

i=L+1

f (ti )20 + 2k

n−1

i=L+1

(3.2.47) ||| uid |||20 .

196 Proper Orthogonal Decomposition Methods for Partial Differential Equations

Applying the discrete Gronwall lemma, i.e., Lemma 1.4.1 to (3.2.47) and using Lemma 3.1.4 as well as combining with (3.2.43), we get (3.2.41) when max{k, h}  1. Since U d ⊂ Uh , subtracting Eqs. of Problem 3.2.6 from those of Problem 3.2.5 by taking vh = vd ∈ U d , we have unh − und = unh − P d unh , n = 1, 2, · · · , L,

(3.2.48)

(unh − und , ∗h vd ) + εa(unh − und , ∗h vd ) + γ ka(unh − und , ∗h vd ) n−1 ∗ ∗ d − un−1 − un−1 = (un−1 h d , h vd ) + εa(uh d , h vd ), ∀vd ∈ U ,

n = L + 1, L + 2, · · · , N.

(3.2.49)

For n = 1, 2, · · · , L, by using Lemma 3.2.4, from (3.2.48), we obtain % & l &

& unh − und 1  'L λj , n = 1, 2, · · · , L.

(3.2.50)

j =d+1

For n = L + 1, L + 2, · · · , N , from (3.2.31), (3.2.49), and Lemmas 3.1.3 and 3.1.4, we have ||| unh − und |||20 +ε∇(unh − und )20 + γ k∇(unh − und )20 = (unh − und , ∗h (unh − P d unh )) + (unh − und , ∗h (P d unh − und )) + εa(unh − und , unh − und ) + γ ka(unh − und , unh − und ) = (unh − unh , ∗h (unh − P d unh )) + (unh − und , ∗h (P d unh − und )) + εa(unh − und , unh − P d unh ) + εa(unh − und , P d unh − und ) + γ ka(unh − und , unh − P d unh ) + γ ka(unh − und , P d unh − und ) = (unh − unh , ∗h (unh − P d unh )) + εa(unh − P d unh , unh − P d unh ) d n n − un−1 + γ ka(unh − P d unh , unh − P d unh ) + εa(un−1 h d , P uh − ud ) ∗ d n n − un−1 + (un−1 h d , h (P uh − ud )).

(3.2.51)

By using (3.2.33) and the Hölder and Cauchy–Schwarz inequalities, we have (unh − unh , ∗h (unh − P d unh )) ||| unh − und |||0 ||| unh − P d unh |||0  Ch ||| unh − und |||0 ∇(unh − P d unh )0  C[h2 ||| unh − und |||20 +∇(unh − P d unh )20 ], ∗ d n n (un−1 − un−1 h d , h (P uh − ud ))

− un−1 |||0 ||| unh − P d unh |||0 ||| un−1 h d + ||| un−1 − un−1 |||0 ||| unh − und |||0 h d

(3.2.52)

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

197

 Ch ||| un−1 − un−1 |||0 ||| ∇(unh − P d unh ) |||0 h d + ||| un−1 − un−1 |||0 ||| unh − und |||0 h d  C[h2 ||| un−1 − un−1 |||20 +∇(unh − P d unh )20 ] h d 1 − un−1 |||20 + ||| unh − und |||20 ]. + [||| un−1 h d 2

(3.2.53)

By using (3.2.31) and the Hölder and Cauchy–Schwarz inequalities again, we have d n n a(un−1 − un−1 h d , P uh − ud ) n−1 d n n n n = a(un−1 − un−1 − un−1 h d , P uh − uh ) + a(uh d , uh − ud ) n−1 d n n n n = a(un−1 − P d un−1 − un−1 h h , P uh − uh ) + a(uh d , uh − ud ) 2 n d n 2 − P d un−1  C(∇(un−1 h h )0 + ∇(uh − P uh )0 ) 1 2 n n 2 − un−1 + [∇(un−1 h d )0 + ∇(uh − ud )0 ]. 2

(3.2.54)

Combining (3.2.52)–(3.2.54) with (3.2.51) yields ||| unh − und |||20 +ε∇(unh − und )20 + γ k∇(unh − und )20 2 n d n 2  C(∇(un−1 − P d un−1 h h )0 + ∇(uh − P uh )0 ) 1 − un−1 |||20 + ||| unh − und |||20 ] + [||| un−1 h d 2 ε 2 n n 2 + [∇(un−1 − un−1 h d ) |||0 +∇(uh − ud )0 ] 2 + Ch2 [||| un−1 − un−1 |||20 + ||| unh − und |||20 ]. h d

(3.2.55)

Further, we get ||| unh − und |||20 +ε∇(unh − und )20  Ch2 [||| unh − und |||20 + ||| un−1 − un−1 |||20 ] h d 2 n d n 2 − P d un−1 + C(∇(un−1 h h )0 + ∇(uh − P uh )0 ) 2 + ||| un−1 − un−1 |||20 +ε∇(un−1 − un−1 h d h d )0 .

(3.2.56)

Summing (3.2.56) from L + 1 to n yields ||| unh − und |||20 +ε∇(unh − und )20 ||| unh − und |||20 +2ε∇(unh − und )20 n n



 Ch2 ||| uih − uid |||20 +C ∇(uih − P d uih )20 i=L L + ||| uL h − ud

i=L

|||20

+ε∇(uL h

2 − uL d ) |||0 .

(3.2.57)

198 Proper Orthogonal Decomposition Methods for Partial Differential Equations

If k = O(h) is sufficiently small such that Ch2  1/2, from the above inequality and (3.2.50), we obtain ||| unh − und |||20 +ε∇(unh − und )20  Ck 2

n−1

||| uih − uid |||20

i=L

+C

n

∇(uih − P d uih )20 + CL

l

λj .

(3.2.58)

j =d+1

i=L

Applying Lemma 1.4.1 to (3.2.58) and using Lemmas 3.1.3 and 3.1.4 and Proposition 3.2.3, we obtain ⎡ ⎤ n l



∇(unh − und )20  C ⎣ ∇(uih − P d uih )20 + L λj ⎦ exp(Ck 2 n). j =d+1

i=L+1

(3.2.59) By Lemma 3.2.4 and Theorem 3.2.2, we have ∇(P d uih − uih )0  ∇(uih − u(ti ))0 + ∇(u(ti ) − P h u(ti ))0 + ∇P h (u(ti ) − uih )0  C[∇(u(ti ) − P h u(ti ))0 + ∇(u(ti ) − uih )0 ]  C(h + k). (3.2.60) Since n  L  N = O(k −1 ) = O(h−1 ), the factor exp(Ck 2 n) is a bounded small positive number. Thus, from (3.2.50) and (3.2.59)–(3.2.60), we obtain (3.2.42). This completes the proof of Theorem 3.2.5. Combining Theorems 3.2.2 with 3.2.5 yields the following result. Theorem 3.2.6. Under the same assumptions of Theorem 3.2.5, we have the error estimates between the solutions for Problem 3.2.2 and the solutions for the PODROEFVE algorithm Problem 3.2.6 with respect to the norm in H01 () as follows: % & l &

& n u(tn ) − ud 1  M(n) + C(k + h) + C 'L λj , n = 1, 2, · · · , N, j =d+1

√ where M(n) = 0 (n = 1, 2, · · · , L), whereas M(n) = C(k + h) n − L (n = L + 1, L + 2, · · · , N). Remark 3.2.4. The following are additional remarks on Theorems 3.2.5 and 3.2.6:

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

199

1. The estimate (3.2.41) in Theorem 3.2.5 shows that the solutions und (n = 1, 2, · · · , N ) for Problem 3.2.6 are numerically stable and continuously dependent on the source term f (x, y, t) and the initial data ϕ0 (x, y) (even when ϕ0 (x, y) is a nonzero function). √ 2. The conditions L  N and L < 5 in Theorem 3.2.5 show the relationship between the number L of snapshots and the number N of all time instants. Therefore, it is unnecessary to take all total transient solutions at every time instant tn as snapshots such as in references [62,63]. 3. The errors in Theorem 3.2.5 suggest how to choose the number of POD  # 1/2 bases: the number d of POD bases should satisfy L lj =d+1 λj  max{k, h}. √ 4. The errors (k + h) n − L (L + 1  n  N ) in Theorems 3.2.5 and 3.2.6 are errors caused by extrapolation. Thus, it serves as a guide for duly renewing the POD basis generated# by the newest solutions of Problem 3.2.6, namely, if (k 2 + h2 )(n − L)  L lj =d+1 λj (L + 1  n  N ), und (L + 1  n  N ) is the sequence of solutions for the PODROEFVE algorithm, #Problem 3.2.6, satisfying a desired accuracy; else, if (k 2 + h2 )(n − L) > L lj =d+1 λj (L + 1  n  N ), it is necessary to renew the POD basis, which is to be generated by some new snapshots Wi = un−i (i = 1, 2, · · · , L). d

3.2.4 The Implementation of the Reduced-Order Algorithm for the 2D Sobolev Equation In what follows, we provide the implementation for solving the PODROEFVE algorithm, i.e., Problem 3.2.7, which consists of seven steps. Step 1. Extraction of snapshots For given ε and γ , the source term f (x, y, t), the initial value function k = O(h), and ϕ0 (x, y), the triangulation parameter h, the time step increment √ the trial function space Uh , solving Problem 3.2.5 at the first L ( L < 5) steps produces an ensemble of snapshots Wi (x, y) = uih (i = 1, 2, · · · , L), which may alternatively be produced by drawing samples from experiments in the physical system’s evolution. Step 2. Formulation of correlation matrix A The correlation matrix A = (Aij )L×L , where Aij = (Wi , Wj )U /L = j (∇uih , ∇uh )/L and (·, ·) is L2 -inner product. Step 3. Computing eigenvalues and eigenvectors of A Find the eigenvalues λ1  λ2  · · ·  λl > 0 (l = dim{W1 , W2 , · · · , WL }) j j j and corresponding eigenvectors v j = (a1 , a2 , · · · , aL )T (j = 1, 2, · · · , l) of the matrix A. Step 4. Determine the number of POD bases ( # The number d of POD bases should satisfy L lj =d+1 λj  max{k, h}.

200 Proper Orthogonal Decomposition Methods for Partial Differential Equations

Step 5. Formulation of POD basis # j i  The POD basis ψj = L i=1 ai uh / Lλj (j = 1, 2, · · · , d). Step 6. Compute the reduced-order FVE solutions Solve the following system of equations which contains only d degrees of freedom at each time level: αjn = (∇unh , ∇ψj ), j = 1, 2, · · · , d, n = 1, 2, · · · , L,   

∂ψj (x, y) ∂ψj (x, y) Cij = ψi (xz , yz ) dx − dy , ∂y ∂x ∂Vz Vz ∈ ∗h 

ψi (xz , yz ) f (x, y, tn )dxdy, Fin = Vz

Vz ∈ ∗h

Bij =



Vz ∈ ∗h d



ψj (x, y)dxdy, i, j = 1, 2, · · · , d,

ψi (xz , yz ) Vz

αjn [Bij + (ε + γ k)Cij ] = Fin +

j =1

d

αjn−1 (Bij + εCij ),

j =1

i = 1, 2, · · · , d, n = L + 1, L + 2, · · · , N. One obtains (α1n , α2n , · · · , αdn )T ∈ R d (n = 1, 2, · · · , N ) and further the reduced# order FVE solutions und = di=1 αin ψi (x, y) (n = 1, 2, · · · , L, L + 1, · · · , N ). Step 7. Renewal of POD basis and circulation or end # If (k 2 + h2 )(n − L)  L lj =d+1 λj (L + 1  n  N ), then und (n = 1, 2, · · · , N) are the solutions for Problem 3.2.7 satisfying one desired accu# racy requirement. Else, namely, if (k 2 + h2 )(n − L) > L lj =d+1 λj (L + 1  n  N ), let Wi = uid (i = n − L, n − L − 1, · · · , n − 1) and return to Step 2.

3.2.5 Numerical Experiments In this subsection, we present some numerical experiments for the 2D Sobolev equation to validate the feasibility and efficiency for the PODROEFVE algorithm. We have performed the calculations using Matlab 7.11 (R2010b) on a computer with an Intel Core i5-3210M CPU of 2.50 GHz and Windows 7. Moreover, by comparing the numerical results of the reduced-order FVE extrapolation algorithm, the classical FVE formulation, the reduced-order FE formulation, the classical FE formulation, the reduced-order FD scheme, and the classical FD scheme for the 2D Sobolev equation, it is shown that the PODROEFVE algorithm is highly effective.

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

201

FIGURE 3.2.1 The left graphic is a reduced-order FVE solution at t = 0.4, while the right one is a classical FVE solution at t = 0.4.

3.2.5.1 Comparison: the POD-Based Reduced-Order Extrapolation Finite Volume Element Solutions Versus the Classical Finite Volume Element Solutions For the sake of convenience and without loss of generality, here we take a simple 2D Sobolev equation as an example. The computational domain consists of a square with dimensions 2 cm × 2 cm, i.e.,  = [0, 2] × [0, 2], ε = 1/(20π 2 ), γ = 1/(2π 2 ), the initial function ϕ0 (x, y) = sin πx sin πy, and the source term f (x, y, t) = 12 sin πx sin πy exp(10t). Though the data set is rather simple, the ideas and approaches here can be easily applied to much more complex situations of numerical computations. We first divide  into small squares with side length x = y = 0.001 cm. Then we link diagonals of squares so as to divide each square into two triangles in direction. Finally, we form a triangulation h with diameter h = √ the same 2 × 10−3 and nodes Nh = 4 × 106 . We take the dual decomposition ∗h as a barycentric dual decomposition, namely, the barycenter of the right triangle K ∈ h is the node of the dual decomposition. In order to satisfy the condition k = O(h) in Theorem 3.2.5, we take the time step increment as k = 10−3 . We choose the first L = 20 solutions unh (n = 1, 2, · · · , 20, i.e., at time t = 0.001, 0.002, · · · , 0.02S) as the snapshots from 400 solutions unh (n = 1, 2, · · · , 400, i.e., at time t = 0.001, 0.002, · · · , 0.4) which are obtained by solving the classical FVE formulation or Problem 3.2.5. This is(achieved by # computing that, when d = 6 and k = 10−3 , the error satisfies 20 20 j =7 λj −3  5 × 10 in Theorem 3.2.5. Thus, it is necessary to take the six major POD basis elements to expand into the subspace U d . When we find the solution at t = 0.4 from the PODROEFVE algorithm Problem 3.2.7 with six optimal POD bases, according to the seven steps of implementation of the PODROEFVE algorithm in Section 3.2.4, it is necessary to renew the POD basis once at t = 0.3. The solution at t = 0.4 obtained from Problem 3.2.6 is depicted on the left part of Fig. 3.2.1. The solution at time t = 0.4 obtained from the classical FVE formulation Problem 3.2.5 is depicted on the right part of Fig. 3.2.1. The two graphics in Fig. 3.2.1 are exhibiting close or “quasiidentical” similarity. How-

202 Proper Orthogonal Decomposition Methods for Partial Differential Equations

ever, due to Problem 3.2.6, the truncated error accumulation is reduced in the computational process. Therefore, the PODROEFVE algorithm solution could be superior to the classical FVE solution.

3.2.5.2 Comparison: the Finite Volume Element Solutions Versus the Finite Element Solutions and the Finite Difference Solutions Using the same time step increment and spatial step discretization size, data, and initial and boundary conditions as in Section 3.2.5.1, we obtain the solution by the following classical FE formulation, Problem 3.2.8, with FE space Uh of piecewise linear polynomials at time t = 0.4. Problem 3.2.8. Find unh ∈ Uh (n = 1, 2, · · · , N) such that ⎧ n n n ⎪ ⎨ (uh , vh ) + ε(∇uh , ∇vh ) + kγ (∇uh , vh ) n−1 = k(f n , vh ) + (un−1 h , vh ) + ε(∇uh , ∇vh ), ∀vh ∈ Uh , ⎪ ⎩ 0 uh = sin πx sin πy. And we also obtain the solution for which is found by the following reducedorder FE extrapolation algorithm, Problem 3.2.9, based on the POD method with six POD bases at time t = 0.4. Problem 3.2.9. Find und ∈ U d (n = 1, 2, · · · , N) such that ⎧ d ⎪

⎪ ⎪ n d n ⎪ = P u = (∇unh , ∇ψj )ψj , n = 1, 2, · · · , L, u ⎪ h ⎨ d j =1

⎪ ⎪ (und , vd ) + ε(∇und , ∇vd ) + kγ (∇und , vd ) = k(f n , vd ) + (un−1 ⎪ ⎪ d , vd ) ⎪ ⎩ n−1 d + ε(∇ud , ∇vd ), ∀vd ∈ U , n = L + 1, L + 2, · · · , N. For Problem 3.2.9, it is also necessary to renew the POD basis once at t = 0.3 with the similar seven steps of implementation of algorithm as given in Section 3.2.4, but the snapshots of Step 1 are substituted with solutions of Problem 3.2.8 and Step 6 is substituted by Problem 3.2.9. The solutions are depicted on the left and right columns in Fig. 3.2.2, respectively. Though both graphics in Fig. 3.2.2 are also manifesting quasiidentical similarity, the reduced-order FE solution (see the left graphic in Fig. 3.2.2) is comparable to or even better than the classical FE solution (see the right graphic in Fig. 3.2.2). Both graphics in Fig. 3.2.2 have known sharper dissipation than those in Fig. 3.2.1, but they have not performed better than the reduced-order FVE solution and the classical FVE solution, respectively. Also using the same time step size and spatial discretization size, data, and initial and boundary conditions as in Section 3.2.5.1, we obtain the solution at time t = 0.4 by the following classical implicit FD scheme:

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

203

FIGURE 3.2.2 The left graphic is a reduced-order FE solution at t = 0.4, while the right one is a classical FE solution when t = 0.4.

FIGURE 3.2.3 The left graphic is a reduced-order FD solution at t = 0.4, while the right one is a classical FD solution at t = 0.4.

Problem 3.2.10. Find uni,j (n = 1, 2, · · · , N) such that uni,j

+ ε(

+ νk(

uni+1,j − 2uni,j + uni−1,j uni+1,j

un−1 i+1,j

x 2 − 2uni,j + uni−1,j

x 2 n−1 − 2ui,j + un−1 i−1,j

+ +

uni,j +1 − 2uni,j + uni,j −1 y 2 uni,j +1 − 2uni,j + uni,j −1

y 2 n−1 n−1 ui,j +1 − 2ui,j + un−1 i,j −1

+ x 2 y 2 i = j = 1, 2, · · · , 2000, n = 1, 2, · · · , N, = ε(

) )

n ) + kfi,j + un−1 i,j ,

and the reduced-order FD solution at time t = 0.4 by the reduced-order FD scheme (which is derived similarly to that in [91]) based on the POD method with six POD bases. Here, it is also necessary to renew the POD basis once at t = 0.3 along the similar seven steps of implementation of algorithm in Section 3.2.4, and are depicted on the left and right columns in Fig. 3.2.3, respectively. However, the left chart and right graphic in Fig. 3.2.3 are also exhibiting close similarity and the reduced-order FD solution (see the left graphic in Fig. 3.2.3) is better than the classical FD solution (see the right graphic in Fig. 3.2.3), as they have more evident dissipation than the graphics in Figs. 3.2.1

204 Proper Orthogonal Decomposition Methods for Partial Differential Equations

FIGURE 3.2.4 The absolute error (log 10) between the classical solutions and the reduced-order solutions with different numbers of POD bases for FVE, FE, and FD formulations at t = 0.4.

TABLE 3.2.1 Comparison of CPU Time at t = 0.4 s on a ThinkPad E530 Laptop Numerical method

CPU time

PODROEFVE algorithm

4s

Classical FVE formulation

728 s

Reduced-order FE extrapolation algorithm

6s

Classical FE formulation

972 s

Reduced-order FD scheme

4.3 s

Classical FD scheme

735 s

and 3.2.2. As solutions, their accuracy is not better than the reduced-order FVE solution nor the classical FVE solution. Fig. 3.2.4 shows the absolute errors (log 10) at time t = 0.4 between the reduced-order FVE solutions with different number of POD bases versus the classical FVE solution, the reduced-order FE solutions with different number of POD bases versus the classical FE solution, and the reduced-order FD solutions with different number of POD bases versus the classical FD solution. It again shows that the numerical results of FVE formulations are the best among the three numerical methods. Table 3.2.1 shows the comparison of CPU time for the implementation of numerical simulation at t = 0.4 for the PODROEFVE algorithm, Problem 3.2.6, containing six POD bases, the classical FVE formulation, Problem 3.2.5, the reduced-order FE extrapolation algorithm, Problem 3.2.9, based on the POD method containing six POD bases, the classical FE formulation, Problem 3.2.8, the reduced-order FD scheme based on POD method containing six POD bases, and the classical FD scheme, Problem 3.2.10. It is shown that the PODROEFVE method is most efficient as far as the CPU time is concerned. Fig. 3.2.4 and Table 3.2.1 have shown that the PODROEFVE algorithm can reduce the degrees of freedom, alleviate the computational load, save CPU time of calculations, and lessen the truncation error accumulation in the com-

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

205

putational process. Therefore, the PODROEFVE algorithm is one of the most effective numerical methods, and the results for numerical experiments are consistent with theoretical estimates since theoretical and numerical errors do not exceed 5 × 10−3 . It is also shown that finding the numerical solutions for the 2D Sobolev equation with the PODROEFVE algorithm, Problem 3.2.6, is computationally most effective.

3.3 POD-BASED REDUCED-ORDER STABILIZED CRANK–NICOLSON EXTRAPOLATION MIXED FINITE VOLUME ELEMENT MODEL FOR THE TWO-DIMENSIONAL NONSTATIONARY INCOMPRESSIBLE BOUSSINESQ EQUATION In this section, we first establish a semidiscretized CN formulation with respect to time for the 2D nonstationary incompressible Boussinesq equation. We formulate a fully discretized stabilized CN mixed FVE (SCNMFVE) scheme based on two local Gaussian quadratures and parameter-free real. This is done by building directly from the semidiscretized CN formulation with respect to time. We then derive the error estimates for the fully discretized SCNMFVE solutions by means of the standard CN MFE (CNMFE) method. This work was published in Luo [86]. Continuing the above, we develop a POD-based reduced-order stabilized CN extrapolation mixed FVE (PODROSCNEMFVE) model with very few degrees of freedom for the 2D nonstationary incompressible Boussinesq equation and derive error estimates of the PODROSCNEMFVE solutions and implement the algorithm of the model. Finally, we present numerical experiments to validate the reliability of the PODROSCNEMFVE model. Moreover, it is shown that the numerical errors are consistent with theoretical analysis and the computing load for the PODROSCNEMFVE model is far lower than that for the classical SCNMFVE formulation. The TE accumulation in the computational process is also far less than that of the classical SCNMFVE formulation. Thus, the advantages of the model for the 2D nonstationary incompressible Boussinesq equation are excellent. The latter part work was published in Luo [87].

3.3.1 Model Background and Survey for the 2D Nonstationary Incompressible Boussinesq Equation The nonstationary incompressible Boussinesq equation is a nonlinear system of PDEs including the velocity vector field and the pressure field as well as the temperature distribution (see [80,105,178]). This model constitutes a nonstationary conduction–convection problem and may be described by the following nonlinear system of PDEs.

206 Proper Orthogonal Decomposition Methods for Partial Differential Equations

Problem 3.3.1. Find u = (u1 , u2 )T , p, and T such that, for tN > 0, ⎧ (x, y, t) ∈  × (0, tN ), ut − νu + (u · ∇)u + ∇p = T j , ⎪ ⎪ ⎪ ⎪ ⎪ ∇ · u = 0, (x, y, t) ∈  × (0, tN ), ⎪ ⎨ −1 (x, y, t) ∈  × (0, tN ), Tt − γ0 T + (u · ∇)T = 0, ⎪ ⎪ ⎪ u(x, y, t) = u0 (x, y, t), T (x, y, t) = ϕ(x, y, t), (x, y, t) ∈ ∂ × (0, tN ), ⎪ ⎪ ⎪ ⎩ (x, y) ∈ , u(x, y, 0) = u0 (x, y), T (x, y, 0) = ψ(x, y), (3.3.1) where  ⊂ R 2 is a bounded, open, and connected domain and u = (u1 , u2 )T represents the fluid velocity vector, p the pressure, T the temperature, tN the terminal total time, and j = (0, 1)T the unit vector in the y-direction, √ √ ν = P r/Re, Re the Reynolds number, P r the Prandtl number, γ0 = ReP r, and u0 (x, y, t), u0 (x, y), ϕ(x, y, t) and ψ(x, y) are given functions. For the sake of convenience and without loss of generality, we may assume that u0 (x, y, t) = u0 (x, y) = 0 and ϕ(x, y, t) = 0 in the following analysis. In general, no analytical solutions are available for Problem 3.3.1 due to the complexity of the nonlinear PDEs. One must rely on numerical solutions (see, e.g., [33,80,105]). Most of the existing literature uses either the FE method or FD scheme as discretization tools. Just as mentioned in Sections 3.1.1 or 3.2.1, compared to FD and FE methods, the FVE method (see [23,24,61,153]) is considered a highly effective discretization approach for PDEs, since it is generally easier to implement and offer flexibility in the handling of complicated domains. An additional advantage is that it can ensure local mass conservation, which is a highly desirable feature. The FVE method has been found particularly effective for fluid dynamics problems (see, e.g., [4,32,49,64,73,86,146,187,188]). A fully discretized FVE formulation without adopting any stabilization (see [72]) and a fully discretized stabilized mixed FVE (SMFVE) formulation (see [105]) for the 2D nonstationary incompressible Boussinesq equation have been proposed, but they have only the first-order time accuracy. Thus, in order to obtain sufficient accuracy with respect to the time variable, there is a need to refine the time steps so that, as we move on, the truncation errors in the computational process would quickly accumulated. Therefore, in [87], we have improved the methods in [72,105] and established a fully discretized SCNMFVE formulation based on two local Gaussian quadratures and parameter-free real for the 2D nonstationary incompressible Boussinesq equation. Thus, although the trial function spaces of the fluid velocity, temperature, and pressure of the SCNMFVE formulation are the same as those in [72,105], the SCNMFVE solution has improved the first-order accuracy in time more than those in [72,105] such that it could alleviate the computational load and decrease the accumulation of truncation errors in the process. To this end, in this section, we first introduce

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

207

in detail the SCNMFVE formulation for the 2D nonstationary incompressible Boussinesq equation. As mentioned in Sections 3.1.1 and 3.2.1, even if the classical SCNMFVE formulation for the nonstationary incompressible Boussinesq equation in [105] is more advantageous than its fully discretized FVE method without adopting any stabilization in [72] (for example, it can avoid the constraint of the Babuška– Brezzi (B-B) condition and its numerical solutions are more stable than those in [72]), it contains large degrees of freedom (namely, unknowns, which are the same as those of the fully discretized FVE method in [72]) in such a way that it generates large computational loads for the practical engineering problems. Therefore, our second important task is to employ the POD technique (see [56,60]) to reduce the degrees of freedom for the fully discretized SCNMFVE formulation in order to alleviate the computational load and decrease the TE accumulation in the computational process, which lead to speedup in an effective way that guarantees sufficiently accurate numerical solutions. As done by us before, the POD method has been used to establish some POD-based reduced-order Galerkin, FE, and FD numerical models for the timedependent PDEs (see [14,28,63,93,118,121,135,137,177]); some reduced-order FD schemes (see [81,156]) and mixed FE formulations (see [89,119]) based on the POD method had been studied for the nonstationary incompressible Boussinesq equation, even though a POD-based reduced-order FVE algorithm for viscoelastic equations (see [71]), a POD-based reduced-order CN FVE formulation for parabolic equations (see [104]), and a reduced-order extrapolation algorithm based on the SFVE method and POD technique for the nonstationary Stokes equation (see [84]) have been presented, because the nonstationary incompressible Boussinesq equation is more complicated than the above-mentioned models, the PODROSCNEMFVE model for the system of nonstationary incompressible Boussinesq equation is far more challenging. Indeed, most of the existing POD-based reduced-order models (see, e.g., [14,28,63,71,89,93, 104,118,119,121,135,137,156,177]) have used the numerical solutions obtained from the classical numerical models on the total time span [0, tN ] to formulate the POD bases and, from them, build the POD-based reduced-order models. Solutions must then be recomputed on the same time span [0, tN ]. Our PODROSCNEMFVE model improves the existing POD-based reduced-order models, by using only the first few SCNMFVE solutions on the very short time span [0, t˜] (t˜  tN ) as snapshots to formulate the POD basis and build the PODROSCNEMFVE model for the nonstationary incompressible Boussinesq equation. We then find the numerical solutions on the time span [t˜, tN ] by extrapolating iteration. Thus, our PODROSCNEMFVE model has the synthesized advantages of the POD technique and the FVE method.

208 Proper Orthogonal Decomposition Methods for Partial Differential Equations

3.3.2 Semidiscretized Crank–Nicolson Formulation About Time for the 2D Nonstationary Incompressible Boussinesq Equation The Sobolev spaces along with their properties used in this context are standard (see Section 2.1.1 or [1]).  the notation ) Here we adapt * as follows. For brevity, we write M = L20 () = q ∈ L2 ();  qdxdy = 0 , X = H01 ()2 , and W = H01 (). A mixed variational formulation for Problem 3.3.1 can be stated as follows. Problem 3.3.2. Find (u, p, T ) ∈ H 1 (0, tN ; X)2 ×L2 (0, tN ; M)×H 1 (0, tN ; W ) such that, for almost all t ∈ (0, tN ), ⎧ (ut , v) + a(u, v) + a1 (u, u, v) − b(p, v) = (T j , v), ∀v ∈ X, ⎪ ⎪ ⎪ ⎨b(q, u) = 0, ∀q ∈ M, (3.3.2) ⎪ , φ) + d(T , φ) + a (u, T , φ) = 0, ∀φ ∈ W, (T t 2 ⎪ ⎪ ⎩ u(x, y, 0) = 0, T (x, y, 0) = ψ(x, y), (x, y) ∈ , where (·, ·) denotes the inner product in L2 ()2 or L2 (), and  a(u, v) = ν ∇u · ∇vdxdy, ∀u, v ∈ X,   b(q, v) = q divvdxdy, ∀v ∈ X, q ∈ M,   1 a1 (u, v, w) = [(u∇v) · w − (u∇w) · v] dxdy, ∀u, v, w ∈ X, 2   1 a2 (u, T , φ) = [(u · ∇T )φ − (u · ∇φ)T ] dxdy, ∀u ∈ X, ∀T , φ ∈ W, 2  d(T , φ) = γ0−1 ∇T · ∇φdxdy, ∀T , φ ∈ W. 

The trilinear forms above a1 (·, ·, ·) and a2 (·, ·, ·) have the following properties (see [72,80]): a1 (u, v, w) = −a1 (u, w, v), a1 (u, v, v) = 0, ∀u, v, w ∈ X,

(3.3.3)

a2 (u, T , φ) = −a2 (u, φ, T ), a2 (u, φ, φ) = 0, ∀u ∈ X, ∀T , φ ∈ W. (3.3.4) The bilinear forms above a(·, ·), d(·, ·), and b(·, ·) have the following properties (see also [72,80]): a(v, v)  ν | v |21 , ∀v ∈ X, | a(u, v) | ν | u |1 | v |1 , ∀u, v ∈ X, d(φ, φ)  γ0−1

|φ

|21 ,

∀φ ∈ W,

| d(T , φ) | γ0−1

(3.3.5)

| T |1 | φ |1 , ∀T , φ ∈ W, (3.3.6)

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

b(q, v)  βq0 , ∀q ∈ M, v∈X | v |1 sup

209

(3.3.7)

where β is a constant. Define a1 (u, v, w) , u,v,w∈X | u |1 · | v |1 · | w |1 a2 (u, T , φ) sup . N˜ 0 = u∈X,(T ,φ)∈W ×W | u |1 · | T |1 · | φ |1

N0 =

sup

(3.3.8)

The following result is classical (see [178, Theorem 1.4.1] or [80, Theorem 5.2]). Theorem 3.3.1. If ψ ∈ L2 (), then Problem 3.3.2 has at least a solution. In addition, the solution is unique provided that ψ20  2ν 2 tN /(2N0 tN−1 exp(tN )+ νγ0 N˜ 02 ). The following a priori estimate holds: u20 + ν∇u2L2 (L2 )  tN2 ψ20 exp(tN ), T 20 + γ0−1 ∇T 2L2 (L2 )  ψ20 . Let N be a positive integer, let k = tN /N denote the time step increment, let tn = nk (0  n  N ), let (un , p n , T n ) be the semidiscretized approximation of (u(t), p, T ) at tn = nk (n = 0, 1, · · · , N ) with respect to time, and let u¯ n = (un + un−1 )/2 and T¯ n = (T n + T n−1 )/2. If the differential quotients ut and Tt in Problem 3.3.2 at time t = tn are, respectively, approximated by means of the backward difference quotients ∂¯t un = (un − un−1 )/k and ∂¯t T n = (T n − T n−1 )/k, then the semidiscretized CN scheme for Problem 3.3.2 with respect to time can be stated as follows. Problem 3.3.3. Find (un , p n , T n ) ∈ X × M × W (1  n  N ) such that ⎧ n (u , v) + ka(u¯ n , v) + ka1 (u¯ n , u¯ n , v) − kb(p n , v) ⎪ ⎪ ⎪ ⎪ n n−1 ⎪ ⎪ ⎨ = k(T¯ j , v) + (u , v), ∀v ∈ X, b(q, un ) = 0, ∀q ∈ M, (3.3.9) ⎪ ⎪ n n n n n−1 ⎪ ¯ ¯ ⎪ (T , φ) + kd(T , φ) + ka2 (u¯ , T , φ) = (T , φ), ∀φ ∈ W, ⎪ ⎪ ⎩ 0 0 u = 0, T = ψ(x, y), (x, y) ∈ . We have the following for Problem 3.3.3. Theorem 3.3.2. Under the assumptions of Theorem 3.3.1, Problem 3.3.3 has a unique sequence of solutions (un , p n , T n ) ∈ X × M × W (n = 1, 2, · · · , N ) such that ( √ (3.3.10) un 20 + k ν∇un 0  ν −1 γ0 ψ0 ,

210 Proper Orthogonal Decomposition Methods for Partial Differential Equations −1/2

T n 0 + kγ0

√ −1/2 ∇T n 0  max{2, 2kγ0 }ψ1 , kp n 0  Cψ1 , (3.3.11)

where C is a generic constant independent of k, but dependent on ψ and the known data. If (u, p, T ) ∈ W 3,∞ (0, tN ; X)2 × H 1 (0, tN ; M) × W 3,∞ (0, tN ; W ) is the exact solution for the Problem 3.3.1, we have the following error estimates: u(tn ) − un 0 + k[∇(u(tn ) − un )0 + p(tn ) − p n 0 ]  Ck 2 , (3.3.12) T (tn ) − T n 0 + k∇(T (tn ) − T n )0  Ck 2 , 1  n  N.

(3.3.13)

Proof. First, we prove that Problem 3.3.3 has a unique sequence of solutions. For 1  n  N , we consider the following linearized auxiliary problem: ⎧ n , v) ⎪ (unm , v) + ka(u¯ nm , v) + ka1 (u¯ nm−1 , u¯ nm , v) − kb(pm ⎪ ⎪ ⎪ ⎪ ⎪ = k(T¯mn j , v) + (un−1 ⎪ m , v), ∀v ∈ X, m = 1, 2, · · · , ⎪ ⎪ ⎨ b(q, un ) = 0, ∀q ∈ M, m = 1, 2, · · · , m n ⎪ (Tm , φ) + kd(T¯mn , φ) + ka2 (u¯ nm−1 , T¯mn , φ) = (Tmn−1 , φ), ⎪ ⎪ ⎪ ⎪ ⎪ ∀φ ∈ W, m = 1, 2, · · · , ⎪ ⎪ ⎪ ⎩ 0 um (x, y) = 0, Tm0 = ψ(x, y), (x, y) ∈ , m = 1, 2, · · · .

(3.3.14)

n n n−1 in the system of Eq. By taking v = unm + un−1 m , q = pm , and φ = Tm + Tm (3.3.14) and by using (3.3.3), (3.3.4), and the Hölder and Cauchy–Schwarz inequalities, we obtain 2 n n−1 2 2(unm 20 − un−1 m 0 ) + kν∇(um + um )0 = 2k(T¯mn j , unm + un−1 m )

 kTmn + Tmn−1 −1 ∇(unm + un−1 m )0 k kν 2  Tmn + Tmn−1 2−1 + ∇(unm + un−1 m )0 2ν 2

(3.3.15)

and 2(Tmn 20 − Tmn−1 20 ) + kγ0−1 ∇(Tmn + Tmn−1 )20 = 0.

(3.3.16)

Summing (3.3.15) and (3.3.16) from 1 to n and simplifying, we get unm 20

+ kν

n

i=1

∇(uim

2 + ui−1 m )0

n k i  Tm + Tmi−1 2−1 2ν i=1

(3.3.17)

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

211

and Tmn 20 + kγ0−1

n

∇(Tmi + Tmi−1 )20  2ψ20 .

(3.3.18)

i=1

By (3.3.18), from (3.3.17), we have unm 20

+ kν

n

2 −1 2 ∇(uim + ui−1 m )0  ν γ0 ψ0 .

(3.3.19)

i=1

Extracting the # √ square root for (3.3.18) and (3.3.19) and using n | b | / n and a + b0  a0 − b0 , we obtain i i=1 −1/2

Tmn 0 + kγ0

#n

√ −1/2 ∇Tmn 0  max{2, 2kγ0 }ψ1 ,

( √ unm 20 + k ν∇unm 0  ν −1 γ0 ψ0 .

 2 1/2 i=1 bi



(3.3.20)

(3.3.21)

Using the first equation of (3.3.9) and (3.3.7), (3.3.20), and (3.3.21), we easily get n 0  Cψ1 . pm

(3.3.22)

Thus, for 1  n  N , if ψ = 0, the system of linear equations (3.3.14) has only the zero solution. Therefore, the system of linear equations (3.3.14) has a unique n , T n ) ∈ X × M × W (m = 1, 2, · · · ). Because sequence of solutions (unm , pm m the spaces X × M × W are weakly and sequentially compact Hilbert spaces, by the sequential compactness of bounded sequence (see [143]), we conclude n , T n ) ∈ X × M × W has a subsequence that the sequence of solutions (unm , pm m n , T n )) that is uniquely (without loss of generality, we still might denote (unm , pm m n n n weakly convergent to (u , p , T ) ∈ X × M × W for Problem 3.3.3, i.e., Problem 3.3.3 has at least a sequence of solutions (un , p n , T n ) ∈ X × M × W (n = 1, 2, · · · , N ). By using the same technique as in the proof of the uniqueness of solution for Problem 3.3.2 (see [178, Theorem 1.4.1] or [80, Theorem 5.2]), we can prove that the sequence of solutions for Problem 3.3.3 is unique. Second, we prove that (3.3.10) and (3.3.11) hold. By taking v = un + un−1 , q = p n , and φ = T n + T n−1 in Problem 3.3.3 and by using (3.3.3), (3.3.4), and the Hölder and the Cauchy–Schwarz inequalities, we obtain 2(un 20 − un−1 20 ) + kν∇(un + un−1 )20 = 2k(T¯ n j , un + un−1 )  kT n + T n−1 −1 ∇(un + un−1 )0 k kν  T n + T n−1 2−1 + ∇(un + un−1 )20 2ν 2

(3.3.23)

212 Proper Orthogonal Decomposition Methods for Partial Differential Equations

and 2(T n 20 − T n−1 20 ) + kγ0−1 ∇(T n + T n−1 )20 = 0.

(3.3.24)

Summing (3.3.23) and (3.3.24) from 1 to n and simplifying, we obtain un 20 + kν

n

∇(ui + ui−1 )20 

i=1

n k i T + T i−1 2−1 2ν

(3.3.25)

i=1

and T n 20 + kγ0−1

n

∇(T i + T i−1 )20  2ψ20 .

(3.3.26)

i=1

By (3.3.26), from (3.3.25), we get un 20 + kν

n

∇(ui + ui−1 )20  ν −1 γ0 ψ20 .

(3.3.27)

i=1

By extracting√the square root of (3.3.26) and (3.3.27) and using # n i=1 | bi | / n and a + b0  a0 − b0 , we obtain −1/2

T n 0 + kγ0

#n

√ −1/2 ∇T n 0  max{2, 2kγ0 }ψ1 ,

( √ un 20 + k ν∇un 0  ν −1 γ0 ψ0 .

 2 1/2 i=1 bi



(3.3.28)

(3.3.29)

By the first equation of (3.3.9) and (3.3.7), (3.3.28), and (3.3.29), we easily get pn 0  Cψ1 .

(3.3.30)

Finally, we prove that the error estimates (3.3.12) and (3.3.13) hold. Set en = u(tn ) − un , θ n = T (tn ) − T n , and ηn = p(tn ) − p n . By subtracting Eqs. of Problem 3.3.3 from those of Problem 3.3.2, by taking t = tn− 1 , v = en + en−1 , 2

ϕ = θ n + θ n−1 , and q = ηn , and using Taylor’s formula, we obtain kν ∇(en + en−1 )20 2 = k((θ n + θ n−1 )j , en + en−1 ) k3 k3ν + (uttt (ξ1n ), en + en−1 ) + (∇utt (ξ2n ), ∇(en + en−1 )) 24 4 k3 − (Ttt (ξ3n )j , en + en−1 ) + , (3.3.31) 4

en 20 − en−1 20 +

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

θ n 20 − θ n−1 20 + =

213

k ∇(θ n + θ n−1 )20 2γ0

k3 (Tttt (ζ1n ), θ n + θ n−1 ) 24 k3 + (∇Ttt (ζ2n ), ∇(θ n + θ n−1 )) + , 4γ0

(3.3.32)

where  = ka1 (u(tn− 1 ), u(tn− 1 ), en + en−1 ) − ka1 (u¯ n , u¯ n , en + en−1 ) and 2 2  = ka2 (u(t 1 ), T (t 1 ), θ n + θ n−1 ) − ka2 (u¯ n , T¯ n , θ n + θ n−1 ) (tn−1  n− 2

n− 2

ξ1n , ξ2n , ξ3n , ζ1n , ζ2n  tn ). By using Taylor’s formula and the Hölder and Cauchy–Schwarz inequalities, there are ξin ∈ [tn−1 , tn ] (i = 4, 5) such that  = ka1 (u(tn− 1 ), u(tn− 1 ), en + en−1 ) − ka1 (u¯ n , u¯ n , en + en−1 ) 2

2

= ka1 (u(tn− 1 ) − u¯ n , u(tn− 1 ), en + en−1 ) 2

2

+ ka1 (u¯ n , u(tn− 1 ) − u¯ n , en + en−1 ) 2

k3 k3 ¯ utt (ξ5n ), en + en−1 ) = a1 (utt (ξ4n ), u(tn− 1 ), en + en−1 ) + a1 (u, 2 16 16 k 5 N02 ∇u(t)2W 2,∞ (t ,t ;L2 ) (∇u(t)2L2,∞ (t ,t ;L2 ) + ∇ u¯ n 20 )  n−1 n n−1 n 128ν kν + ∇(en + en−1 )20 . (3.3.33) 12 Again, by the Hölder and Cauchy–Schwarz inequalities, we have k3 k3ν (uttt (ξ1n ), en + en−1 ) + (∇utt (ξ2n ), ∇(en + en−1 )) 24 4 k3 − (Ttt (ξ3n )j , en + en−1 ) | 4 kν k5  ∇(en + en−1 )20 + u2W 3,∞ (t ,t ;H −1 ) n−1 n 12 64ν 5 9νk 5 9k + ∇u2W 2,∞ (t ,t ;L2 ) + T 2W 3,∞ (t ,t ;H −1 ) . n−1 n n−1 n 16 16ν |

(3.3.34)

Combining (3.3.33) and (3.3.34) with (3.3.31), using the Hölder and Cauchy– Schwarz inequalities, and simplifying, we get en 20 − en−1 20 + 

kν ∇(en + en−1 )20 4

3k n θ + θ n−1 2−1 + C˜ 2 k 5 , ν

(3.3.35)

214 Proper Orthogonal Decomposition Methods for Partial Differential Equations

where 1 9ν C˜ 2 = u2W 3,∞ (t ,t ;H −1 ) + ∇u2W 2,∞ (t ,t ;L2 ) n n−1 n−1 n 64ν 16 9 + T 2W 3,∞ (t ,t ;H −1 ) n−1 n 16ν 2 N + 0 ∇u(t)2W 2,∞ (t ,t ;L2 ) (∇u(t)2L2,∞ (t ,t ;L2 ) + ∇ u¯ n 20 ). n−1 n n−1 n 128ν By using the same techniques as those in the proof of (3.3.33), we can prove  = ka2 (u(tn− 1 ), T (tn− 1 ), θ n + θ n−1 ) − ka2 (u¯ n , T¯ n , θ n + θ n−1 ) 2

2

k 5 N˜ 02 γ0 k  ∇(θ n + θ n−1 )20 + ∇T (t)2W 2,∞ (t ,t ;L2 ) ∇ u¯ n 20 n−1 n 8γ0 64 k 5 N˜ 02 γ0 + ∇u(t)2W 2,∞ (t ,t ;L2 ) ∇T (t)2L2,∞ (t ,t ;L2 ) . (3.3.36) n−1 n n−1 n 64 By using the Hölder and Cauchy–Schwarz inequalities, we have |

k3 k3 (∇Ttt (ζ2n ), ∇(θ n + θ n−1 )) | (Tttt (ζ1n ), θ n + θ n−1 ) + 24 4γ0 k k 5 γ0  ∇(θ n + θ n−1 )20 + T 2W 3,∞ (t ,t ;H −1 ) n−1 n 8γ0 144 5 k + ∇T 2W 2,∞ (t ,t ;L2 ) . (3.3.37) n−1 n 4

Combining (3.3.36) and (3.3.37) with (3.3.32) and simplifying, we get θ n 20 − θ n−1 20 +

k ∇(θ n + θ n−1 )20  Cˆ 2 k 5 , 4γ0

(3.3.38)

where N˜ 02 γ0 ∇u(t)2W 2,∞ (t ,t ;L2 ) ∇T (t)2L2,∞ (t ,t ;L2 ) n−1 n n−1 n 64 1 γ 0 + T 2W 2,∞ (t ,t ;H 1 ) + T 2W 3,∞ (t ,t ;H −1 ) n−1 n n−1 n 4 144 N˜ 02 γ0 + ∇T (t)2W 2,∞ (t ,t ;L2 ) ∇ u¯ n 20 . n−1 n 64

Cˆ 2 k 5 =

Summing (3.3.38) from 1 to n yields θ n 20

n k

+ ∇(θ i + θ i−1 )20  Cˆ 2 nk 5 . 4γ0 i=1

(3.3.39)

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

215

By extracting the square root of (3.3.39) and using a + b0  a0 − b0 and  #n #n √ 2 1/2  i=0 bi i=0 | bi | / n, we obtain √ k θ n 0 + √ ∇θ n 0  Cˆ tN k 2 , 2 γ0

(3.3.40)

which yields (3.3.13). By (3.3.35) and (3.3.40), we have en 20 − en−1 20 +

kν ∇(en + en−1 )20  Cˆ 02 k 5 , 4

(3.3.41)

where Cˆ 02 = 12γ0 Cˆ 2 T∞ ν −1 k 3 + C˜ 2 . Summing (3.3.41) from 1 to n, we get en 20 +

n kν

∇(ei + ei−1 )20  Cˆ 02 nk 5 . 4

(3.3.42)

i=1

By extracting the square root of (3.3.42) and using a + b0  a0 − b0 and  #n #n √ 2 1/2  i=0 bi i=0 | bi | / n, we obtain √ √ k ν ∇en 0  Cˆ 02 tN k 2 . e 0 + 2

(3.3.43)

n

Using Taylor’s formula, we have ξin ∈ [tn−1 , tn ] (i = 5, 6, 7, 8, 9, 10, 11) satisfying b(p(tn ) − p n , v) 1 ν = (en − en−1 , v) + (∇(en + en−1 ), ∇v) k 2 1 n n−1 + [a1 (e + e , u(tn− 1 ), v) 2 2 + a1 (u¯ n , en + en−1 , v) − ((θ n + θ n−1 )j , v)] −

k2 (uttt (ξ5n ), v) 48

νk 2 νk 2 k2 (uttt (ξ6n ), v) − (∇utt (ξ7n ), ∇v) − (∇utt (ξ8n ), ∇v) 48 16 16 k2 k + b(v, pt (ξ9n )) − a1 (utt (ξ10n ), u(tn− 1 ), v) 2 2 16 2 k k2 − a1 (u¯ n , utt (ξ11n ), v) − (Ttt (ξ3n )j , v), ∀v ∈ X. (3.3.44) 16 16

−

Then, with (3.3.40), (3.3.43), (3.3.44), (3.3.7), and the Hölder and Cauchy– Schwarz inequalities, we can prove b(p(tn ) − p n , v)  Ck. ∇v0 v∈X

p(tn ) − p n 0  β −1 sup

(3.3.45)

216 Proper Orthogonal Decomposition Methods for Partial Differential Equations

Combining (3.3.45) and (3.3.43) yields (3.3.12). This completes the proof of Theorem 3.3.2.

3.3.3 Fully Discretized Stabilized Crank–Nicolson Mixed Finite Volume Element Formulation for the 2D Nonstationary Incompressible Boussinesq Equation In order to formulate the fully discretized SCNMFVE scheme for Problem 3.3.2, it is necessary to introduce an FVE approximation for the spatial variables of Problem 3.3.3 such as in Section 3.1.3 (for more details see [12,65]). Let h = {K} be a quasiuniform triangulation of  with h = max hK , where hK is the diameter of the triangle K ∈ h , and let ∗h be a dual partition based on h whose elements are called the control volumes as described in Section 3.1.3. The trial function spaces Xh , Wh , and Mh of the velocity, temperature, and pressure are, respectively, defined as follows:

Xh = v h ∈ X ∩ C()2 ; v h |K ∈ [P1 (K)]2 , ∀K ∈ h , * ) Wh = wh ∈ W ∩ C(); wh |K ∈ P1 (K), ∀K ∈ h , Mh = {qh ∈ M; qh |K ∈ P1 (K), ∀K ∈ h } ,

(3.3.46)

where P1 (K) is the space of linear functions on K. It is obvious that Xh ⊂ X = H01 ()2 and Wh ⊂ W = H01 (). For (u, T ) ∈ X × W , let (h u, ρh T ) be the interpolation projection of (u, T ) onto the trial function spaces Xh × Wh . With these in place, by the interpolation theory in Sobolev spaces (see Theorem 2.1.10 or [21,31,65,72,80]), we have the following error estimates: | u − h u |m  Ch2−m | u |2 , ∀u ∈ H 2 ()2 , m = 0, 1,

(3.3.47)

| T − ρh T |m  Ch2−m | T |2 , ∀T ∈ H 2 (), m = 0, 1,

(3.3.48)

where C is a generic constant independent of the spatial mesh size h and temporal mesh size k. The test function spaces X˜ h and W˜ h of the velocity and temperature are, respectively, chosen as follows: X˜ h = v h ∈ L2 ()2 ; v h |Vz ∈ [P0 (Vz )]2 (Vz ∩ ∂ = ∅), * v h |Vz = 0 (Vz ∩ ∂ = ∅), ∀Vz ∈ ∗h , W˜ h = wh ∈ L2 (); wh |Vz ∈ P0 (Vz ) (Vz ∩ ∂ = ∅), * wh |Vz = 0 (Vz ∩ ∂ = ∅), ∀Vz ∈ ∗h ,

(3.3.49)

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

217

where P0 (Vz ) is the space of constant functions on the control volume Vz . In fact, they can be spanned by the following basis functions:  1, (x, y) ∈ Vz , (3.3.50) φz (x, y) = ∀z ∈ Zh◦ . 0, elsewhere, For (u, w) ∈ X × W , let (∗h u, ρh∗ w) be the interpolation projection of (u, w) onto the test function spaces X˜ h × W˜ h , i.e.,

(∗h u, ρh∗ w) = (u(z), w(z))φz . (3.3.51) z∈Zh◦

By the interpolation theory in Sobolev spaces (see Theorem 2.1.10 or [21,31,65, 72,80]), we have u − ∗h u0  Ch | u |1 ,

w − ρh∗ w0  Ch | w |1 .

(3.3.52)

By using the same ideas as in references [72,105], the fully discretized SCNMFVE formulation for Problem 3.3.2 can now be stated as follows. Problem 3.3.4. Find (unh , phn , Thn ) ∈ Xh × Mh × Wh (1  n  N ) such that ⎧ ⎪ (∂¯t unh , ∗h v h ) + ah (u¯ nh , ∗h v h ) + a1h (u¯ nh , u¯ nh , ∗h v h ) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ + bh (phn , ∗h v h ) = (T¯hn j , ∗h v h ), ∀v h ∈ Xh , ⎪ ⎪ ⎨ b(qh , unh ) + Dh (phn , qh ) = 0, ∀qh ∈ Mh , (3.3.53) ⎪ ⎪ (∂¯t Thn , ρh∗ ϕh ) + dh (T¯hn , ρh∗ ϕh ) + a2h (u¯ nh , T¯hn , ρh∗ ϕh ) ⎪ ⎪ ⎪ ⎪ = 0, ∀ϕh ∈ Wh , ⎪ ⎪ ⎪ ⎩ 0 uh = 0, Th0 = ρh ψ(x, y), (x, y) ∈ , where ah (unh , ∗h v h ) = −ν bh (qh , ∗h v h ) =

 Vz ∈ ∗h



∂Vz

(v h (z)∇unh ) · nds,

(3.3.54)

 v h (z)

qh nds,

(3.3.55)

∂Vz

Vz ∈ ∗h

a1h (unh , w nh , ∗h v h ) = ((unh · ∇)w nh , ∗h v h ) + ((divunh )w nh , ∗h v h )/2, (3.3.56) 

wh (z) ∇Thn · nds, (3.3.57) dh (Thn , ρh∗ wh ) = −γ0−1 Vz ∈ ∗h

∂Vz

a2h (unh , Thn , ρh∗ ϕh ) = ((unh · ∇)Thn , ρh∗ ϕh ) + ((divunh )Thn , ρh∗ ϕhn )/2,

(3.3.58)

218 Proper Orthogonal Decomposition Methods for Partial Differential Equations

Dh (phn , qh ) = ε

+ K∈ h

K,2

,

 phn qh dxdy

− K,1

phn qh dxdy

;

 here ε is a positive real number, i.e., parameter-free real, and K,i g(x, y)dxdy indicates an appropriate Gaussian quadrature on K that is exact for polynomials of degree i (i = 1, 2) for all g(x, y) = ph qh , which is a polynomial of degree not more than i (i = 1, 2). Thus, for all test functions qh ∈ Mh , the trial function ph ∈ Mh must be piecewise constant when i = 1. Consequently, we define the L2 -projection operator h : L2 () → Wˆ h such that, for any p ∈ L2 (), (p, qh ) = (h p, qh ), ∀qh ∈ Wˆ h ,

(3.3.59)

where Wˆ h ⊂ L2 () denotes the space of piecewise constants associated with h . The projection operator h has the following properties (see Theorem 2.1.19 or [21,80]): h p0  p0 , ∀p ∈ L2 (),

(3.3.60)

p − h p0  Chp1 , ∀p ∈ H 1 ().

(3.3.61)

Now, using the definition of h , we can rewrite the bilinear form Dh (·, ·) as follows: Dh (ph , qh ) = ε(ph − h ph , qh ) = ε(ph − h ph , qh − h qh ).

(3.3.62)

3.3.4 Existence, Stability, and Convergence of the Stabilized Crank–Nicolson Mixed Finite Volume Element Solutions for the 2D Nonstationary Incompressible Boussinesq Equation In order to discuss the existence, uniqueness, stability, and error estimates of the solutions for the fully discretized SCNMFVE formulation with second-order time accuracy, i.e., Problem 3.3.4, we introduce some preliminary lemmas. First, we have the following lemma (whose proof is also available in [12,65, 73,80]). Lemma 3.3.3. The following properties hold: i. ah (uh , ∗h v h ) = a(uh , v h ), a1h (v h , uh , ∗h uh ) = 0, ∀uh , v h , wh ∈ Uh ; ii. dh (Th , ρh∗ φh ) = d(Th , φh ), a2h (uh , Th , ρh∗ Th ) = 0, ∀Th , φh ∈ Wh , ∀uh ∈ Xh ; iii. bh (ph , ∗h v h ) = −b(ph , v h ), ∀v h ∈ Xh , ∀ph ∈ Mh .

Further, ah (uh , ∗h vh ) and dh (Th , ρh∗ wh ) are all symmetric, bounded, and positive definite, that is,

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

219

iv. ah (uh , ∗h v h ) = ah (v h , ∗h uh ), ∀uh , v h ∈ Uh ; v. dh (Th , ρh∗ wh ) = dh (wh , ρh∗ Th ), ∀Th , wh ∈ Wh , and there exist positive constants h0  h > 0 such that vi. ah (uh , ∗h uh ) = ν | uh |21 , | ah (uh , ∗h v h ) | νuh 1 v h 1 , ∀uh , v h ∈ Xh ; vii. dh (Th , ρh∗ Th ) = γ0−1 | Th |21 , | dh (Th , ρh∗ wh ) | C˜ 0 Th 1 wh 1 , ∀Th , wh ∈ Wh . Proof. Equation (i) can be proved by using the definitions of (3.3.55) and (3.3.56) and the same method as in the proof of Lemma 3.1.3, while (ii) and (iv)–(vii) are straightforward extensions of Lemma 3.1.3 to the vector form. Thus, it is only necessary to prove (iii). For any v h ∈ Xh and any ph ∈ Mh , it follows from the definition (3.3.55) of bh (·, ·), Lemma 3.1.2, and Green’s formula that

 bh (ph , ∗h v h ) = ∇ph ∗h v h dxdy =

K∈ h

=

Vz ∈ ∗h





K

Vz

∇ph ∗h v h dxdy ∇ph v h dxdy

K∈ h K

=−



ph divv h dxdy

K∈ h K

= −b(ph , v h ), which completes the proof of equation (iii). Next, by using the same approach as in the proof of Lemma 3.1.4, we have the following results in vector form (whose proof is also available in [12,65,73, 80]). Lemma 3.3.4. We have (uh , ∗h v h ) = (v h , ∗h uh ),

∀uh , v h ∈ Xh .

Furthermore, for any u ∈ H m ()2 (m = 0, 1) and v h ∈ Xh , | (u, v h ) − (u, ∗h vh ) | Chm+n um v h n , n = 0, 1. Set ||| uh |||0 = (uh , ∗h uh )1/2 . Then ||| · |||0 is equivalent to ·0 on Xh , namely, there exist two positive constants C1 and C2 such that C1 uh 0 ||| uh |||0  C2 uh 0 , ∀uh ∈ Xh .

(3.3.63)

220 Proper Orthogonal Decomposition Methods for Partial Differential Equations

Remark 3.3.1. For scalar functions, i.e., if uh and v h in Xh are, respectively, substituted with wh and Th in Wh , the results of Lemma 3.3.4 remain valid (see [65, Theorem 3.2.1 and Lemma 5.1.5]). The following results are concerned with the existence, uniqueness, and stability of the solutions for Problem 3.3.4. Theorem 3.3.5. Under the hypotheses of Theorems 3.3.1 and 3.3.2, there exists a unique sequence of solutions (unh , phn , Thn ) (n = 1, 2, · · · , N) for the fully discretized SCNMFVE formulation with second-order time accuracy, i.e., Problem 3.3.4, satisfying √ unh 0 + Thn 0 + k∇unh 0 + k∇Thn 0 + kphn 0  C(ψ0 + k∇ψ |||0 ),

(3.3.64)

which shows that the sequence of solutions of Problem 3.3.4 is stable. Proof. First, we demonstrate that Problem 3.3.4 has a unique sequence of solutions. Because the finite-dimensional subspaces Xh × Mh × Wh are also weakly and sequentially compact Hilbert spaces, by using the same techniques as in the proof that Problem 3.3.3 has a unique sequence of solutions, we can prove that Problem 3.3.4 has a unique sequence of solutions (unh , phn , Thn ) ∈ Xh ×Mh ×Wh . Next, we prove (3.3.64). By taking v h = u¯ nh in the first equation of Problem 3.3.4 and qh = phn in the second equation of Problem 3.3.4 and using Lemmas 3.3.3 and 3.3.4 and the Hölder and Cauchy–Schwarz inequalities, we obtain 1 |||20 ) + kν∇ u¯ nh 20 + kε(phn 20 − ρh phn 20 ) (||| unh |||20 − ||| un−1 h 2 = k(T¯hn j , ∗h u¯ nh ) k kν  T¯hn 2−1 + ∇ u¯ nh 20 . (3.3.65) 2ν 2 It follows from (3.3.65) that |||20 +2kν∇ u¯ nh 20 + 2kε(phn 20 − ρh phn 20 )2 ||| unh |||20 − ||| un−1 h (3.3.66)  kν −1 T¯ n 2−1 . h

Summing (3.3.66) from 1 to n, we have ||| unh |||20 +

n kν

2 n 2 n 2 ∇(uih + ui−1 h )0 + 2kε(ph 0 − ρh ph 0 ) 2 i=1

 kν −1

n

i=1

T¯hi 2−1 .

(3.3.67)

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

221

If phn = 0, then it is easily seen that phn 0 > ρh phn 0 from (3.3.62). Therefore, there exists a constant δ ∈ (0, 1) such that δphn 0 = ρh phn 0 . By extracting  #n #n √ 2 1/2  the square root for (3.3.67), using i=1 bi i=1 | bi | / n, a + b0  a0 − b0 , and Lemma 3.3.4, and then simplifying, we get unh 0

+ k∇unh 0

! n "1/2

√ n i 2 ¯ + kph 0  C k Th −1 .

(3.3.68)

i=1

By taking ϕh = T¯hn in the third equation of Problem 3.3.4 and using Lemmas 3.3.3 and 3.3.4 and the Hölder and Cauchy–Schwarz inequalities, we obtain ||| Thn |||20 − ||| Thn−1 |||20 +

k ∇(Thn + Thn−1 )20 = 0. 2γ0

(3.3.69)

Summing (3.3.69) from 1 to n, we have ||| Thn |||20 +

n k

∇(Thi + Thi−1 )20 =||| ψ |||20 . 2γ0

(3.3.70)

i=1

By extracting the square root of (3.3.70) and using a + b0  a0 − b0 ,  #n #n √ 2 1/2  i=0 bi n=0 | bi | / n, and Lemma 3.3.4, we obtain Thn 0 + k∇Thn 0  C(ψ0 + k∇ψ |||0 ).

(3.3.71)

From (3.3.68) and (3.3.71), using Lemma 3.3.4, we obtain unh 0 + k∇unh 0 +

√ kphn 0  Cψ0  C(ψ0 + k∇ψ |||0 ). (3.3.72)

Combining (3.3.71) with (3.3.72) yields (3.3.64). If phn = 0, (3.3.64) is correct, which completes the proof of Theorem 3.3.5. Set A((Sh u¯ n , Qh p n ); (v h , qh )) = a(Sh u¯ n , v h ) + a1 (Sh u¯ n , Sh u¯ n , v h ) − b(Qh p n , v h ) + b(qh , Sh un ), A((u¯ n , p n ); (v h , qh )) = a(u¯ n , v h ) + a1 (u¯ n , u¯ n , v h ) − b(p n , v h ) + b(qh , un ).

(3.3.73)

By using the stabilized CN mixed FE (SCNMFE) methods (for example, see Section 2.1.1 or [21,80,105]) for the nonstationary Navier–Stokes equation, we immediately obtain the following lemma.

222 Proper Orthogonal Decomposition Methods for Partial Differential Equations

Lemma 3.3.6. Let (Sh un , Qh p n ) be the Navier–Stokes projection of the solutions (un , p n ) for Problem 3.3.3 on Uh × Mh , that is, for the solutions (un , p n ) ∈ U × M of Problem 3.3.3, there exist (Sh un , Qh p n ) (n = 1, 2, · · · , N ) such that kA((Sh u¯ n , Qh p n ); (v h , qh )) + (Sh un − Sh un−1 , v h ) + kDh (Qh p n , qh ) = kA((u¯ n , p n ); (v h , qh )) + (un − un−1 , v h ), ∀(v h , qh ) ∈ Uh × Mh , n = 1, 2, · · · , N,

(3.3.74)

Sh u0 = h u0 (x, y), u0 = u0 (x, y),

(3.3.75)

(x, y) ∈ .

Then, we have Sh un 1 + Qh p n 0  C(un 1 + p n 0 ), 1  n  N. (3.3.76) If k = O(h) and the solution (un , p n ) ∈ H 2 ()2 × H 1 () (n = 1, 2, · · · , N ) for Problem 3.3.3, then we have the following error estimates: un − Sh un 0 + k∇(un − Sh un )0 + kp n − Qh p n 0  Ch2 , n = 1, 2, · · · , N. (3.3.77) Remark 3.3.2. In fact, (3.3.74) and (3.3.75) are the system of error equations between the standard SCNMFE formulation and the semidiscretized CN formulation with respect to time for the nonstationary Navier–Stokes equation, thus (3.3.76) and (3.3.77) are directly obtainable from the SCNMFE method (see, e.g., [21,105,80]). By the FE methods for elliptic equations (see Theorem 2.1.16 or [21,31,80]), we have the following. Lemma 3.3.7. Let Rh : W → Wh be a generalized Ritz projection, i.e., for given unh ∈ Xh , T n−1 ∈ W , and Thn−1 ∈ Wh and for any T n ∈ W (n = 1, 2, · · · , N ), there exist Rh T n ∈ Wh (n = 1, 2, · · · , N ) such that (Rh T n , wh ) + kd(Rh T¯ n , wh ) + ka2 (u¯ nh , Rh T¯ n , wh ) − (Rh T n−1 , wh ) = (T n , wh ) + kd(T¯ n , wh ) + ka2 (u¯ n , T¯ n , wh ) − (T n−1 , wh ), ∀wh ∈ Wh , n = 1, 2, · · · , N.

(3.3.78)

If (un , p n , T n ) (n = 1, 2, · · · , N ) are the solutions of Problem 3.3.3 and T n ∈ H 2 () ∩ W , then we have the following estimates: Rh T n 0 + k∇Rh T n 0  C∇T n 0 , n = 0, 1, 2, · · · , N,

(3.3.79)

Rh T n − T n 0 + k∇(Rh T n − T n )0  Ch2 ψ2 , n = 0, 1, 2, · · · , N.

(3.3.80)

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

223

The convergence of Problem 3.3.4 comes from the following error estimates. Theorem 3.3.8. Let (u, p, T ) be the solution to Problem 3.3.2 and (unh , phn , Thn ) the sequence of solutions of the fully discretized SCNMFVE formulation with second-order time accuracy (i.e., Problem 3.3.4). Then, under the same hypotheses of Theorems 3.3.2 and 3.3.5, if ph0 = p 0 = 0 (or ph0 = Qh p 0 ), h = O(k), N0 ν −1 ∇ u¯ nh 0  1/4, and ψ ∈ H 1 (), we have the following error estimates: u(tn ) − unh 0 + T (tn ) − Thn 0 + k[p(tn ) − phn 0 + ∇(u(tn ) − unh )0 + ∇(T (tn ) − Thn )0 ]  C(h2 + k 2 )ψ1 , n = 1, 2, · · · , N.

(3.3.81)

Proof. By subtracting Eqs. of Problem 3.3.4 from those of Problem 3.3.3 by taking v = v h , q = qh , and ϕ = ϕh and by using Lemmas 3.3.3 and 3.3.4, we obtain the following system of error equations: ⎧ ⎪ (un − unh , v h ) + (unh − ∗h unh , v h − ∗h v h ) + ka(u¯ n − u¯ nh , v h ) ⎪ ⎪ ⎪ ⎪ ⎪ + ka1 (u¯ n , u¯ n , v h ) − ka1h (u¯ nh , u¯ nh , ∗ v h ) − kb(p n − phn , v h ) ⎪ ⎪ ⎪ ⎪ ⎪ = k((T¯ n − T¯hn )j , v h ) − k(T¯hn j , ∗h v h − v h ) + (un−1 − un−1 ⎪ h , vh) ⎪ ⎪ ⎪ n−1 ∗ ⎪ − (uh , h v h − v h ), ∀v h ∈ Xh , n = 1, 2, · · · , N, ⎪ ⎪ ⎪ ⎨ b(q , un − un ) − ε(p n −  p n , q −  q ) h h h h h h h h (3.3.82) = 0, ∀qh ∈ Mh , n = 1, 2, · · · , N, ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ (T n − Thn , ϕh ) − (Thn , ρh∗ ϕh − ϕh ) + kd(T¯ n − T¯hn , ϕh ) ⎪ ⎪ ⎪ ⎪ ⎪ + ka2 (u¯ n , T¯ n , ϕh ) − ka2h (u¯ nh , T¯hn , ρh∗ ϕh ) = (T n−1 − Thn−1 , ϕh ) ⎪ ⎪ ⎪ ⎪ ⎪ − (Thn−1 , ρh∗ ϕh − ϕh ), ∀ϕh ∈ Wh , n = 1, 2, · · · , N, ⎪ ⎪ ⎪ ⎩ 0 u − u0h = 0, T 0 − Th0 = ψ(x, y) − ρh ψ(x, y), (x, y) ∈ . ¯ = Sh u¯ n − u¯ n . Using (3.3.74) Set ζ n = Qh p n − phn , E n = Sh un − unh , and E h and the system of error equations (3.3.82), we have n

1 n 2 ¯ n |21 E 0 + kν | E 2 ¯ n ) + ka(Sh u¯ n − u¯ n , E ¯ n) = (Sh un − un , E ¯ n ) + ka(u¯ n − u¯ n , E ¯ n ) − 1 (E n−1 , E n ) + (un − unh , E h 2 n−1 n−1 ¯ n n n ¯n ¯ n) = (Sh u − u , E ) + kb(Qh p − p , E ) + ka1 (u¯ n , u¯ n , E ¯ n ) − ka1 (u¯ n , u¯ n , E ¯ n ) + ka1h (u¯ n , u¯ n , ∗h E ¯ n) − ka1 (Sh u¯ n , Sh u¯ n , E h h ¯ n ) − (un − ∗h un , E ¯ n − ∗h E ¯ n) + kb(p n − p n , E h

h

h

1 ∗ ¯ n−1 n ¯ + (un−1 − ∗h un−1 − un−1 h h , E − h E ) − (Sh u h ,E ) 2 ∗ ¯n ¯n ¯n ¯n ¯n ¯n ¯n + (un−1 − un−1 h , E ) + k((T − Th )j , E ) − k(Th j , h E − E ) n

n

224 Proper Orthogonal Decomposition Methods for Partial Differential Equations

¯ n ) − k(T¯ n j , ∗h E ¯ n ) − ka1 (E ¯ nh , u¯ n , E ¯ n) ¯n−E = k((T¯ n − T¯hn )j , E h h ¯ n − ∗h E ¯ n) ¯ n ) + 1 (E n−1 , E n−1 ) + kb(ζ n , E + ka1h (u¯ nh , u¯ nh , E 2 n−1 ¯ n ∗ n ¯ n ). − (unh − un−1 −  (u − u ), E − ∗h E (3.3.83) h h h h By Lemma 3.3.4 and the Hölder and Cauchy–Schwarz inequalities, we have ¯ n ) − k(T n j , ∗h E ¯ n) | ¯n−E | k((T¯ n − T¯hn )j , E h ¯ n 0 + Ckh2 ∇ T¯ n 0 ∇ E ¯ n 0  CkT¯ n − T¯ n 0 E h

h

νk ¯ n 20 . (3.3.84)  CkT¯ n − T¯hn 20 + CkE¯n 20 + Ckh4 + ∇ E 8 If k = O(h), by using the inverse error estimate and Taylor’s formula, we get ∗ ¯ ¯ − ∗h (unh − un−1 | (unh − un−1 h h ), E − h E ) | ¯ n |1  Ch2 un − un−1 1 | E n

h

n

h

 Ch3 (∇E n 20 + ∇(Sh un − Sh un−1 )20 + ∇E n−1 20 ) +  ChE n 20 + Ck 2 h3 + ChE n−1 20 +

kν ¯ n 20 . ∇ E 8

kν ¯ n 20 ∇ E 8 (3.3.85)

Noting that b(qh , Sh un − un ) = −kε(Qh p n − ρh (Qh p n ), qh − ρh qh ), by the properties of the operator h and the second equation of (3.3.82), we have ¯ n ) = b(ζ n , Sh u¯ n − u¯ n ) + b(ζ n , u¯ n − u¯ n ) b(ζ n , E h ε n−1 ε n n n n − h ζ n−1 , ζ n − h ζ n ) = − (ζ − h ζ , ζ − h ζ ) − (ζ 2 2 ε ε  − ζ n − h ζ n 20 + ζ n−1 − h ζ n−1 20 . (3.3.86) 4 4 If N0 ν −1 ∇ u¯ nh 0  1/4 (n = 1, 2, · · · , N ), by Lemma 3.3.4, (3.3.3), and (3.3.8), we have ¯ n − ∗h E ¯ n , u¯ n , E¯ n ) | ¯ n ) − a1 (E k | a1h (u¯ nh , u¯ nh , E h kν n 2 4 ¯  Ckh + ∇ E 0 . 4

(3.3.87)

Combining (3.3.84)–(3.3.87) with (3.3.83), we have ¯ n 20 + kε ζ n − h ζ n 20 − kε ζ n−1 − h ζ n−1 20 E n 20 + kν∇ E 2 2 4 2 3 n−1 2 n−1 2  Ckh + Ck h + E 0 + ChE 0 + ChE n 20 + CkT¯ n − T¯ n 20 . (3.3.88) h

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

225

By summing (3.3.88) from 1 to n, if h is sufficiently small such that Ch  1/2 in (3.3.88) and ph0 = p 0 = 0 (or ph0 = Qh p 0 ), we obtain E n 20 + 2kν

n

¯ i 20 + kεζ n − h ζ n 20 E

i=1

 Cnkh + Ck 4

n

T¯ i − T¯hi 20 + Ck

i=1

n−1

E i 20 .

(3.3.89)

i=0

Applying the discrete Gronwall lemma, i.e., Lemma 1.4.1, to (3.3.89), we have E n 20 + k

n

¯ i 20 + kζ n − h ζ n 20 ∇ E

i=1

!

C h +k 4

n

" ¯i

T

− T¯hi 20

exp(Ckn).

(3.3.90)

i=1

By extracting the square root of (3.3.90) and using a + b0  a0 − b0 and  #n #n √ 2 1/2  i=0 bi i=0 | bi | / n, we obtain ¯ n 0 + k(ζ n 0 − h ζ n 0 ) E n 0 + k∇ E "1/2 ! n

4 i i 2 ¯ ¯  T − T 0 . C h +k h

(3.3.91)

i=1

If ζ n = 0, then ζ n 0 > h ζ n 0 . Thus, there is a constant ω ∈ (0, 1) such that ωζ n 0 = h ζ n 0 . Then, by using the triangle inequality, (3.3.91), and Lemma 3.3.6, we have un − unh 0 + k[∇(un − unh )0 + p n − phn 0 ] ! "1/2 n

4 i i 2 ¯ ¯ C h +k  T − T 0 . h

(3.3.92)

i=1

Set en = Rh T n − Thn . On the one hand, by using the system of error equations (3.3.82) and Lemma 3.3.7, we obtain 1 1 en 20 + kγ0−1 ∇ e¯n 20 = (en , e¯n ) + kd(e¯n , e¯n ) − (en , en−1 ) 2 2 = [(Rh T n − T n , e¯n ) + kd(Rh T¯ n − T¯ n , e¯n )] 1 + [(T n − Thn , e¯n ) + kd(T¯ n − T¯hn , e¯n )] − (en , en−1 ) 2 = [(Rh T n−1 − T n−1 , e¯n ) + ka2 (u¯ n , T¯ n , e¯n ) − ka2 (u¯ nh , Rh T¯ n , e¯n )] + [(Thn , ρh∗ e¯n − e¯n ) + ka2h (u¯ nh , T¯hn , ρh∗ e¯n ) − ka2 (u¯ n , T¯ n , e¯n )

226 Proper Orthogonal Decomposition Methods for Partial Differential Equations

1 + (T n−1 − Thn−1 , e¯n ) − (Thn−1 , ρh∗ e¯n − e¯n )] − (en−1 , en ) 2 1 = (en−1 , en ) + (Thn − Thn−1 , ρh∗ e¯n − e¯n ) 2 + ka2h (u¯ nh , T¯hn , ρh∗ e¯n ) − ka2 (u¯ nh , Rh T¯ n , e¯n ). (3.3.93) By Lemma 3.3.4 and the Hölder and Cauchy–Schwarz inequalities, we obtain | (Thn − Thn−1 , ρh∗ e¯n − e¯n ) |  Ch(en 0 + Rh T n − T n 0 + hT n − T n−1 1 + T n−1 − Rh T n−1 0 + en−1 0 )∇ e¯n 0  k  Ch h4 + k 2 h2 + en 20 + en−1 20 + ∇ e¯n 20 , (3.3.94) 4γ0 1 1 1 1 (en−1 , en )  en−1 0 en 0  en−1 20 + en 20 . (3.3.95) 2 2 4 4 If N0 ν −1 ∇ u¯ nh 0  1/4 (n = 1, 2, · · · , N ), by using Lemmas 3.3.3 and 3.3.4, (3.3.4), and the Hölder and Cauchy–Schwarz inequalities, we have k k | a2h (u¯ nh , T¯hn , ρh∗ e¯n ) − a2 (u¯ nh , Rh T¯ n , e¯n ) | ∇ e¯n 20 + Ckh4 . (3.3.96) 4γ0 If k = O(h), by combining (3.3.94)–(3.3.96) with (3.3.93), we have  en 20 + k∇ e¯n 20  Ck h4 + en 20 + en−1 20 + en−1 20 .

(3.3.97)

Summing (3.3.97) from 1 to n and using Lemma 3.3.7 and (3.3.48), we obtain en 20 + kγ0−1

n

∇ e¯i 20  Cnkh4 + e0 20 + Ck

i=1

 Ch4 + Ck

n

n

ei 20

i=1

ei 20 + CRh ψ − ψ20 + Cψ − ρh ψ20

i=1

 Ch4 + Ck

n

ei 20 .

(3.3.98)

i=1

If k is sufficiently small such that Ck  1/2 in (3.3.98), we obtain en 20

+ kγ0−1

n

i=1

∇ e¯i 20

 Ch + Ck 4

n−1

i=0

ei 20 .

(3.3.99)

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

227

Application of the discrete Gronwall lemma, i.e., Lemma 1.4.1, to (3.3.99) yields en 20 + kγ0−1

n

∇ e¯i 20  Ch4 exp(Cnk)  Ch4 .

(3.3.100)

i=1

By extracting the square root of (3.3.100) and using a + b0  a0 − b0  #n #n √ 2 1/2  and i=1 bi i=1 | bi | / n, we get en 0 + k∇en 0  Ch2 .

(3.3.101)

By using the triangle inequality, (3.3.101), and Lemma 3.3.7, we obtain T n − Thn 0 + k∇(T n − Thn )0  Ch2 .

(3.3.102)

Combining (3.3.92) with (3.3.102), we have arrived at un − unh 0 + k[∇(un − unh )0 + p n − phn 0 ]  Ch2 .

(3.3.103)

Combining (3.3.102) and (3.3.103) with Theorem 3.3.2 gives (3.3.81). If ζ n = 0, (3.3.81) is also valid. This completes the proof of Theorem 3.3.8. Remark 3.3.3. It is known from Theorem 3.3.5 and its proof that if ψ1 is sufficiently small, then the conditions N0 ν −1 ∇ u¯ nh 0  1/4 (n = 1, 2, · · · , N ) in Theorem 3.3.8 hold. Remark 3.3.4. Problem 3.3.4 is formulated directly from the semidiscretized scheme with respect to time so that the discussion for the semidiscretized SCNMFVE method with respect to spatial variables can be bypassed and becomes simpler (see, e.g., [49,188]). As long as the Reynolds number Re, the Prandtl coefficient P r, ψ, the time step k, the spatial grid parameter h, and the trail function spaces Xh , Mh , and Wh are provided, a sequence of solutions (unh , phn , Thn ) ∈ Xh × Mh × Wh (1  n  N ) is obtained by solving Problem 3.3.4. The first L solutions (unh , phn , Thn ) (1  n  L, in general, L  N √ and L < 5, for example, L = 20, N = 5000) are taken from N solutions (unh , phn , Thn ) (1  n  N ) as snapshots.

3.3.5 Formulations of POD Bases and the Reduced-Order Model for the 2D Nonstationary Incompressible Boussinesq Equation In this subsection, we employ the POD method from Sections 2.2.1 and 2.2.3 to formulate the POD basis (for more details, see [63,93,118,145]) and build the PODROSCNEMFVE model for the 2D nonstationary incompressible Boussinesq equation.

228 Proper Orthogonal Decomposition Methods for Partial Differential Equations

For (unh , phn , Thn ) (n = 1, 2, · · · , L) in Section 3.3.4, set U i = (unh , phn , Thn ) (n = 1, 2, · · · , L) and V = span{U 1 , U 2 , · · · , U L },

(3.3.104)

referred to as the subspace spanned by the snapshots {U i }L i=1 , at least one of which is a nonzero vector-valued function. Let {ψ j }lj =1 represent an orthonormal basis of V with l = dimV. Then each vector of the set V can be denoted by Ui =

l

(U i , ψ j )Xˆ ψ j , i = 1, 2, · · · , L,

(3.3.105)

j =1

where Xˆ = X × M × W , (U i , ψ j )Xˆ = (∇uih , ∇ψ uj ) + (phi , ψpj ) + (∇Thi , ∇ψTj ), ψ j = (ψ uj , ψpj , ψTj ), and (·, ·) is the L2 -inner product. Definition 3.3.1 (POD technique and POD basis). The POD technique consists in finding an orthonormal basis ψ j (j = 1, 2, · · · , l) satisfying

1

(U i , ψ j )U ×M ψ j min Ui − {ψ j }dj =1 L L

d

j =1

i=1

2

,

(3.3.106)

Xˆ

subject to (ψ r , ψ j )Xˆ = δrj , 1  r  d, 1  j  r,

(3.3.107)

where U i 2ˆ = ∇uih 20 + phi 20 + ∇Thi 20 . A set of solutions {ψ j }dj =1 of X (3.3.106)–(3.3.107) is referred to as a POD basis with rank d. We compose a correlation matrix A = (Aij )L×L ∈ R L×L by Aij = j j j [(∇uih , ∇uh ) + (phi , ph ) + (∇Thi , ∇Th )]/L. Due to the fact that the matrix A with rank l is positive semidefinite, the set of solutions {ψ j }dj =1 for (3.3.106)–(3.3.107) can be found, with the following properties (see Propositions 2.2.4 and 2.3.8 or [63,93,118]). Lemma 3.3.9. Let λ1  λ2  · · ·  λl > 0 be the positive eigenvalues of A and v 1 , v 2 , · · · , v l the corresponding orthonormal eigenvectors. Then a set of POD bases is determined by 1 ψi = √ (U 1 , U 2 , · · · , U L ) · v i , Lλi

1  i  d  l.

(3.3.108)

Moreover, we have the following error equation: L d l



1

U i − (U i , ψ j )Xˆ ψ j 2Xˆ = λj . L i=1

j =1

j =d+1

(3.3.109)

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

229

) ) * * Set X d = span ψ u1 , ψ u2 , · · · , ψ ud , M d = span ψp1 , ψp2 , · · · , ψpd , and W d = span {ψT 1 , ψT 2 , · · · , ψT d }. For uh ∈ Xh , ph ∈ Mh , and Th ∈ Wh , we define three operators P d : Xh → X d , Qd : Mh → M d , and R d : Wh → W d by, respectively, the following relations: (∇P d uh , ∇wd ) = (∇uh , ∇w d ), ∀wd ∈ X d ,

(3.3.110)

(Qd ph , pd ) = (ph , qd ), ∀qd ∈ M d ,

(3.3.111)

(∇R d Th , ∇φd ) = (∇Th , ∇φd ), ∀φd ∈ W d .

(3.3.112)

Thus, by functional analysis (see, e.g., [143]), there exist three extensions P h : X → Xh , Qh : M → Mh , and R h : W → Wh of P d , Qd , and R d satisfying P h |Xh = P d : Xh → X d , Qh |Mh = Qd : Mh → M d , and R h |Wh = R d : Wh → W d defined, respectively, by (∇P h u, ∇wh ) = (∇u, ∇w h ), ∀wh ∈ Xh ,

(3.3.113)

(Qh p, ph ) = (p, qh ), ∀qh ∈ Mh ,

(3.3.114)

(∇R h T , ∇φh ) = (∇T , ∇φh ), ∀φh ∈ Wh ,

(3.3.115)

where (u, p, T ) ∈ X × M × W . Thanks to (3.3.113), (3.3.114), and (3.3.115), the operators P h , Qh , and R h all are bounded, i.e., ∇(P h u)0  ∇u0 , ∀u ∈ X,

(3.3.116)

Qh p0  p0 , ∀p ∈ M,

(3.3.117)

∇(P h T )0  ∇T 0 , ∀T ∈ W.

(3.3.118)

Furthermore, the following properties (see [80,93,118]) are satisfied: u − P h u0  Ch∇(u − P h u)0 , ∀u ∈ X,

(3.3.119)

u − P h u−1  Chu − P h u0 , ∀u ∈ X,

(3.3.120)

T − R h T 0  Ch∇(T − R h T )0 , ∀T ∈ W,

(3.3.121)

T − R h T −1  ChT − R h T 0 , ∀T ∈ W.

(3.3.122)

230 Proper Orthogonal Decomposition Methods for Partial Differential Equations

In addition, we have the following consequences (see Theorems 2.1.19 and 2.1.20 or [63,80,93,118]). Lemma 3.3.10. For every d (1  d  l), the operators P d , Qd , and R d satisfy, respectively, l L 

1  i λj , (3.3.123) uh − P d uih 20 + h2 ∇(uih − P d uih )20  Ch2 L j =d+1

i=1

1 L

L

phi − Qd phi 20 

l

(3.3.124)

λj ,

j =d+1

i=1

l L 

1  i λj , Th − R d Thi 20 + h2 ∇(Thi − R d Thi )20  Ch2 L

(3.3.125)

j =d+1

i=1

where (uih , phi , Thi ) ∈ Xh × Mh × Wh (i = 1, 2, · · · , L) are the sequence of solutions to Problem 3.3.4. Moreover, assume that (u, p, T ) ∈ H 2 (0, tN ; H 2 ()2 )2 × H 1 (0, tN ; H m ()) × H 2 (0, tN ; H 2 ()) is the solution to Problem 3.3.2 and that these and their projections by P h , Qh , and R h satisfy, respectively, the following error estimates: u(tn ) − P h u(tn )−1 + hu(tn ) − P h u(tn )0 + h2 ∇(u(tn ) − P h u(tn ))0  Ch3 ,

n = 1, 2, · · · , N,

(3.3.126)

p(tn ) − Qh p(tn )s  Chm−s , n = 1, 2, · · · , N, s = −1, 0, m = 1, 2,

(3.3.127)

T (tn ) − P h T (tn )−1 + hT (tn ) − P h T (tn )0 + h2 ∇(T (tn ) − P h T (tn ))0  Ch3 ,

n = 1, 2, · · · , N.

(3.3.128)

Thus, with regard to the function spaces setting X d × M d × W d , the PODROSCNEMFVE model for the 2D nonstationary incompressible Boussinesq equation is described as follows. Problem 3.3.5. Find (und , pdn , Tdn ) ∈ X d × M d × W d (n = 1, 2, · · · , N ) satisfying (und , pdn , Tdn ) =

d

((∇ψ uj , ∇unh )ψ uj , (ψpj , phn )ψpj , (∇ψTj , ∇Thn )ψTj ),

j =1

n = 1, 2, · · · , L,

(3.3.129)

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

231

(∂t und , ∗h v d ) + ah (u¯ nd , ∗h v d ) + a1h (u¯ nd , u¯ nd , ∗h v d ) + bh (pdn , ∗h v d ) = (T¯ n j , ∗h v d ), ∀v d ∈ X d , L + 1  n  N, (3.3.130) d

b(und , qd ) + D(pdn , qd ) = 0, ∀qd ∈ M d , L + 1  n  N,

(3.3.131)

(∂t Tdn , ρh∗ φh ) + dh (T¯dn , ρh∗ φd ) + a2h (u¯ nd , T¯dn , ρh∗ φd ) = 0, ∀φd ∈ W d , L + 1  n  N,

(3.3.132)

where (uih , phi , Thi ) ∈ Xh × Mh × Wh (i = 1, 2, · · · , L) are the first L solutions to Problem 3.3.4. Remark 3.3.5. It is easy to see that Problem 3.3.4 at each time node contains 4Nh (where Nh is the number of vertices of triangles in h ; see [21,31,80]) degrees of freedom, but Problem 3.3.5 at the same time node contains only 4d (d  l  L  N  Nh ) degrees of freedom. For real-world engineering problems, the number Nh of vertices of triangles in h is easily more than hundreds of thousands or even a few hundred millions, whereas d is only the number of a few eigenvalues and is very small (in the numerical example to be given in Section 3.3.8, d = 6, but Nh = 136×104 ). Problem 3.3.5 is thus the PODROSCNEMFVE model with very few degrees of freedom for the 2D nonstationary incompressible Boussinesq equation. In particular, Problem 3.3.5 uses only the first few known L solutions of Problem 3.3.4 and from them we then seek other (N − L) solutions; we do not have to do repetitive computations. In other words, the first L PODROSCNEMFVE solutions are obtained by projecting the first L classical SCNMFVE solutions onto a POD basis, while the rest (N − L) PODROSCNEMFVE solutions at (N − L) time instants are obtained by extrapolating and iterating equations (3.3.130), (3.3.131), and (3.3.132). Therefore, it is completely different from the existing POD-based reduced-order formulations (see, e.g., [14,28,63,71,89,93,104,118,119,121,135,137,156,177]).

3.3.6 Existence, Uniqueness, Stability, and Convergence of the Reduced-Order Solutions for the 2D Nonstationary Incompressible Boussinesq Equation In the following, we employ the classical FVE method to discuss the existence, uniqueness, stability, and errors of the solutions for the PODROSCNEMFVE model, i.e., Problem 3.3.5, of the 2D nonstationary incompressible Boussinesq equation. We have the following main results of the existence, uniqueness, stability, and convergence. Theorem 3.3.11. Under the same hypotheses of Theorems 3.3.1, 3.3.2, and 3.3.8, there exists a unique sequence of solutions (und , pdn , Tdn ) ∈ X d × M d × W d to Problem 3.3.5 satisfying √ und 0 + Tdn 0 + k∇und 0 + k∇Tdn 0 + kpdn 0  Cψ0 , n = 1, 2, · · · , L, L + 1, · · · , N. (3.3.133)

232 Proper Orthogonal Decomposition Methods for Partial Differential Equations

Therefore, the sequence of solutions (und , pdn , Tdn ) (n = 1, 2, · · · , N ) to Problem 3.3.5 is stable. If k = O(h) and N0 ν −1 ∇ u¯ nd 0  1/4 (n = L + 1, L + 2, · · · , N ), we have the following error estimates: unh − und 0 + Thn − Tdn 0 + k∇(unh − und )0 + k∇(Thn − Tdn )0 ⎛ ⎞1/2 l

√ + kphn − pdn 0  CLk ⎝ λj ⎠ , n = 1, 2, · · · , L, (3.3.134) j =d+1

unh − und 0 + k[∇(unh − und )0 + phn − pdn 0 ] ⎛ ⎞1/2 l

 λj ⎠ ,  C(k 2 + h2 ) k(n − L) + C(k 2 + h2 ) + CLk ⎝ j =d+1

n = L + 1, L + 2, · · · , N, ⎛ Thn − Tdn 0 + k∇(Thn − Tdn )0  C(k 2 + h2 ) + CLk ⎝

l

(3.3.135) ⎞1/2 λj ⎠

,

j =d+1

n = L + 1, L + 2, · · · , N.

(3.3.136)

Proof. (1) We first prove the result (3.3.133) implying the existence, uniqueness, and stability for Problem 3.3.5. If (3.3.133) holds, then when ψ = 0, we must have und = 0, pdn = 0, and n Td = 0 (n = 1, 2, · · · , N ). Thus, there exists a unique sequence of solutions (und , pdn , Tdn ) ∈ X d × M d × W d to Problem 3.3.5. If n = 1, 2, · · · , L, by (3.3.116)–(3.3.118) and Theorem 3.3.5, (3.3.133) holds when n = 1, 2, · · · , L. If n = L + 1, L + 2, · · · , N , by taking v d = u¯ nd in (3.3.130) and qd = pdn in (3.3.131) and by using Lemmas 3.3.3 and 3.3.4 and the Hölder and the Cauchy– Schwarz inequalities, we obtain 1 |||20 ) + kν∇ u¯ nd 20 + kεpdn − h phd 20 (||| und |||20 − ||| un−1 d 2 kν  CkT¯dn 2−1 + ∇ u¯ nd 20 . (3.3.137) 2 It follows from (3.3.137) that |||20 +2kν∇ u¯ nd 20 + 2kεpdn − h pdn 20 ||| und |||20 − ||| un−1 d  CkT¯ n 2−1 . (3.3.138) d

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

233

By summing (3.3.138) from L to n, we get ||| und |||20 +

n kν

∇ u¯ id 20 + 2kεpdn − h pdn 20 2 i=L

||| uL−1 |||20 +Ck d

n

T¯di 2−1 .

(3.3.139)

i=L

If pdn = 0, then it is easy to see that pdn 0 > h pdn 0 from (3.3.62). Therefore, there exists a constant δ ∈ (0, 1) such that δpdn 0  h pdn 0 . By extracting the #n  #n √ 2 1/2  square root of (3.3.139), using i=1 bi i=1 | bi | / n, a − b0  a0 − b0 , and Lemma 3.3.4 and then simplifying, we obtain √ und 0 + k∇und 0 + kpdn 0 ! "1/2 n

L−1 2 2 L 2 i 2 ¯  C u 0 + k ud 1 + k T −1 . (3.3.140) d

d

i=L

Taking φd = T¯dn in (3.3.132) and using Lemmas 3.3.3 and 3.3.4 and the Hölder and Cauchy–Schwarz inequalities, we obtain ||| Tdn |||20 +

2k ∇ T¯dn 20 ||| Tdn−1 |||20 . γ0

(3.3.141)

Sum (3.3.141) from L to n to get 2k

∇ T¯di 20 ||| T L−1 |||20 . γ0 n

||| Tdn |||20 +

(3.3.142)

i=L

#n  2 1/2 By# extracting the√square root of (3.3.142) and using Lemma 3.3.4, i=1 bi  ni=1 | bi | / n, a − b0  a0 − b0 , and (3.3.133) when n = L, we have Tdn 0 + k∇Tdn 0  Cψ0 , n = L, L + 1, . . . , N.

(3.3.143)

Because Tdn −1  CTdn 0 , from (3.3.140) and (3.3.143) as well as Lemma 3.3.4, we obtain und 0 + k∇und 0 +

√ kpdn 0  Cψ0 , n = L, L + 1, ..., N. (3.3.144)

Combining (3.3.143) with (3.3.144), we arrive at (3.3.133) for n = L + 1, L + 2, · · · , N . If pdn = 0, (3.3.133) remains valid. (2) Now, we prove the error estimates (3.3.134) and (3.3.135). If n = 1, 2, · · · , L, by Lemma 3.3.10, we immediately obtain (3.3.134).

234 Proper Orthogonal Decomposition Methods for Partial Differential Equations

If n = L + 1, L + 2, · · · , N , by subtracting Eqs. of Problem 3.3.5 from those of Problem 3.3.4, taking v h = v d , qh = qd , and φh = φd , and applying Lemmas 3.3.3 and 3.3.4, we obtain the following system of error equations: (unh − und , ∗h v d ) + ka(u¯ nh − u¯ nd , v d ) + ka1h (u¯ nh , u¯ nh , ∗h v d ) − ka1h (u¯ nd , u¯ nd , ∗h v d ) − kb(phn − pdn , v d ) = k((T¯ n − T¯ n )j , ∗h v d ) + (un−1 − un−1 , ∗h v d ), ∀v d ∈ X d , h

d

h

d

n = L + 1, L + 2, · · · , N,

(3.3.145)

b(qd , unh − und ) + ε(phn − pdn − h (phn − pdn ), qd − h qd ) = 0, ∀qd ∈ M d , n = L + 1, L + 2, · · · , N, (3.3.146) (Thn − Tdn , ρh∗ φd ) + kd(T¯hn − T¯dn , φd ) + ka2h (u¯ nh , T¯hn , ρh∗ φd ) − ka2h (u¯ n , T¯ n , ρh∗ φd ) = (T n−1 − T n−1 , ρh∗ φd ), ∀φd ∈ W d , d

d

h

d

n = L + 1, L + 2, · · · , N.

(3.3.147)

Set en = P d unh − und , ρ n = unh − P d unh , e¯ n = P d u¯ nh − u¯ nd , ρ¯ n = u¯ nh − P d u¯ nh , = Qd phn − pdn , and ξ n = phn − Qd phn . By (3.3.110) and Green’s formula, we have b(ξ n , e¯ n ) = (phn − Qd phn , div¯en ) = 0. Thus, on the one hand, using (3.3.110), the error equations (3.3.145) and (3.3.146), and Lemma 3.3.4, we have  1 ||| en |||20 − ||| en−1 |||20 + kν∇ e¯ n 20 2 ∗ n ¯ ) + ka(P d u¯ nh − u¯ nd , e¯ n ) = (P d unh − und − (P d un−1 − un−1 h d ),  e

ηn

∗ n ¯ ) = (ρ n−1 − ρ n , ∗ e¯ n ) + (unh − und , ∗ e¯ n ) − (un−1 − un−1 h d , e n n n + ka(u¯ h − u¯ d , e¯ )

= (ρ n−1 − ρ n , e¯ n ) + (ρ n − ρ n−1 , e¯ n − ∗h e¯ n ) + kb(phn − pdn , e¯ n ) + ka1h (u¯ nd , u¯ nd , ∗h e¯ n ) − ka1h (u¯ nh , u¯ nh , ∗h e¯ n ) + k((T¯hn − T¯dn )j , ∗h e¯ n ) = (ρ n−1 − ρ n , e¯ n ) + (ρ n − ρ n−1 , e¯ n − ∗h e¯ n ) + kb(phn − pdn , e¯ n ) + ka1h (u¯ nd , u¯ nd , ∗h e¯ n ) − ka1h (u¯ nh , u¯ nh , ∗h e¯ n ) + k((T¯hn − T¯dn )j , ∗h e¯ n ) + kb(ξ n , e¯ n ) − kε(phn − pdn − h (phn − pdn ), ηn − h η)  C(k −1 ρ n−1 − ρ n 2−1 + kρ n − ρ n−1 20 ) kν + ∇ e¯ n 20 − kεηn − h η20 + ka1h (u¯ nd , u¯ nd , ∗h e¯ n ) 8 − ka1h (u¯ nh , u¯ nh , ∗h e¯ n ) + k((T¯hn − T¯dn )j , ∗h e¯ n ).

(3.3.148)

On the other hand, if N0 ν −1 ∇ u¯ nh 0  1/4 and N0 ν −1 ∇ u¯ nd 0  1/4 (n = 1, 2, · · · , N ), using the properties of a1h (·, ·, ·), the Hölder and Cauchy–Schwarz

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

235

inequalities, and Lemmas 3.3.3 and 3.3.4, we have | ka1h (u¯ nd , u¯ nd , ∗h e¯ n ) − ka1h (u¯ nh , u¯ nh , ∗h e¯ n )) | kν  Ckρ n − ρ n−1 20 + ∇ e¯ n 20 . 4

(3.3.149)

Using the Hölder and Cauchy–Schwarz inequalities, we have kν | k((T¯hn − T¯dn )j , ∗h e¯ n ) | CkT¯hn − T¯dn 2−1 + ∇ e¯ n 20 . 8

(3.3.150)

By combining (3.3.149) and (3.3.150) with (3.3.148) and by using (3.3.120), we have ||| en |||20 − ||| en−1 |||20 +kν∇ e¯ n 20 + 2kεηn − h ηn 20  Ck(ρ n 20 + ρ n−1 20 ) + CkT¯ n − T¯ n 2−1 . (3.3.151) h

d

Summing (3.3.151) from L to n and using Lemma 3.3.4, we obtain en 20 + kν

n

∇ e¯ i 20 + 2kεηn − h ηn 20

i=L

 CeL−1 20

+ Ck

n

(ρ i 20 + T¯hi − T¯di 2−1 ).

(3.3.152)

i=L

By extracting the square root of (3.3.152), using a − b0  a0 − b0 and  #n #n √ 2 1/2  i=1 bi n=1 | bi | / n, and then simplifying, we obtain  en 0 + k∇en 0 + k 1/2 (ηn 0 − h ηn 0 )  C eL−1 0 + k∇eL 0 .1/2 - n

i 2 i i 2 (ρ 0 + T¯ − T¯ −1 ) . (3.3.153) +C k h

d

i=L

If ηn = 0, then ηn 0 > h ηn 0 (n = L, L + 1, · · · , N ). Thus, there are α1 ∈ (0, 1) such that α1 ηn 0  h ηn 0 . We can drop the term h ηn 0 in (3.3.153) and obtain  en 0 + k∇en 0 + k 1/2 ηn 0  C eL−1 0 + k∇eL 0 .1/2 - n

i 2 i i 2 (ρ 0 + T¯ − T¯ −1 ) . (3.3.154) +C k h

d

i=L

Moreover, with Lemma 3.3.10 and Theorem 3.3.5, we get ρ i 0  uih − u(ti )0 + u(ti ) − P h u(ti )0 + P h (u(ti ) − uih )0  C[u(ti ) − P h u(ti )0 + u(ti ) − uih 0 ]  C(h2 + k 2 ).

(3.3.155)

236 Proper Orthogonal Decomposition Methods for Partial Differential Equations

By combining (3.3.155) with (3.3.154) and using Lemma 3.3.10 and (3.3.134), we have  en 0 + k∇en 0 + k 1/2 ηn 0  C(k 2 + h2 ) k(n − L) ⎞1/2 ⎛ - n .1/2 l



i i 2 ¯ ¯ ⎠ ⎝ λj +C k T − T −1 (. 3.3.156) + CLk h

j =d+1

d

i=L

Let En = Qd Thn − Tdn and E¯ n = Qd T¯hn − T¯dn . On the one hand, with the error Eq. (3.3.147) and Lemmas 3.3.3 and 3.3.4, we obtain  1 ||| En |||20 − ||| En−1 |||20 + kγ0−1 ∇ E¯ n 20 2 = (En − En−1 , ρ ∗ E¯ n ) + kd(E¯ n , E¯ n ) = [(Qd T n − T n − (Qd T n−1 − T n−1 ), ρ ∗ E¯ n ) + kd(Qd T¯ n − T¯ n , E¯ n )] h

h

h

h

h

h

+ [(Thn − Tdn − (Thn−1 − Tdn−1 ), ρ ∗ E¯ n ) + kd(T¯hn − T¯dn , E¯ n )] = (Qd T n − T n − (Qd T n−1 − T n−1 ), ρ ∗ E¯ n ) h

h

h

h

+ ka2h (u¯ nd , T¯dn , ρh∗ E¯ n ) − ka2h (u¯ nh , T¯hn , ρh∗ E¯ n ) = (Qd T n − T n − (Qd T n−1 − T n−1 ), E¯ n ) + ka2h (u¯ n , T¯ n , ρh∗ E¯ n ) h

h d d h h n ¯n ∗ ¯ d n n d n−1 − ka2h (u¯ h , Th , ρh En )) + (Q Th − Th − (Q Th − Thn−1 ), ρ ∗ E¯ n  Ck −1 (Qd Thn − Thn 2−1 + Qd Thn−1 − Thn−1 2−1 )

+ Ck(Qd Thn − Thn 20 + Qd Thn−1 − Thn−1 20 ) +

− E¯ n )

k ∇ E¯ n 20 4γ0

+ ka2h (u¯ nd , T¯dn , ρh∗ E¯ n ) − ka2h (u¯ nh , T¯hn , ρh∗ E¯ n ).

(3.3.157)

On the other hand, if N0 ν −1 ∇und 0  1/4 (n = 1, 2, · · · , N ), using Lemmas 3.3.3 and 3.3.4, (3.3.4), and the Hölder and Cauchy–Schwarz inequalities, we have | ka2h (u¯ nd , T¯dn , ρh∗ E¯ n ) − ka2h (u¯ nh , T¯hn , ρh∗ E¯ n ) | k  ∇ E¯ n 20 + Ckh4 . 4γ0

(3.3.158)

By combining (3.3.157) with (3.3.158) and using (3.3.122), Lemma 3.3.10, and Theorem 3.3.5, we obtain ||| En |||20 +kγ0−1 ∇ E¯ n 20  Ck(h4 + k 4 )+ ||| En−1 |||20 . Summing (3.3.159) from L to n yields En 20

+ kγ0−1

n

i=L

∇ E¯ i 20

(3.3.159)

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

237

 Cnk(h4 + k 2 ) + CEL−1 20  C(h4 + k 4 ) + CEL−1 20 .

(3.3.160)

Extracting the #n √square root of (3.3.160) and using (3.3.134) and | b | / n and a − b0  a0 − b0 , we get i i=1 ⎛ En 0 + k∇Ei 0  C(h2 + k 2 ) + CLk ⎝

l

 2 1/2 i=1 bi

#n



⎞1/2 λj ⎠

.

(3.3.161)

j =d+1

Using the triangle inequality, (3.3.161), and Lemma 3.3.10, we obtain Thn − Tdn 0 + k∇(Thn − Tdn )0 ⎛ ⎞1/2 l

 C(h2 + k 2 ) + CLk ⎝ λj ⎠ .

(3.3.162)

j =d+1

This yields (3.3.136). By combining (3.3.156) with (3.3.162) and using Lemma 3.3.10, we arrive at √ unh − und 0 + k∇(unh − und )0 + kphn − pdn 0  C(k 2 + h2 ) ⎛ ⎞1/2 l

 λj ⎠ . (3.3.163) + C(k 2 + h2 ) k(n − L) + CLk ⎝ j =d+1

This yields (3.3.135). If ηn = 0, (3.3.135) also remains valid. This completes the proof of Theorem 3.3.11. Remark 3.3.6. It is easy to see from (3.3.64) in Theorem 3.3.5, (3.3.133) in Theorem 3.3.11, and their proofs that the conditions N0 ν −1 ∇ u¯ nh 0  1/4 and N0 ν −1 ∇ u¯ nd 0  1/4 (n = 1, 2, · · · , N ) hold when ψ0 is sufficiently small. Combining Theorem 3.3.5 with Theorem 3.3.11, we have the following. Theorem 3.3.12. Under the same assumptions of Theorems 3.3.5 and 3.3.11, we have the error estimates between the solution (u, p, T ) to Problem 3.3.2 and the solutions (und , pdn , Tdn ) to Problem 3.3.5 as follows: u(tn ) − und 0 + T (tn ) − Tdn 0 + k[∇(u(tn ) − und )0 + ∇(T (tn ) − Tdn )0 + p(tn ) − pdn 0 ] ⎛ ⎞1/2 l

λj ⎠ , n = 1, 2, · · · , L,  C(k 2 + h2 ) + CLk ⎝ j =d+1

238 Proper Orthogonal Decomposition Methods for Partial Differential Equations

u(tn ) − und 0 + k[∇(u(tn ) − und )0 + p(tn ) − pdn 0 ] ⎛ ⎞1/2 l

 λj ⎠ ,  C(k 2 + h2 ) k(n − L) + C(k 2 + h2 ) + CLk ⎝ j =d+1

n = L + 1, L + 2, · · · , N, T (tn ) − Tdn 0 + k∇(T (tn ) − Tdn )0 ⎛ ⎞1/2 l

 C(k 2 + h2 ) + CLk ⎝ λj ⎠ , n = L + 1, L + 2, · · · , N. j =d+1

# 1/2 l Remark 3.3.7. The factor Lk in Theorems 3.3.11 and 3.3.12 j =d+1 λj is contributed by POD-based order reduction for Problem 3.3.4; it can serve as a guide for choosing the number of POD basis – error control requires d to be # taken such that k 2 L2 lj =d+1 λj = O(k 2 , h4 ).

3.3.7 Algorithm Implementation of the Reduced-Order Model for the 2D Nonstationary Incompressible Boussinesq Equation The algorithm implementation of the PODROSCNEMFVE model for the 2D nonstationary incompressible Boussinesq equation consists of the following seven steps. Step 1. Extraction of snapshots Extract snapshots U i (x, y) = (uih , phi , Thi ) (i = 1, 2, · · · , L  N ) from the classical SCNMFVE solutions for Problem 3.3.4 or from the samples of actual physical system trajectories. Step 2. Formulation of correlation matrix A j

The correlation matrix A = (Aij )L×L , where Aij = [(∇uih , ∇uh ) + j j (phi , ph ) + (∇Thi , ∇Th )]/L. Step 3. Computing eigenvalues and eigenvectors of A Find the eigenvalues λ1  λ2  · · ·  λl > 0 (l = dim{U 1 , U 2 , · · · , U L }) j j j and their associated eigenvectors v j = (a1 , a2 , · · · , aL )τ (j = 1, 2, · · · , l) of A. Step 4. Determination of the number of POD bases For given spatial grid diameter h and time step k, with error δ needed, deter# mine the number d of POD bases to satisfy k 2 + h4 + L2 k 2 lj =d+1 λj  δ 2 .

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

239

Step 5. Formulation of POD basis The POD basis is obtained by ψ j (x, y) = (ψ uj (x, y), ψpj (x, y), ψTj (x, y)) L  

j = ai uih , phi , Thi / Lλj , j = 1, 2, · · · , d. i=1

Step 6. Compute the PODROSCNEMFVE solutions Md = Take X d = span{ψ u1 (x, y), ψ u2 (x, y), · · · , ψ ud (x, y)}, d span{ψp1 (x, y), ψp2 (x, y), · · · , ψpd (x, y)}, and W = span{ψT 1 (x, y), ψT 2 (x, y), · · · , ψT d (x, y)} and solve Problem 3.3.5, obtaining the PODROSCNEMFVE solutions (und , pdn , Tdn ) (n = 1, 2, · · · , L, L + 1, · · · , N ). Step 7. Renewal of POD basis and circulation or end n−1 − und 0  und − un+1 −pdn 0  pdn − pdn+1 0 , and If un−1 d d 0 , pd Tdn−1 − Tdn 0  Tdn − Tdn+1 0 (n = L, L + 1, · · · , N − 1), then (und , pdn , Tdn ) (n = 1, 2, · · · , N ) are the solutions for the PODROSCNEMFVE model satisfying a given accuracy requirement. Else, namely, if un−1 − und 0 < d n+1 n−1 n+1 n−1 n n n n ud − ud 0 or pd − pd 0 < pd − pd 0 or Td − Td 0 < Tdn − Tdn+1 0 (n = L, L + 1, · · · , N − 1), let U i = (uid , pdi , Tdn ) (i = n − L, n − L − 1, · · · , n − 1), return to Step 2.

3.3.8 Numerical Experiments for the 2D Nonstationary Incompressible Boussinesq Equation In the following, we provide some numerical experiments to show the feasibility and efficiency for the PODROSCNEMFVE model for the 2D nonstationary incompressible Boussinesq equation and to validate that the numerical results are consistent with theoretical ones. ¯ consists of a channel of width 6 and length 20 The computational domain  and two rectangular protrusions at the top and bottom of the channel as illustrated in Fig. 3.3.1. The two rectangular protrusions have width 2 and length 4. It is first divided into small squares with side length x = y = 0.01, and then each square is divided by a diagonal along the same direction into two trian√ gles. This constitutes the triangularization h with h = 2 × 10−2 . The dual decomposition ∗h is taken in the barycentric form, namely, the barycenter of the right triangle K ∈ h is taken as the node of the dual decomposition. Take Re = 1000, P r = 7, and ε = 1. Except for the inflow from the left boundary with a velocity of u = (0.1(y − 4.5)(5.5 − y), 0)T (x = 0, 4.5  y  5.5) and outflow on the right boundary with a velocity of u = (u1 , u2 )T satisfying u2 = 0 and u1 (x, y, t) = u1 (19, y, t) (19  x  20, 2  y  8, 0  t  tN ), all initial and boundary value conditions are taken as 0. In order to make k = O(h) satisfied, the time step increment is taken as k = 0.01.

240 Proper Orthogonal Decomposition Methods for Partial Differential Equations

FIGURE 3.3.1 The computational domain and the boundary conditions of flow.

FIGURE 3.3.2 The top chart is the classical SCNMFVE solution and the bottom one the PODROSCNEMFVE solution of the velocity u with six POD bases for Re = 1000 and P r = 7 at time t = 50.

First, 20 numerical solutions (unh , phn , Thn ) (n = 1, 2, · · · , 20) obtained from the classical SCNMFVE formulation Problem 3.3.4 are used to formulate the snapshots U i = (uih , phi , Thi ) (i = 1, 2, · · · , 20). Then 20 eigenvalues in decreasing order and the associated 20 eigenvectors are found by using Step 3 in # 1/2 20 Section 3.3.7. By computing, we reckon out the error factor Lk  j =7 λj

3 × 10−2 for k = 0.01 and L = 20. Thus, we only need to take the first six POD bases (ψ uj , ψpj , ψTj ) (j = 1, 2, · · · , 6) for the construction of subspaces Xd × M d × W d and then find the numerical solutions (und , pdn , Tdn ) (n = 5 × 103 , namely, at t = 50) by the PODROSCNEMFVE model according to the seven steps in Section 3.3.7, whose velocity, pressure, and temperature numerical solutions are depicted at the bottom charts in Figs. 3.3.2, 3.3.3, and 3.3.4, while the numerical solutions of velocity, pressure, and temperature obtained by the classical SCNMFVE formulation Problem 3.3.4 are depicted at the top charts in Figs. 3.3.2, 3.3.3, and 3.3.4 at t = 50, respectively. Every pair of charts in Figs. 3.3.2, 3.3.3, and 3.3.4 has shown close similarity to one another.

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

241

FIGURE 3.3.3 The top chart is the classical SCNMFVE solution and the bottom one the PODROSCNEMFVE solution of the pressure p with six POD bases for Re = 1000 and P r = 7 at time t = 50.

FIGURE 3.3.4 The top chart is the classical SCNMFVE solution and the bottom one the PODROSCNEMFVE solution of the temperature T for Re = 1000 and P r = 7 at time t = 50.

Even though the relative errors of the PODROSCNEMFVE solutions during the initial starting time span are slightly larger than those of the classical SCNMFVE solutions, the PODROSCNEMFVE model at each time node contains only 24 degrees of freedom while the classical SCNMFVE formulation has 4 × 136 × 104 degrees of freedom. The degrees of freedom for the PODROSCNEMFVE model are far fewer than those for the classical SCNMFVE formulation. Therefore, the PODROSCNEMFVE model can greatly reduce the TE accumulation in the computational process, alleviate the computational load, save the CPU time, and improve the accuracy of calculation. Therefore, after some time duration, the relative numerical errors of the PODROSCNEMFVE solutions have gone down in comparison with those of the classical SCNMFVE

242 Proper Orthogonal Decomposition Methods for Partial Differential Equations

FIGURE 3.3.5 The two curves show the growing relative errors of the PODROSCNEMFVE solution and the classical SCNMFVE solution of velocity u for Re = 103 and P r = 0.1 on 0  t  50.

solutions (see Figs. 3.3.5, 3.3.6, and 3.3.7). In fact, Figs. 3.3.5, 3.3.6, and 3.3.7 have shown the cases of the TE accumulation on 0  t  5. It is shown that the relative errors of the classical SCNMFVE solutions are much larger than those of the PODROSCNEMFVE solutions obtained from the PODROSCNEMFVE model. According to the growing trends of the relative errors of the classical SCNMFVE solutions, the classical SCNMFVE formulation would appear to not converge after some computational steps, while the error accumulation of the PODROSCNEMFVE model is rather slow such that it could continuously simulate the development of fluid flow. Thus, in this connection, the solutions obtained by the PODROSCNEMFVE model are better than the classical SCNMFVE solutions. In particular, the PODROSCNEMFVE solutions of pressure and temperature are much more accurate than the classical SCNMFVE solutions. It is also shown that the PODROSCNEMFVE model is computationally effective for finding the numerical solutions of the 2D nonstationary incompressible Boussinesq equation. Fig. 3.3.8 shows the absolute errors between solutions obtained from the PODROSCNEMFVE model Problem 3.3.5 with different numbers of POD bases and solutions obtained from the classical SCNMFVE formulation Problem 3.3.4 when t = 50, for P r = 7 and Re = 1000. It shows that the numerical computing results are consistent with those obtained from theory because the theoretical and numerical errors are all O(10−4 ). Further, comparing the classical SCNMFVE formulation Problem 3.3.4 with the PODROSCNEMFVE model Problem 3.3.5 containing six POD bases when t = 50, for P r = 7 and Re = 1000, we have found that the classical SCNMFVE formulation Problem 3.3.4 had 4 × 136 × 104 unknown quantities at each time node and the required computing time is 240 minutes, while the PODROSCNEMFVE model Problem 3.3.5 with six POD bases only had 4 × 6 unknown quantities at the same time node, and the associated computing time is less than 60 seconds, so the computing time of the classical SCNMFVE formulation Problem 3.3.4 is about 240 times more. Thus, the PODROSCNEMFVE model

Reduced-Order Extrapolation Finite Volume Element Methods Chapter | 3

243

FIGURE 3.3.6 The two curves show the growing relative errors of the PODROSCNEMFVE solution and the classical SCNMFVE solution of pressure p for Re = 103 and P r = 0.1 on 0  t  50.

FIGURE 3.3.7 The two curves show the growing relative errors of the PODROSCNEMFVE solution and the classical SCNMFVE solution of temperature T for Re = 103 and P r = 0.1 on 0  t  50.

FIGURE 3.3.8 Absolute errors between the classical SCNMFVE solution and the PODROSCNEMFVE solutions with different numbers of POD bases for P r = 7 and Re = 1000 at the time level t = 50.

244 Proper Orthogonal Decomposition Methods for Partial Differential Equations

can significantly save the time of calculations, alleviate the computational load, and also reduce the TE accumulation in the computational process. The feasibility and effectiveness are obvious.

3.4 CONCLUDING REMARKS In this Chapter 3, we have introduced the basic principles and methods of the PODROEFVE and PODROSCNEMFVE algorithms. The given concrete numerical examples have demonstrated the viabilities, advantages, and effectiveness of the PODROEFVE and PODROSCNEMFVE algorithms. For more examples, see references [71,84,85,87,98,99,104,106,108,114,120,162,163].