A reduced-order extrapolation algorithm based on SFVE method and POD technique for non-stationary Stokes equations

Applied Mathematics and Computation 247 (2014) 976–995

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

A reduced-order extrapolation algorithm based on SFVE method and POD technique for non-stationary Stokes equations q Zhendong Luo School of Mathematics and Physics, North China Electric Power University, Beijing 102206, China

a r t i c l e

i n f o

a b s t r a c t

Keywords: Proper orthogonal decomposition method Stabilized finite volume element method Non-stationary Stokes equations Error estimate Numerical simulation

In this paper, a reduced-order extrapolation algorithm (ROEA) based on proper orthogonal decomposition (POD) technique and classical stabilized finite volume element (SFVE) method for non-stationary Stokes equations is established. The error estimates between the ROEA solutions and the classical SFVE solutions and the implementation for solving ROEA are provided. Some numerical examples are used to verify that the results of numerical computation are consistent with theoretical conclusions. Moreover, it is shown that ROEA based on SFVE formulation and POD method is feasible and efficient for solving non-stationary Stokes equations. Ó 2014 Elsevier Inc. All rights reserved.

1. Introduction Let X  R2 be a bounded and connected polygonal domain. Consider the following incompressible non-stationary Stokes equations. Problem I. Find u ¼ ðu1 ; u2 ÞT and p such that, for T > 0,

8 ut  mDu þ rp ¼ f ; > > > < r  u ¼ 0; > uðx; y; tÞ ¼ uðx; y; tÞ; > > : uðx; y; 0Þ ¼ u0 ðx; yÞ;

ðx; y; tÞ 2 X  ð0; T; ðx; y; tÞ 2 X  ð0; T; ðx; y; tÞ 2 @ X  ð0; T;

ð1:1Þ

ðx; yÞ 2 X;

T

where u ¼ ðu1 ; u2 Þ is the velocity vector, p the pressure, T the total time, m ¼ 1=Re; Re the Reynolds number, f ðx; y; tÞ the given body force vector, and uðx; y; tÞ and u0 ðx; yÞ are two given velocity vector functions. For the sake of convenience, without loss of generality, we may as well assume that uðx; y; tÞ ¼ 0 in the following theoretical analysis. The non-stationary Stokes equations constitute an important system of equations in fluid dynamics and are widely and successfully applied to a lot of practical engineering fields (see for instance, mentioned in [1–6]). Since their computational fields for practical engineering problems are usual irregular, it is not easy to find their analytic solutions, namely, exact solutions. The efficient approaches are to find their numerical solutions.

q

This work is jointly supported by National Science Foundation of China (11271127). E-mail address: [email protected]

http://dx.doi.org/10.1016/j.amc.2014.09.057 0096-3003/Ó 2014 Elsevier Inc. All rights reserved.

Z. Luo / Applied Mathematics and Computation 247 (2014) 976–995

977

The finite volume element (FVE) method (see [7–9]) is regarded as one of the most effective numerical methods due to its following advantages. Firstly, it can keep the conservation law of mass or energy. Secondly, its accuracy is higher than that of finite difference (FD) method and it is suitable for computational field with complicated boundary exceed FD method. Thirdly, it has the same accuracy as the finite element (FE) method but is simpler and more convenient than FE method. The FVE method is also known as a box method in an early reference (see [10]) where one discretizes the integral form of conservation law of differential equation by choosing linear (or bilinear) FE space as trial space, or known as a generalized difference method (see [11,12]). The FVE method has been widely applied to finding numerical solutions of different types of partial differential equations, for example, second order elliptic equations, parabolic equations, and Stokes equations (see [7–18]). However, some classical fully discrete FVE formulations for non-stationary Stokes equations generally include many degrees of freedom and the combination of velocity trial function space with pressure one needs to satisfy an important convergence stability condition, namely, the Brezzi–Babuška (BB) inequality (see [18,13]) such that they hold many difficulties in practical engineering computing. Thus, an important problem is how to circumvent the constraint of the BB inequality and alleviate the computational load by saving time-consuming calculations in the computational process in a way that guarantees sufficiently accurate numerical solutions. To circumvent the constraint of BB inequality in FVE formulations for non-stationary Stokes equations, some stabilized finite volume element (SFVE) formulations based on two local Gauss integrals or parameter-free have developed (see [17,19]). These methods are stabilized with the lowest equal-order (namely, P 1  P 1 ) elements by the residual of these local integrals on each triangular element and are free of stabilization parameters, which do not require any calculation of highorder derivatives or edge-based data structures and can be implemented at the element level. It has been shown that the proper orthogonal decomposition (POD) method by combining with some numerical methods for partial differential equations can provide an efficient means of generating reduced-order formulations and alleviating the computational load and memory requirements (see [20]). The POD method has been widely and successfully applied to numerous fields, including signal analysis and pattern recognition (see [21]), statistics (see [22]), geophysical fluid dynamics and meteorology (also see [22]). The POD method essentially provides an orthogonal basis for representing the given data in a certain least squares optimal sense, that is, it provides a way to find optimal lower dimensional approximations to the given data. In early times, the POD method was mainly used to perform principal component analysis in computations of statistics and search the main behavior of dynamic systems (see [20–29]), until the method of snapshots was introduced by Sirovich (see [30]) and was then widely applied to reducing the order of the POD eigenvalue problem. Until ten years ago, some Galerkin POD methods for parabolic problems and a general equation in fluid dynamics have not been presented (see [31,32]), and the singular value decomposition approach combined with POD technique has not been used to treat the Burgers equation (see [33]) and the cavity flow problem (see [34]). More recently, some reduced-order FD formulations and FE (or mixed FE or least-square mixed FE) formulations and error estimates for the upper tropical Pacific Ocean formulation, parabolic problems, Burgers equations, the non-stationary Navier–Stokes equations, the non-stationary conduction–convection problems, and CVD equations based on POD method were presented (see [39–47]). Though a reduced-order FVE formulation based on POD method for parabolic equations has been presented (see [48]), to the best of our knowledge, there are no published results addressing the case that any reduced-order extrapolation algorithm (ROEA) based on the FVE method and the POD technique for finding solutions to the non-stationary Stokes equations is established or providing error estimates between classical SFVE solutions and the ROEA solutions or implementation for solving ROEA. In this paper, we extend the developments in [48], namely, apply the POD method to a classical SFVE formulation based on two local Gauss integrals and parameter-free for non-stationary Stokes equations, establish an ROEA with lower dimensions and sufficiently high accuracy based on SFVE method and POD technique for non-stationary Stokes equations, and provide the error estimates between the ROEA solutions and the classical SFVE solutions and the implementation for solving ROEA. Thus, we could provide some scientific theoretic basis and a new computing technology for service applications. Some numerical examples are provided for illustrating the fact that the results of numerical computation are consistent with theoretical conclusions. Moreover, it is shown that ROEA based on SFVE method and POD technique is feasible and efficient for finding numerical solutions for non-stationary Stokes equations. The plan of this paper is organized as follows. Section 2 is to derive first a semi-discrete formulation with respect to time for Problem I and the existence and error estimates of its solutions, and then start directly from semi-discrete formulation with respect to time to establish a classical fully discrete SFVE formulation based on two local Gauss integrals and parameter-free and to do error analysis. Thus, we could circumvent the semi-discrete SFVE formulation with respect to space variable such that our method is simpler and more convenient than existing methods (see, e.g., [13,17]). In Section 3, the optimal orthonormal POD bases are reconstructed from the first fewer given numerical solutions for the classical SFVE formulation on very short time span ½0; T 0  (T 0  T) and an ROEA with lower dimensional number based on the classical SFVE formulation and POD method for non-stationary Stokes equations on total time span ½0; T is developed. In Section 4, the error estimates between the classical SFVE solutions and the ROEA solutions and the implementation for solving ROEA are provided. In Section 5, some numerical examples are provided for illustrating that the errors between the ROEA approximate solutions and the classical SFVE solutions are consistent with previously obtained theoretical results, thus validating the feasibility and efficiency of ROEA based on the SFVE formulation and POD technique. Section 6 provides main conclusions and future tentative ideas.

978

Z. Luo / Applied Mathematics and Computation 247 (2014) 976–995

2. Semi-discrete formulation with respect to time and fully discrete SFVE formulation for non-stationary Stokes equations In this section, we first give a semi-discrete formulation with respect to time for Problem I and the existence and error estimates of its solutions, and then start directly from semi-discrete formulation with respect to time to establish a classical fully discrete SFVE formulation based on two local Gauss integrals and parameter-free and to do error analysis. Thus, we could circumvent the semi-discrete FVE formulation with respect to space variable such that our method is simpler and more convenient than existing other methods. 2.1. Semi-discrete formulation about time for non-stationary Stokes equations spaceso used in this context are standard (see [49]). Let U ¼ H10 ðXÞ2 n The Sobolev R q 2 L2 ðXÞ; X qdxdy ¼ 0 . Then, the variational formulation for Problem I can be stated as follows.

and

M ¼ L20 ðXÞ ¼

Problem II. Find ðuðtÞ; pðtÞÞ : ½0; T ! U  M such that

8 > < ðut ; v Þ þ maðu; v Þ þ bðv ; pÞ ¼ ðf ; v Þ 8v 2 U; bðu; qÞ ¼ 0 8q 2 M; > : uðx; y; 0Þ ¼ u0 ðx; yÞ; ðx; yÞ 2 X; @u1 @x @u2 @x

where aðu; v Þ ¼ ðru; rv Þ; ru ¼ ðru1 ; ru2 ÞT ¼ 2

2

2

22

ð2:1Þ @u1 @y @u2 @y

! ; bðv ; qÞ ¼ ðr  v ; qÞ, ð; Þ denotes the inner product in L2 ðXÞ or

L ðXÞ or L ðXÞ . Introducing the generalized bilinear form on U  M by

Bððu; pÞ; ðv ; qÞÞ ¼ maðu; v Þ þ bðv ; pÞ  bðu; qÞ;

ð2:2Þ

Problem II can be rewritten in the following compact form. Problem III. Find ðuðtÞ; pðtÞÞ : ½0; T ! U  M such that



ðut ; v Þ þ Bððu; pÞ; ðv ; qÞÞ ¼ ðf ; v Þ; uðx; y; 0Þ ¼ u0 ðx; yÞ; ðx; yÞ 2 X:

8ðv ; qÞ 2 U  M;

Put V ¼ fv 2 U; divv ¼ 0g and DðAÞ ¼ V \ H2 ðXÞ2 . f ðx; y; tÞ 2 L2 ð0; T; L2 ðXÞ2 Þ satisfy (see [4] or [19])

ku0 k2 þ

Z 0

T

ðkf k21 þ kf t k20 þ kf tt k20 Þdt

1=2

ð2:3Þ Let

the

initial

velocity

u0 2 DðAÞ

and

the

6 C;

body

force

ð2:4Þ

where C as well as occurrences in context is a positive constant independent of u; p, and f , but dependent on other datum of Problem I. The following estimates for bilinear term Bðð; Þ; ð; ÞÞ hold (see [4] or [19]):

jBððu; pÞ; ðu; pÞÞj ¼ mkruk20 ;

8ðu; pÞ 2 U  M;

ð2:5Þ

jBððu; pÞ; ðv ; qÞÞj 6 Cðkuk1 þ kpk0 Þðkv k1 þ kqk0 Þ; jBððu; pÞ; ðv ; qÞÞj P b0 ðkuk1 þ kpk0 Þ; kv k1 þ kqk0 ðv ;qÞ2UM sup

8ðu; pÞ;

ðv ; qÞ 2 U  M;

8ðu; pÞ 2 U  M;

ð2:6Þ ð2:7Þ

where b0 > 0 is a constant. The following result concerning existence, uniqueness, and regularity of a global strong solution to the non-stationary Stokes equations is well known (see [4] or [5]). Lemma 1. Assume that (2.4)–(2.7) hold. Then, for any given T > 0, there exists a unique solution ðu; pÞ to Problem II satisfying the following regularity

n o Z T sup kuðtÞk22 þ kpðtÞk21 þ kut ðtÞk20 þ kutt ðtÞk20 þ ðkruk20 þ jruj22 þ krut k22 þ kpt k21 þ krutt k20 Þdt 6 C:

0
ð2:8Þ

0

Let k denote the time step size and un be the semi-discrete approximation of uðtÞ at t n ¼ nk ðn ¼ 0; 1; . . . ; N ¼ ½T=kÞ. If the differential quotient ut in Problem III is approximated with the backward difference quotient @t un ¼ ðun  un1 Þ=k at time t ¼ tn , then the semi-discrete approximation scheme of Euler backward one step of Problem III about time reads as follows.

979

Z. Luo / Applied Mathematics and Computation 247 (2014) 976–995

Problem IV. Find ðun ; pn Þ 2 U  M ð1 6 n 6 NÞ such that

(

ð@t un ; v Þ þ Bððun ; pn Þ; ðv ; qÞÞ ¼ ðf ðt n Þ; v h Þ; u0 ¼ uðx; y; 0Þ;

8ðv ; qÞ 2 U  M;

ð2:9Þ

ðx; yÞ 2 X;

or, is equivalently read as follows. Problem V. Find ðun ; pn Þ 2 U  M ð1 6 n 6 NÞ such that

(

ðun ; v Þ þ kBððun ; pn Þ; ðv ; qÞÞ ¼ kðf ðt n Þ; v Þ þ ðun1 ; v Þ; 0

u ¼ uðx; y; 0Þ;

8ðv ; qÞ 2 U  M;

ð2:10Þ

ðx; yÞ 2 X:

There hold the following results for Problem IV. Theorem 2. Assume that (2.4)–(2.7) hold. Then, there exists a unique series of solution ðun ; pn Þ ð1 6 n 6 NÞ to Problem IV or Problem V satisfying the following error estimates

kuðtn Þ  un k0 þ k

n X 

 kpðt i Þ  pi k0 þ kuðt i Þ  ui k1 6 Ck;

n ¼ 1; 2; . . . ; N:

ð2:11Þ

i¼1

Proof. Due to assumptions (2.4)–(2.7), there exists a unique series of solution ðun ; pn Þ to Problem IV or Problem V for given ðun1 ; pn1 Þ ð1 6 n 6 NÞ (see [4] or [5]). Taking t ¼ tn in Problem II and using Taylor’s formula yield that

ðuðt n Þ; v Þ þ kBððuðtn Þ; pðtn ÞÞ; ðv ; qÞÞ ¼ kðf ðt n Þ; v Þ þ ðuðtn1 Þ; v Þ þ ðR; v Þ; uðt 0 Þ ¼ uðx; y; 0Þ;

8ðv ; qÞ 2 U  M;

ðx; yÞ 2 X;

ð2:12Þ ð2:13Þ

2

where R1 ¼ k utt ðnn Þ=2 and nn 2 ½tn1 ; tn . Put en ¼ uðt n Þ  un and sn ¼ pðtn Þ  pn (n ¼ 0; 1; 2; . . . ; N). Subtracting (2.10) from (2.12) and (2.13) yields that

ðen ; v Þ þ kBððen ; sn Þ; ðv ; qÞÞ ¼ ðen1 ; v Þ þ ðR; v Þ;

8ðv ; qÞ 2 U  M;

n ¼ 1; 2; . . . ; N;

ð2:14Þ

e0 ¼ 0: Taking

ð2:15Þ

v ¼ en and q ¼ sn in (2.14) and using Hölder inequality and Cauchy inequality yield that  k 1 km ken1 k20 þ ken k20 þ kutt ðnn Þk1 þ kren k20 ; 2 2m 2 3

ken k20 þ kmkren k20 6

n ¼ 1; 2; . . . ; N;

ð2:16Þ

namely,

ken k20 þ kmkren k20 6 ken1 k20 þ k

3 1

m kutt ðnn Þk1 ; n ¼ 1; 2; . . . ; N:

ð2:17Þ

Summing (2.17) from 1 to n and using (2.15) and (2.8) yield that

ken k20 þ km

n n X X 3 2 krei k20 6 k m1 kutt ðni Þk1 6 Ck ; i¼1

n ¼ 1; 2; . . . ; N:

ð2:18Þ

i¼1

Summing (2.14) from 1 to n and using (2.7), (2.15), (2.18), and (2.8) yield that

b0 k

n n X X 2 ðkei k1 þ ksi k0 Þ 6 ken k0 þ k kutt ðni Þk1 6 Ck: i¼1

ð2:19Þ

i¼1

Combining (2.18) with (2.19) yield (2.11), which completes the proof of Theorem 2.

h

2.2. Fully discrete SFVE formulation for non-stationary Stokes equations In order to get the numerical solution of SFVE for Problem II, it is necessary to introduce a FVE approximation for the spatial variables of Problem IV or Problem V. Firstly, let Ih ¼ fKg be a quasi-uniform triangulation of X with h ¼ max hK , where hK is the diameter of the triangle K 2 Ih (see [50] or [51]). In order to describe the SFVE formulation, we introduce a dual partition Ih based on Ih whose elements are called the control volumes. We construct the control volume in the same way as in [11]. Let zK ¼ ðxK ; yk Þ be the barycenter of K 2 Ih . We connect zK with line segments to the midpoints of the edges of K, thus partitioning K into three quadrilaterals K z

980

Z. Luo / Applied Mathematics and Computation 247 (2014) 976–995

S (z ¼ ðx; yÞ 2 Z h ðKÞ, where Z h ðKÞ are the vertices of K). Then with each vertex z 2 Z h ¼ K2Ih Z h ðKÞ we associate a control volume V z , which consists of the union of the sub-regions K z , sharing the vertex z. Finally, we obtain a group of control volumes covering the domain X, which is called a barycenter-type dual partition Ih of the triangulation Ih (see Fig. 1). We denote the set of interior vertices of Z h by Z h . The partition Ih is known as regular or quasi-uniform, if there exist two positive constants C 1 and C 2 , being independent of the spatial mesh size h and temporal mesh size k, such that 2

2

8V z 2 Ih :

C 1 h 6 mesðV z Þ 6 C 2 h ;

ð2:20Þ

The barycenter-type dual partition can be introduced for any FE triangulation Ih and leads to relatively simple calculations. Besides, if the FE triangulation Ih is quasi-uniform, then the dual partition Ih is also quasi-uniform (see [11] or [51]). The trial function space U h chosen as the linear element space related to Ih is the set of all the functions uh satisfying the following conditions:  Þ2 ; uh j ¼ 0; (i) uh 2 CðX @X (ii) uh jK 2 P 1 ðKÞ2 , namely uh is a linear vector function of xand yon each triangular element K 2 Ih , determined only by its values on the three vertexes. It is obvious that U h  U ¼ H10 ðXÞ2 . For u 2 U ¼ H10 ðXÞ2 , let Ph u be the interpolation projection of u onto the trial function space U h . By the interpolation theory of Sobolev spaces (see [2] or [11] or [51]), we have, if u 2 H2 ðXÞ2 , that

ju  Ph ujm 6 Ch

2m

m ¼ 0; 1;

juj2 ;

ð2:21Þ

where C in this context indicates a positive constant which is possibly different at different occurrences, being independent of the spatial mesh size h and temporal mesh size k. ~ h is chosen as the piecewise constant vector function space with respect to I but is zero vector on any The test space U h boundary dual element V z 2 Ih (V z \ @ X – ;), namely,

~h ¼ U

n

v h 2 L2 ðXÞ2 ; v h jV

z

2 P 0 ðV z Þ2 ðV z \ @ X ¼ ;Þ; v h jV z ¼ 0 ðV z \ @ X – ;Þ; 8V z 2 Ih

o

spanned by the following basis functions: for any point z 2 Z h ,

 /z ðx; yÞ ¼ For any vector

vh ¼

X

1; ðx; yÞ 2 V z ;

ð2:22Þ

0; elsewhere:

v h 2 U~ h ,

v h ðzÞ/z :

ð2:23Þ

z2Z h

~ h , namely, For w 2 U, let Ph w be the interpolation projection of w onto the test space U

Ph w ¼

X

wðzÞ/z :

ð2:24Þ

z2Z h

By the interpolation theory of Sobolev spaces (see [2] or [11] or [51]), we have

kw  Ph wk0 6 Chjwj1 : 0

ð2:25Þ 2

Moreover, due to U h  C ðXÞ \ [52]).

H10 ð

2

 h

XÞ , the interpolation projection P satisfies the following properties (see [11] or

Fig. 1. Left-hand chart is a triangle K partitioned into three sub-regions K z . Right-hand chart is a sample region with dotted lines indicating the corresponding control volume V z .

Z. Luo / Applied Mathematics and Computation 247 (2014) 976–995

Lemma 3. If

Z K

981

v h 2 U h , then

ðv h  Ph v h Þdxdy ¼ 0;

kv h  Ph v h kLr ðXÞ 6 Chkv h kW 1;r ðXÞ ;

K 2 Ih ;

1 6 r 6 1:

ð2:26Þ

~ h å U h . As in the case of nonconforming Though the trial function space U h satisfies U h  U like FE method, the test space U ~ h on the boundary of two neighboring elements. So the FE method, this is due to the loss of continuity of the functions in U bilinear forms aðu; v Þ and bðv ; pÞ must be revised accordingly. Note that

Z

X Z

Du  v dxdy ¼

X

V z 2Ih

Du  v dxdy ¼ 

Vz

X Z Vz

V z 2Ih

X Z

rurv dxdy þ

V z 2Ih

ðrunÞ  v ds

ð2:27Þ

@V z

and

Z

X Z

rp  v dxdy ¼ 

X

V z 2Ih

pdiv

v dxdy þ

Vz

X Z V z 2Ih

pv  nds;

ð2:28Þ

@V z

R where @Vz denotes the line integrals, with the counter clockwise direction, on the boundary @V z of the dual element; n ¼ ðn1 ; n2 ÞT is the unit outer normal vector to @V z . So bilinear forms aðu; v Þ and bðv ; pÞ of (2.1) are rewritten as

aðu; v Þ ¼

X Z

rurv dxdy  Vz

V z 2Ih

ðrunÞ  v ds

Z

ð2:29Þ

@V z

and

bðv ; pÞ ¼

X Z

pv  nds 

Z

:

ð2:30Þ

Vz

@V z

V z 2Ih

v dxdy

pdiv

~ h is the piecewise constant function space with the characteristic functions of the dual elements V z as the basis Since U functions, we have

aðu; v Þ ¼ 

X Z

ðrunÞ  v ds;

bðv ; pÞ ¼

@V z

V z 2Ih

X Z V z 2Ih

pv  nds;

~ h: 8v 2 U

ð2:31Þ

@V z

 Let M h ¼ qh 2 M; qh jK 2 P 1 ðKÞ; 8K 2 Ih . Then a fully discrete SFVE formulation for Problem II is written as follows. Problem VI. Find ðunh ; pnh Þ 2 U h  M h ð1 6 n 6 NÞ such that

8    n n  n  > < ð@ t uh ; Ph v h Þ þ mah ðuh ; Ph v h Þ þ bh ðPh v h ; ph Þ ¼ ðf ðtn Þ; Ph v h Þ; bðunh ; qh Þ  Dh ðpnh ; qh Þ ¼ 0; 8qh 2 Mh ; > : 0 uh ¼ Ph u0 ðx; yÞ; ðx; yÞ 2 X;

8v h 2 U h ; ð2:32Þ

where

ah ðunh ; Ph v h Þ ¼

X

v h ðzÞ  aðunh ; /z Þ;

ðunh ; /z Þ ¼  a

X

v h ðzÞ 

X Z

K2Ih

ð2:33Þ

Z

qh nds;

ð2:34Þ

@V z

z2Z h

Dh ðpnh ; qh Þ ¼ e

runh nds;

@V z

z2Z h

bh ðPh v h ; qh Þ ¼

Z

K;2

pnh qh dxdy 

Z K;1

pnh qh dxdy ;

ph ; qh 2 Mh ;

ð2:35Þ

R here e is a positive real number and K;i gðx; yÞdxdy indicate an appropriate Gauss integral over K that is exact for polynomials of degree i (i ¼ 1; 2Þ and gðx; yÞ ¼ ph qh is a polynomial of degree not more than i (i ¼ 1; 2Þ. Thus, for all test functions qh 2 M h , the trial function ph 2 M h must be piecewise constant when i ¼ 1. Consequently, we define the L2 -projection operator qh : L2 ðXÞ ! W h such that 8p 2 L2 ðXÞ satisfying

ðp; qh Þ ¼ ðqh p; qh Þ;

8qh 2 W h ;

ð2:36Þ

2

where W h  L ðXÞ denotes the piecewise constant space associated with Ih . The projection operator qh has the following properties (see [53] or [51]).

kqh pk0 6 Ckpk0 ;

8p 2 L2 ðXÞ;

kp  qh pk0 6 Chkpk1 ;

8p 2 H1 ðXÞ:

ð2:37Þ ð2:38Þ

982

Z. Luo / Applied Mathematics and Computation 247 (2014) 976–995

Now, using the definition of qh , we can rewrite the bilinear form Dh ð; Þ as follows:

Dh ðph ; qh Þ ¼ eðph  qh ph ; qh Þ ¼ eðph  qh ph ; qh  qh qh Þ:

ð2:39Þ

From [11] or [19], we have the following two lemmas. Lemma 4. The bilinear form ah ðuh ; Ph v h Þ is symmetric, bounded, and positive definite, namely,

að uh ; Ph v h Þ ¼ ah ðv h ; Ph uh Þ ¼ aðuh ; v h Þ;

8uh ;

v h 2 Uh ;

and there hold that

aðuh ; uh Þ P kruh k20 ;

aðuh ; v h Þ 6 kruh k0 krv h k0 ; The bilinear form bh ðP  h

bh ðP

 h

8 uh ; v h 2 U h :

v h ; qh Þ satisfies

v h ; qh Þ ¼ bðv h ; qh Þ; 8v h 2 Uh ;

qh 2 M h :

Lemma 5. There is the following statement:

 h Þ ¼ ðu  h ; Ph uh Þ; ðuh ; Ph u

8 uh ;

For any u 2 Hm ðXÞ2 ðm ¼ 0; 1Þ and  h

jðu; v h Þ  ðu; P Set jkuh kj0 ¼ ðuh ; P

v h 2 Uh ,

mþn

v h Þj 6 Ch

1=2  , h uh Þ

 h 2 Uh : u

kukm kv h kn ;

n ¼ 0; 1:

then jk  kj0 is equivalent to k  k0 on U h , namely, there exist two positive constants C 3 and C 4 such that

C 3 kuh k0 6 jkuh kj0 6 C 4 kuh k0 ;

8 uh 2 U h :

Put

Aððunh ; pnh Þ; ðv h ; qh ÞÞ ¼ maðunh ; v h Þ þ bðv h ; pnh Þ  bðunh ; qh Þ þ Dh ðpnh ; qh Þ:

ð2:40Þ

Then Problem VI can be rewritten as follows in a compact form. Problem VII. Find ðunh ; pnh Þ 2 U h  M h ð1 6 n 6 NÞ such that

(

ð@t unh ; Ph v h Þ þ Aððunh ; pnh Þ; ðv h ; qh ÞÞ ¼ ðf ðt n Þ; Ph v h Þ; u0h ¼ Ph u0 ðx; yÞ;

8ð v h ; q h Þ 2 U h  M h ;

ðx; yÞ 2 X;

ð2:41Þ

or, is equivalently read as follows. Problem VIII. Find ðunh ; pnh Þ 2 U h  M h ð1 6 n 6 NÞ such that

(

 ðunh ; Ph v h Þ þ kAððunh ; pnh Þ; ðv h ; qh ÞÞ ¼ kðf ðtn Þ; Ph v h Þ þ ðun1 h ; Ph v h Þ;

u0h

¼ uh ðx; y; 0Þ ¼ Ph u0 ðx; yÞ;

8ðv h ; qh Þ 2 U h  Mh ;

ðx; yÞ 2 X:

ð2:42Þ

Problem VIII is also referred to as Euler backward one step fully discrete SFVE formulation. From Theorem 4.1 in [19], there holds the following Lemma 6. Lemma 6. The bilinear form Að; Þ on ðU h ; M h Þ  ðU h ; M h Þ holds the following continuity and weak coercivity

Aððuh ; ph Þ; ðv h ; qh ÞÞ 6 Cðkuh k1 þ kph k0 Þðkv h k1 þ kqh k0 Þ; sup ðv h ;qh Þ2U h M h

Aððuh ; ph Þ; ðv h ; qh ÞÞ P bðkuh k1 þ kph k0 Þ; kv h k1 þ kqh k0

8ðuh ; ph Þ; ðv h ; qh Þ 2 U h  Mh ;

8ðuh ; ph Þ 2 U h  Mh ;

ð2:43Þ ð2:44Þ

where b is independent of h. For given solutions ðun ; pn Þ 2 U  M (n ¼ 1; 2; . . . ; N) of Problem V, define the stabilized Stokes projection ðSh ; Q h Þ : U  M ! U h  M h such that Sh u0 ¼ Ph u0 and

ðSh un  Sh un1 ; v h Þ þ kAððSh un ; Q h pn Þ;

8ðv h ; qh Þ 2 U h  Mh ;

n ¼ 1; 2; . . . ; N:

ðv h ; qh ÞÞ ¼ ðun  un1 ; v h Þ þ kBððun ; pn Þ; ðv h ; qh ÞÞ; ð2:45Þ

In fact, (2.45) is the error equation between Problem V and its standard stabilized FE formulation with solutions ðSh un ; Q h pn Þ. Thus, from standard stabilized FE method (see, e.g., [2] or [4] or [51]), there hold the following results.

983

Z. Luo / Applied Mathematics and Computation 247 (2014) 976–995

Lemma 7. The stabilized Stokes projection ðSh ; Q h Þ : U  M ! U h  M h satisfies

kun  Sh un k1 þ kpn  Q h pn k0 6 Cðkun k1 þ kpn k0 Þ; 2

2

n ¼ 1; 2; . . . ; N:

ð2:46Þ

1

If ðun ; pn Þ 2 H ðXÞ  H ðXÞ, there holds 2

2

kun  Sh un k0 þ hðkun  Sh un k1 þ kpn  Q h pn k0 Þ 6 Ch ðkun k2 þ kpn k1 Þ 6 Ch :

ð2:47Þ

Theorem 8. Under the hypotheses of Theorem 2, there exists a unique series of solutions ðunh ; pnh Þ ðn ¼ 1; 2; . . . ; NÞ of fully discrete SFVE formulation, namely, Problem VI satisfying

kunh k20

! n n X X 2 i 2 i 2 þ k jruh j0 6 C ku0 k0 þ k kf k1 ; i¼1

2

kkpnh k0

ð2:48Þ

i¼1

n X 6 C 4ku0 k0 þ k kf i k21

!1=2 3 5;

ð2:49Þ

i¼1

where C is the constant independent to h and k, which shows that the series of solutions of Problem VI is stable and continuously dependence on given the initial value vector u0 and body force vector f ðx; y; tÞ.   n Proof. Define Fðv h Þ ¼ ðun1 h ; Ph v h Þ þ kðf ; Ph v h Þ and

 n n n n n n n ~ Aððu h ; ph Þ; ðv h ; qh ÞÞ ¼ ðuh ; Ph v h Þ þ kmaðuh ; v h Þ þ kbðv h ; ph Þ  kbðuh ; qh Þ þ kDh ðph ; qh Þ:

ð2:50Þ

Then, by Lemmas 5 and 6, there holds

~ ðunh ; Ph unh Þ þ kAððuh ; ph Þ; ðv h ; qh ÞÞ Aððu h ; ph Þ; ðv h ; qh ÞÞ P sup kv h k1 þ kqh k0 kv h k1 þ kqh k0 ðv h ;qh Þ2U h M h ðv h ;qh Þ2U h Mh sup

sup

P

ðv h ;qh Þ2U h Mh

kAððuh ; ph Þ; ðv h ; qh ÞÞ P bkðkuh k1 þ kph k0 Þ; kv h k1 þ kqh k0

8ðuh ; ph Þ 2 U h  Mh :

ð2:51Þ

~ Þ is a bounded bilinear function and FðÞ is a bounded linear function for given un1 . Thus, Problem VII, Obvious, Að; h namely, Problem VI has a unique series of solutions ðunh ; pnh Þ ðn ¼ 1; 2; . . . ; NÞ from Lax–Milgram Theorem (see, e.g., [50] or [51]). Taking v h ¼ unh in the first equation of Problem VI and qh ¼ pnh in the second equation of Problem IV, by using Lemma 5, Hölder inequality, and Cauchy inequality, we obtain  n  n n kjunh jk20 þ kmjrunh j20 þ kekpnh  qh pnh k20 ¼ ðun1 h ; Ph uh Þ þ kðf ; Ph uh Þ

6

1 1 Ck km kjunh jk20 þ kjunh jk20 þ kf n k21 þ kunh k21 : 2 2 2m 2

ð2:52Þ

Further, we obtain 2 kjunh jk20 þ kmkrunh k20 6 Ckm1 kf n k21 þ kjun1 h jk0 :

ð2:53Þ

Summing (2.53) from 1 to n yields that

kjunh jk20 þ km

n n n X X X juih j21 6 kju0h jk20 þ Ckm1 kf i k21 6 Ckju0 jk20 þ Ckm1 kf i k21 ; i¼1

i¼1

ð2:54Þ

i¼1

which yields (2.48) by using Lemma 5. By the first equation of Problem VIII and using (2.51), Hölder inequality, and (2.48), we obtain  ~ kðf n ; Ph v h Þ þ ðun1 Aððu h ; ph Þ; ðv h ; qh ÞÞ n h ; Ph v h Þ 6 sup 6 Cðkun1 h k0 þ kkf k1 Þ kv h k1 þ kqh k0 kv h k1 v h 2U h ðv h ;qh Þ2U h Mh 2 !1=2 3 n X i 2 4 5; 6 C ku0 k0 þ k kf k1

bkkpnh k0 6

sup

ð2:55Þ

i¼1

which yields (2.49) and completes the proof of Theorem 8.

h

In order to derive the error estimates of the solutions to Problem VIII, it is necessary to introduce the following discrete Gronwall Lemma (see, e.g., [4] or [51]).

984

Z. Luo / Applied Mathematics and Computation 247 (2014) 976–995

Lemma 9. (discrete Gronwall Lemma). If fan g and fbn g are two positive sequences, fcn g is non-decreasing positive sequence, and they satisfy n1 X an þ bn 6 cn þ k ai ;

k > 0;

a0 þ b0 6 c0 ;

ð2:56Þ

i¼0

then

an þ bn 6 cn expðnkÞ;

n P 0:

ð2:57Þ

There holds yet the following main results for Problem VIII. Theorem 10. Under the hypotheses of Theorem 2, let ðu; pÞ be the solution for Problem III and ðunh ; pnh Þ the solution of fully discrete 2

SFVE scheme, namely, Problem VII, respectively. Then, if k ¼ OðhÞ; u0 2 H2 ðXÞ2 , and f 2 L1 ðð0; TÞ; H1 ðXÞÞ , there holds that n   X   2 k kpðt i Þ  pih k0 þ kuðt i Þ  uih k1 þ kuðtn Þ  unh k0 6 C h þ k ;

n ¼ 1; 2; . . . ; N:

ð2:58Þ

i¼1

Proof. Taking equation:

v ¼ vh

and q ¼ qh in Problem III, subtracting with Problem IV, using Lemmas 4 and 3, we obtain the error

8 n ðu  unh ; v h Þ þ ðunh  Ph unh ; v h  Ph v h Þ þ kaðun  unh ; v h Þ þ kbðv h ; pn  pnh Þ ¼ kðf n  Ph f n ; v h  Ph v h Þ > > >  < n1  Ph un1 8v h 2 U h ; n ¼ 1; 2; . . . ; N; þðun1  un1 h ; v h Þ þ ðuh h ; v h  Ph v h Þ; n n n n > bðu  u ; q Þ  e ðp  q p ; q  q q Þ ¼ 0; 8 q 2 Mh ; n ¼ 1; 2; . . . ; N; > h h h h h h h h h > : 0 0 u  uh ¼ u0 ðx; yÞ  Ph u0 ðx; yÞ; ðx; yÞ 2 X:

ð2:59Þ

By using the error Eq. (2.59), (2.45), and Lemmas 4 and 5, we obtain

kSh un  unh k20 þ kmjSh un  unh j21 ¼ ðSh un  un ; Sh un  unh Þ þ kaðSh un  un ; Sh un  unh Þ þ ðun  unh ; Sh un  unh Þ þ kaðun  unh ; Sh un  unh Þ ¼ ðSh un1  un1 ; Sh un  unh Þ  kbðSh un  unh ; Q h pn  pn Þ  n n n n  kbðSh un  unh ; pn  pnh Þ þ ðun1  Ph un1 h h ; Sh u  uh  Ph ðSh u  uh ÞÞ

 ðunh  Ph unh ; Sh un  unh  Ph ðSh un  unh ÞÞ þ kðf n  Ph f n ; Sh un  unh  Ph ðSh un  unh ÞÞ  n  n n n n n n n þ ðun1  un1 h ; Sh u  uh Þ ¼ kðf  Ph f ; Sh u  uh  Ph ðSh u  uh ÞÞ  n n n n  Ph ðunh  un1  kbðSh un  unh ; Q h pn  pnh Þ  ðunh  un1 h h Þ; Sh u  uh  Ph ðSh u  uh ÞÞ:

ð2:60Þ Using Lemma 5, Hölder inequality, and Cauchy inequality, if h ¼ OðkÞ, we obtain 2

4

jkðf n  Ph f n ; Sh un  unh  Ph ðSh un  unh ÞÞj 6 Ckh kf n k1 krðSh un  unh Þk0 6 Ckh kf n k21 þ

mk 4

jSh un  unh j21 ;

ð2:61Þ

With inverse error estimate and Taylor’s formula, there holds that 2

 n n n n n n1 n n 2 jðunh  un1  Ph ðunh  un1 h h Þ; Sh u  uh  Ph ðSh u  uh ÞÞj 6 Ch kuh  uh k1 krðSh u  uh Þk0 3

3

3

2 6 Ch krðunh  Sh un Þk20 þ Ch krðSh un  Sh un1 Þk20 þ Ch krðSh un1  un1 h Þk0 þ 2 3

2 6 Chkunh  Sh un k20 þ Ck h kut k2L1 ðH1 Þ þ ChkSh un1  un1 h k0 þ

km jSh un  unh j21 8

km jSh un  unh j21 : 4

ð2:62Þ

Noting that bðSh un  un ; qh Þ ¼ keðQ h pn  qh ðQ h pn Þ; qh  qh ðqh ÞÞ, by the property of operator qh and the second equation of (2.59), we have

bðSh un  unh ; Q h pn  pnh Þ ¼ bðSh un  un ; Q h pn  pnh Þ  bðun  unh ; Q h pn  pnh Þ ¼ eðQ h pn  pnh  qh ðQ h pn  pnh Þ; Q h pn  pnh  qh ðQ h pn  pnh ÞÞ  eðQ h pn1  pn1 h n n n n  qh ðQ h pn1  pn1 h Þ; Q h p  ph  qh ðQ h p  ph ÞÞ

e

e

2 6  kQ h pn  pnh  qh ðQ h pn  pnh Þk20 þ kQ h pn1  pn1  qh ðQ h pn1  pn1 h h Þk0 : 2 2

ð2:63Þ

Combining (2.61), (2.62), (2.63) and (2.60), we obtain 2 kSh un  unh k20 þ kmjSh un  unh j21 þ ekQ h pn  pnh  qh ðQ h pn  pnh Þk20  ekQ h pn1  pn1  qh ðQ h pn1  pn1 h h Þk0 4

2 3

2 2 n1 n n 2 6 Ckh kf n k21 þ Ck h kut k2L1 ðH1 Þ þ kSh un1  un1  un1 h k0 þ ChkSh u h k0 þ Chkuh  Sh u k0 :

ð2:64Þ

985

Z. Luo / Applied Mathematics and Computation 247 (2014) 976–995

Summing (2.64) from 1 to n and noting that Sh u0  u0h ¼ 0, if h is sufficiently small such that Ch 1=2 in (2.64) and p0h ¼ p0 ¼ 0 (or p0h ¼ Q h p0 ), we obtain

kSh un  unh k20 þ km

n n n1 X X X 4 2 3 jSh ui  uih j21 6 Ckh kf i k21 þ Cnk h kut k2L1 ðH1 Þ þ Ch kSh ui  uih k20 : i¼1

i¼1

ð2:65Þ

i¼0

By Gronwall Lemma 9, we obtain n

kSh u 

unh k20

" # n n X X 4 2 3 4 2 i i 2 i 2 þ km jSh u  uh j1 6 C kh kf k1 þ nk h kut kL1 ðH1 Þ expðcnhÞ 6 Ch : i¼1

Noting that

Pn

2 i¼1 ai

P

ð2:66Þ

i¼1

Pn

i¼1 ai

2

=n, by using (2.21) and triangle inequality, we obtain

n X 2 kun  unh k0 þ k jui  uih j1 6 Ch :

ð2:67Þ

i¼1

From (2.51), (2.45), error Eq. (2.59), (2.67), and Lemmas 5 and 7, we have

bkkQ h pn  pnh k0 6 ¼

~ h un  uh ; Q h  pn p Þ; ðv h ; q ÞÞ AððS h h kv h k1 þ kqh k0 ðv h ;qh Þ2U h M h sup

 kðf n ; v h  Ph v h Þ  ðSh un ; v h  Ph v h Þ þ ðun1 ; v h  Ph v h Þ þ ðun1  un1 h ; Ph v h Þ kv h k1 þ kqh k0 ðv h ;qh Þ2U h Mh

sup 2

2

2

2

2

6 Cðkh kf k1 þ h kun k1 þ h kun1 k1 þ h Þ 6 Ch :

ð2:68Þ

Applying triangle inequality and Lemma 7 to (2.68) yields that 2

kkpn  pnh k0 6 Ch :

ð2:69Þ

Combining (2.67) and (2.69) with Theorem 2 yields (2.58), which completes Theorem 10.

h

If body force vector f ðx; y; tÞ, initial condition vector u0 ðx; yÞ, the triangulation parameter h, the time step increment k, trial function spaces U h and M h , and parameter-free e are given, by solving Problem VI, we can obtain a solution ensemble fðunh ; pnh ÞgNn¼1 for Problem VI. And then we choose the first L (in general, L  N, for example, L ¼ 20, N ¼ 2000 in Section 5) instantaneous solutions ðuih ðx; yÞ; pih ðx; yÞÞ (1 6 i 6 L) (which are useful and of interest for us) from N instantaneous solutions fðunh ; pnh ÞgNn¼1 for Problem VI, which are referred to as snapshots. Remark 1. When one computes actual problems, he may obtain the ensemble of snapshots from physical system trajectories by drawing samples from experiments with interpolation (or data assimilation), then reconstruct the POD optimal basis for the ensemble of snapshots by using the following POD method, and finally the trial function spaces U h  M h are replaced with the subspaces generated with POD basis in order to derive a reduced-order dynamical system with lower dimensions. Thus, the future change of physical phenomenon representing actual problems can be quickly simulated, which is a result of major importance for real-life applications.

3. Generation of POD basis and an ROEA based on SFVE method and POD technique For ðuih ðx; yÞ; pih ðx; yÞÞ (1 6 i 6 L) in Section 2, let U i ðx; yÞ ¼ ðuih ðx; yÞ; pih ðx; yÞÞ ð1 6 i 6 LÞ and

V ¼ spanfU 1 ; U 2 ; . . . ; U L g

ð3:1Þ fU i gLi¼1

and refer to V as the space generating by the snapshots at least one of which is assumed to be non-zero. Let fwj glj¼1 denote an orthonormal basis of V with l ¼ dim V. Then each member of the ensemble can expressed as

Ui ¼

l X ðU i ; wj ÞU^ wj ;

i ¼ 1; 2; . . . ; L;

ð3:2Þ

j¼1

^ ¼ U  M; ðU i ; w Þ ¼ ðrui ; rw Þ þ ðpi ; w Þ, and w and w are orthogonal bases corresponding to u and p, where U ^ j U uj pj uj pj h h respectively. Definition 1. The POD method consists in finding the orthonormal basis wj ði ¼ 1; 2; . . . ; dÞ such that for every d ð1 6 d 6 lÞ the mean square error between the elements U i ð1 6 i 6 LÞ and corresponding dth partial sum of (3.2) is minimized on average

2   L  d X 1X   ðU i ; wj ÞU^ wj  U i   fwj gdj¼1 L i¼1  j¼1 min

^ U

ð3:3Þ

986

Z. Luo / Applied Mathematics and Computation 247 (2014) 976–995

subject to

ðwi ; wj ÞU^ ¼ dij

1 6 i 6 d; 1 6 j 6 i;

ð3:4Þ

where kU i k2U^ ¼ kruih k20 þ kpih k20 . A sequence of solution fwj gdj¼1 of (3.3) and (3.4) is known as a POD basis of rank d. We introduce the correlation matrix A ¼ ðAij ÞLL 2 RLL corresponding to the snapshots fU i gLi¼1 by

Aij ¼

1 ðU i ; U j ÞU^ : L

ð3:5Þ

The matrix A is positive semi-definite and has rank l. Thus, the solution of (3.3) and (3.4) can be found, in addition, the following results hold (see [30] or [38]). Proposition 11. Let k1 P k2 P . . . P kl > 0 denote the positive eigenvalues of A and eigenvectors. Then a POD basis of rank d 6 l is given by L 1 X wi ¼ pffiffiffiffiffiffi ðv i Þj U j ; Lki j¼1

v 1; v 2; . . . ; v l

the associated orthonormal

1 6 i 6 d 6 l;

ð3:6Þ

where ðv i Þj denote the jth component of the eigenvector

v i . Furthermore, the following error formula holds

2   L  d l X X 1X   ðU i ; wj ÞU^ wj  ¼ kj : U i   L i¼1  j¼1 j¼dþ1 ^

ð3:7Þ

U

Let V d ¼ spanfw1 ; w2 ; . . . ; wd g and U d  M d ¼ V d with U d ¼spanfwu1 ; wu2 ; . . . ; wud g  U and M d ¼ span fwp1 ; wp2 ; . . . ; wpd g  M. Define the Ritz-projection P d : U h ! U d and L2 -projection qd : M h ! M d g as follows, respectively:

ðrPd uh ; rwd Þ ¼ ðruh ; rwd Þ; ðqd p; qd Þ ¼ ðph ; qd Þ;

8uh 2 U h ; 8wd 2 U d ;

ð3:8Þ

8ph 2 Mh ; 8qd 2 Md :

ð3:9Þ

Therefore, they are well-defined and bounded, namely,

krðPd uh Þk0 6 kruh k0 ;

8uh 2 U h ;

kqd ph k0 6 kph k0 ;

8ph 2 Mh :

ð3:10Þ h

d

h

d

Then from functional analysis theory (see [54]), there are two extensions P : U ! U h of P such that P jUh ¼ P : U h ! U d and qh : M ! M h of qd such that qh jMh ¼ qd : M h ! M d , respectively, denoted by

ðrPh u; rwh Þ ¼ ðru; rwh Þ; ðqh p; qh Þ ¼ ðp; qh Þ;

8u 2 U;

8p 2 M; h

8wh 2 U h ;

ð3:11Þ

8q h 2 M h :

ð3:12Þ

h

Thus, the linear operators P and q are also well-defined and bounded

krðPh uÞk0 6 kruk0 ;

8u 2 U;

kqh pk0 6 kpk0 ;

8p 2 M:

ð3:13Þ

Lemma 12. For every dð1 6 d 6 lÞ, the projection operators Pd and qd respectively satisfy that L l X 1X krðuih  Pd uih Þk20 6 kj ; L i¼1 j¼dþ1

ð3:14Þ

L l X 1X 2 kuih  Pd uih k20 6 Ch kj ; L i¼1 j¼dþ1

ð3:15Þ

L l X 1X kpih  qd pih k20 6 kj ; L i¼1 j¼dþ1

ð3:16Þ

where ðuih ; pih Þ 2 U h  M h ð1 6 i 6 LÞ are the first L solutions of Problem VI. Proof. For any u 2 U, using (3.11), we obtain

krðu  Ph uÞk20 ¼ ðrðu  Ph uÞ; rðu  Ph uÞÞ ¼ ðrðu  Ph uÞ; rðu  wh ÞÞ 6 krðu  Ph uÞk0 krðu  wh Þk0 ;

8wh 2 U h : ð3:17Þ

987

Z. Luo / Applied Mathematics and Computation 247 (2014) 976–995

Therefore, we get that

krðu  Ph uÞk0 6 krðu  wh Þk0 ;

8w h 2 U h :

Pd

uih ,

i ~ wuj j¼1 ðuh ; wuj ÞU

ð3:18Þ d

If u ¼ then taking wh ¼ 2 U  U h (where wuj is the component of wj corresponding to u) in (3.18), from (3.7), we obtain (3.14). In order to prove (3.12), we consider the following variational problem

ðrw; ruÞ ¼ ðu  Ph u; uÞ;

8u 2 U:

ð3:19Þ

Since u  P h u 2 U, Eq. (3.19) has a unique solution w 2 H10 ðXÞ2 \ H2 ðXÞ2 such that kwk2 6 Cku  P h uk0 . Taking u ¼ u  Ph u in (3.19) and using (3.11), we get that

ku  Ph uk20 ¼ ðrw; rðu  P h uÞÞ ¼ ðrðw  wh Þ; rðu  P h uÞÞ 6 krðw  wh Þk0 krðu  Ph uÞk0 ;

8wh 2 U h :

ð3:20Þ

Taking wh ¼ Ph w as interpolation function of w in U h and using (2.21) or interpolation theory (see, e.g., [51] or [50]) and (3.20), we obtain

ku  Ph uk20 6 Chkwk2 krðu  Ph uÞk0 6 Chku  Ph uk0 krðu  Ph uÞk0 :

ð3:21Þ

Therefore, we get

ku  Ph uk0 6 Chkrðu  Ph uÞk0 :

ð3:22Þ

Thus, if u ¼ uih ; Ph is restricted to Ritz-projection from U h to U d such that Ph jU h ¼ P d : U h ! U d , namely, Ph uih ¼ P d uih 2 U d , then from (3.22) and (3.14) we derive (3.15). Using Hölder inequality and (3.9) yields

kpih  qd pih k20 ¼ ðpih  qd pih ; pih  qd pih Þ ¼ ðpih  qd pih ; pih  qd Þ 6 kpih  qd pih k0 kpih  qd k0 ;

8q d 2 M d :

Consequently,

kpih  qd pih k0 6 kpih  qd k0 ; Taking qd ¼

Pd

i j¼1 ðph ; wpj Þ0 wpj

8q d 2 M d :

ð3:23Þ

(where wpj is the component of wj corresponding to p) in (3.23), from (3.7) we can obtain

(3.16), which completes the proof of Lemma 12.

h

From (3.11), (3.12), and standard FE method (see, e.g., [2] or [51]), we immediately obtain the following corollary. Corollary 13. If ðu; pÞ 2 Hm ðXÞ2  Hm ðXÞðm ¼ 1; 2Þ, there are the following results for the Ritz-projection P h and L2 -projection qh defined by (3.11) and (3.12), respectively, m

ku  Ph uk0 þ hkrðu  Ph uÞk0 6 Ch jujm ;

m

kp  qh pk0 6 Ch jpjm :

ð3:24Þ

Thus, by using U d and M d , an ROEA based on SFVE method and POD technique for non-stationary Stokes equations can be established as follows. Problem IX. Find ðund ; pnd Þ 2 U d  M d ð1 6 n 6 NÞ such that

8 d d X X > > < un ¼ P d un ¼ ðrunh ; rwui Þwui ; pnd ¼ qd pnh ¼ ðpnh ; wpi Þwpi ; n ¼ 1; 2; . . . ; L; d h i¼1 i¼1 > > : n   8ðv d ; qd Þ 2 U d  Md ; ðud ; Ph v d Þ þ kAððund ; pnd Þ; ðv d ; qd ÞÞ ¼ kðf ðt n Þ; Ph v d Þ þ ðun1 d ; Ph v d Þ;

n ¼ L; L þ 1; . . . ; N; ð3:25Þ

where

ðunh ; pnh Þ

(n ¼ 1; 2; . . . ; L  N) are the first L solutions to Problem VIII.

Remark 2. Since Ih is a uniformly regular triangulation and U h and M h are respectively the spaces of piecewise linear vector function and piecewise linear function, the total degrees of freedom for Problem VI, namely, the number of unknown quantities is 3N h (where N h is the number of vertices of triangles in Ih , see [50] or [51]), while the number of total degrees of freedom for Problem IX is 3d ðd  l 6 L  NÞ. For scientific engineering problem, the number of vertices of triangles in Ih is more than ten of thousands or even can exceed more than a hundred million, while d is only the number of few maximal eigenvalues which is chosen L snapshots from the N snapshots so that it is very small (for

988

Z. Luo / Applied Mathematics and Computation 247 (2014) 976–995

example, in Section 5, d ¼ 6, while N h P 56  106 ). Therefore, Problem IX is an ROEA based on SFVE method and POD technique for non-stationary Stokes equations. Especially, it uses given the first fewer L solutions of Problem VI to extrapolate n > L all other solutions and has no re-computation, it is completely different from existing reduced-order formulation based on POD technique (see, e.g., [31,32,38–40,44,45,48]). Moreover, since the development and change of future fluid flow are closely related to previous results and ROEA based on POD method for non-stationary Stokes equations is established with existing results as snapshots to structure POD basis, it may truly capture laws of change of future fluid flow.

4. Error estimates of solutions for Problem IX and implementation of its algorithm 4.1. Error estimates of solutions for Problem IX In this subsection, we recur to the classical SFVE method to derive the error estimates for Problem IX. To this end, it is still necessary to introduce a preparative lemma. Lemma 14. It holds that 1

1

1

1

~ 2 2 2 2 jAððu h ; ph Þ; ðv h ; qh ÞÞj 6 Cðkuh k0 þ k kruh k0 þ k kph k0 Þðkv h k0 þ k krv h k0 þ k kqh k0 Þ;

8ðuh ; ph Þ;

ðv h ; qh Þ 2 U h  M h ; ð4:1Þ

ðud ; pd ÞÞ ¼ kjud jk20 þ kmkrud k20 þ kekpd k20  kekqh pd k20  0;

~ Aððu d ; pd Þ;

8ðud ; pd Þ 2 U d  Md

ð4:2Þ

and if and only if ud ¼ 0 and pd ¼ 0, the second equality sign in (4.2) holds. Proof. It follows from (2.37) and Lemmas 4 and 5 that  ~ jAððu h ; ph Þ; ðv h ; qh ÞÞ 6 jðuh ; Ph v h Þj þ kmjah ðv h ; uh Þj þ kjbh ðv h ; ph Þj þ kjbðuh ; qh Þj þ kjDðph ; qh Þj 1

1

1

1

6 Cðkuh k0 þ k2 kruh k0 þ k2 kph k0 Þðkv h k0 þ k2 krv h k0 þ k2 kqh k0 Þ;

8ðuh ; ph Þ;

ðv h ; qh Þ 2 U h  M h ; ð4:3Þ

namely, the continuity result (4.1) holds. For each pd 2 M d  M h , it follows from (2.36) and (2.39) that

ðpd  qh pd ; pd Þ ¼ ðpd  qh pd ; pd  qh pd Þ ¼ ðpd  qh pd ; pd þ qh pd Þ:

ð4:4Þ

Further, there holds that

ðqh pd ; pd Þ ¼ kqh pd k20 ;

ðpd  qh pd ; pd  qh pd Þ ¼ kpd k20  kqh pd k20  0:

ð4:5Þ

Thus, it follows from (2.37) and Lemmas 4 and 5 that  ~ Aððu d ; pd Þ; ðud ; pd ÞÞ ¼ ðud ; Ph ud Þ þ kmaðud ; ud Þ þ eðpd  qh pd ; pd  qh pd Þ

¼ kjud jk20 þ kmkrud k20 þ kekpd k20  kekqh pd k20 ; which completes the proof of Lemma 14.

8ðud ; pd Þ 2 U d  Md ;

ð4:6Þ

h

In the following, we derive the following existence, stabilization, and convergence of the series of solutions to Problem IX. Theorem 15. Under hypotheses of Theorem 10, Problem IX has a unique series of solutions ðund ; pnd Þ 2 U d  M d such that

kund k20 þ k

! n n X X jruid j20 6 C ku0 k20 þ k kf i k21 ; i¼1

2 kkpnd k0

6 C 4ku0 k0 þ k

n ¼ 1; 2; . . . ; L;

ð4:7Þ

i¼1

n X kf i k21

!1=2 3 5;

n ¼ 1; 2; . . . ; L;

ð4:8Þ

i¼1

kund k20

! n n X X 2 2 i 2 i 2 i 2 þ k ½kud k1 þ kkpd k0  kkqh pd k0  6 C ku0 k0 þ k kf ðt i Þk1 ; i¼1

L þ 1 6 n 6 N;

ð4:9Þ

i¼1

which show that the series of solutions ðund ; pnd Þ (n ¼ 1; 2; . . . ; N) to Problem IX is stabilized. If k ¼ OðhÞ, then the following error estimates hold

989

Z. Luo / Applied Mathematics and Computation 247 (2014) 976–995

kunh



und k0

þk

n X 

kr

ðuih



uid Þk0

þ

kpih



pid k0



6 CkL

i¼1

kunh



und k0

n X

þk

l X

!1=2 kj

;

n ¼ 1; 2; . . . ; L;

ð4:10Þ

j¼dþ1

ðuih

½kr



uid Þk0

þ

kpih



pid k0 

6 Ck

1=2

hðn  LÞ þ CkL

i¼Lþ1

l X

!1=2 kj

;

L þ 1 6 n 6 N:

ð4:11Þ

j¼dþ1

Proof. Due to Lemma 14, there exists a unique series of solutions ðund ; pnd Þ (n ¼ 1; 2; . . . ; N) to Problem IX. Noting that it follows from standard FE method (see, e.g., [51]) and krP h uk0 6 Ckruk0 (8u 2 H10 ðXÞ2 ) that kPh uk0 6 kuk0 . Thus, for n ¼ 1; 2; . . . ; L, by (3.10) and Theorem 8, there hold that

kund k20

! n n n X X X 2 i 2 n 2 i 2 i 2 þ k krud k0 6 Ckuh k0 þ k kruh k0 6 C ku0 k0 þ k kf k1 ; i¼1

i¼1

ð4:12Þ

i¼1

that is just (4.7). But by (3.10) and Theorem 8, there yields (4.8) immediately. Taking v d ¼ und and qd ¼ pnd in Problem IX, it follows from Lemma 5 and (4.5) that  n kjund jk20 þ kmkrund k20 þ kekpd k20  kekqh pnd k20 ¼ kðf n ; Ph und Þ þ ðun1 d ; Ph ud Þ;

L þ 1 6 n 6 N:

ð4:13Þ

By using Hölder inequality and Cauchy inequality, we get that

kjund jk20 þ kmkrund k20 þ kekpd k20  kekqh pnd k20 6

Ckm1 n 2 km 1 1 2 n 2 kf k1 þ krud k20 þ kjun1 d jk0 þ kjud jk0 : 2 2 2 2

ð4:14Þ

Moreover, we get that 2 kjund jk20 þ kmkrund k20 þ kekpd k20  kekqh pnd k20 6 Ckm1 kf n k21 þ kjun1 d jk0 ;

L þ 1 6 n 6 N:

ð4:15Þ

Summing (4.15) from L þ 1 to n and using Lemma 5 and (4.7) yield (4.9). For n ¼ 1; 2; . . . ; L, if h 6 k or h ¼ OðkÞ, (3.22) and Lemma 12, there hold that

kunh  und k0 þ k

n X   krðuih  uid Þk0 þ kpih  pid k0 i¼1

¼ kunh  Pd unh k0 þ k

n h X

n h i i X krðuih  Pd ui Þk0 þ kpih  qd pih k0 6 Chkrðunh  Pd unh Þk0 þ k krðuih  Pd ui Þk0 þ kpih  qd pih k0

i¼1

6 Ck

n h X

kr

ðuih  P d ui Þk0 þ kpih

i¼1

6 CkL

l X

i¼1

q

!1=2 kj

d i ph k0

i

n h i 1X 6 CkL krðuih  Pd ui Þk20 þ kpih  qd pih k20 L i¼1

!1=2

; n ¼ 1;2;. .. ;L;

ð4:16Þ

j¼dþ1

that is just (4.10). Since U d  M d  U h  M h , subtracting Problem IX from Problem VIII taking ðv h ; qh Þ ¼ ðv d ; qd Þ 2 U d  M d yields from Lemma 4 that

eðpnh  pnd  qh ðpnh  pnd Þ; qd Þ þ ðdivðunh  und Þ; qd Þ ¼ 0; 8qd 2 Md ; L þ 1 6 n 6 N;  ðunh  und ; Ph v d Þ þ kmðrðunh  und Þ; rPh v d Þ  kðdivPh v d ; pnh  pnd Þ ¼ ðun1  un1 h d ; Ph v d Þ;

ð4:17Þ

8v d 2 U d ;

L þ 1 6 n 6 N: ð4:18Þ

We have from Lemmas 4 and 5, (4.18), (3.11), (4.17), (3.22), Hölder inequality, and Cauchy inequality that

kjunh  und jk20 þ kmkrðunh  und Þk20 þ kkpnh  pnd k20 ¼ ½ðunh  und ; Ph ðunh  Pd unh ÞÞ þ ðunh  und ; Ph ðPd unh  und ÞÞ þ kmðrðunh  Pd unh Þ; rðunh  Pd unh ÞÞ þ kmðrðunh  und Þ; rðPd unh  und ÞÞ þ kðpnh  pnd ; pnh  pnd Þ ¼ ðunh  und ; Ph ðunh  Pd unh ÞÞ þ kmðrðunh  Pd unh Þ; rðunh  Pd unh ÞÞ þ kðdivðPd unh  und Þ;pnh  pnd Þ d n  n n n n n þ ðun1  un1 h d ; Ph ðP uh  ud ÞÞ þ kðph  pd ; ph  pd Þ  d n n  un1 ¼ ðunh  und ; Ph ðunh  Pd unh ÞÞ þ kmðrðunh  Pd unh Þ; rðunh  Pd unh ÞÞ þ ðun1 h d ; Ph ðP uh  ud ÞÞ

þ kðdivðunh  und Þ;pnh  qd pnh Þ þ kðdivðunh  und Þ; qd pnh  pnd Þ þ kðdivðPd unh  unh Þ;pnh  pnd Þ þ kðpnh  pnd ; pnh  pnd Þ d n  n ¼ ðunh  und ; Ph ðunh  Pd unh ÞÞ þ kmðrðunh  Pd unh Þ; rðunh  Pd unh ÞÞ þ ðun1  un1 h d ; Ph ðP uh  ud ÞÞ

þ kðdivðPd unh  unh Þ;pnh  pnd Þ þ kðdivðunh  und Þ; pnh  qd pnh Þ þ kðpnh  pnd ;pnh  pnd Þ  keðpnh  pnd  qh ðpnh  pnd Þ; qd pnh  pnd Þ

990

Z. Luo / Applied Mathematics and Computation 247 (2014) 976–995

6 Chkjunh  und jk0 krðunh  Pd unh Þk0 þ kC½krðunh  Pd unh Þk20 þ kpnh  qd pnh k20  þ

km k krðunh  und Þk20 þ kpnh  pnd k20 4 4

1 1 n 2 n 2 n n n n n n n n d n n þ kjun1  un1 d jk0 þ kjuh  ud jk0 þ kðph  pd ; ph  pd Þ  keðph  pd  qh ðph  pd Þ; q ph  pd Þ 2 h 2 6 kC½krðunh  Pd unh Þk20 þ kpnh  qd pnh k20  þ kðpnh  pnd ; pnh  pnd Þ  keðpnh  pnd  qh ðpnh  pnd Þ; qd pnh  pnd Þ þ

km k 1 1 n 2 n 2 krðunh  und Þk20 þ kpnh  pnd k20 þ kjun1  un1 d jk0 þ kjuh  ud jk0 : 4 2 h 2 2

ð4:19Þ

On the other hand, due to pnh  pnd – qh ðpnh  pnd Þ, by (2.36) and (2.39), we have

ðpnh  pnd  qh ðpnh  pnd Þ; pnh  pnd Þ ¼ ðpnh  pnd  qh ðpnh  pnd Þ; pnh  pnd  qh ðpnh  pnd ÞÞ ¼ ðpnh  pnd  qh ðpnh  pnd Þ; pnh  pnd þ qh ðpnh  pnd ÞÞ ¼ kpnh  pnd k20  kqh ðpnh  pnd Þk20 > 0: ð4:20Þ Therefore, there exists real number

kpnh

j



pnd k20

n h ðph

¼ kq



j 2 ð0; 1Þ such that

pnd Þk20 :

Thus, by taking proper positive real number that

ð4:21Þ

e such that ð1  jÞe P 3=4 and using (2.37), (4.20), and (4.21), we have

ðpnh  pnd ; pnh  pnd Þ  eðpnh  pnd  qh ðpnh  pnd Þ; qd pnh  pnd Þ ¼ ðpnh  pnd ; pnh  pnd Þ  eðpnh  pnd  qh ðpnh  pnd Þ; pnh  pnd Þ  eðpnh  pnd  qh ðpnh  pnd Þ; qd pnh  pnh Þ ¼ ðpnh  pnd ; pnh  pnd Þ  eðpnh  pnd  qh ðpnh  pnd Þ; qd pnh  pnh Þ  eðpnh  pnd  qh ðpnh  pnd Þ; pnh  pnd  qh ðpnh  pnd ÞÞ ¼ kpnh  pnd k20  ekpnh  pnd k20 þ ekqh ðpnh  pnd Þk20  eðpnh  pnd  qh ðpnh  pnd Þ; qd pnh  pnh Þ 1 6 kpnh  pnd k20 þ Ckqd pnh  pnh k20 : 4

ð4:22Þ

If k ¼ OðhÞ, combining (4.22) with (4.19) and using Cauchy inequality yield that 2 d n 2 n n d n 2 kjunh  und jk20 þ kmkrðunh  und Þk20 þ kkpnh  pnd k20 6 kjun1  un1 h d jk0 þ Ck½krðuh  P uh Þk0 þ kph  q ph k0 :

ð4:23Þ

Summing (4.23) from L þ 1 to n ðL þ 1 6 n 6 NÞ and using Lemma 5 yield that

kunh  und k20 þ k

n X

½krðuih  uid Þk20 þ kpih  pid k20  6 Ck

i¼Lþ1

n X

½krðuih  Pd uih Þk20 þ kpih  qd pih k20  þ kuLh  uLd k20 :

ð4:24Þ

i¼Lþ1

Further, we get that

kunh  und k0 þ k

n X

½krðuih  uid Þk0 þ kpih  pid k0  6 CkuLh  uLd k0 þ Ck

i¼Lþ1

1=2

n X

½krðuih  Pd uih Þk0 þ kpih  qd pih k0 :

ð4:25Þ

i¼Lþ1

P Noting that, if k ¼ OðhÞ, by Corollary 13, (3.13), and (2.67) in proof of Theorem 10, there hold that ni¼Lþ1 krðui  P d ui Þk0 6 P n h i h h i i i i i i i i i i¼Lþ1 ½krðuh  u Þk0 þ krðu  P u Þk0 þ krP ðu  uh Þk0  6 C i¼Lþ1 ðkrðu  P u Þk0 þ krðu  uh Þk0 Þ 6 C½h þ hðn  LÞ 6 Pn Pn Pn i d i i i i h i h i i i i i h i Chðn  LÞ and i¼Lþ1 kph  q ph k0 6 i¼Lþ1 ½kph  p k0 þ kp  q p k0 þ kq ðp  ph Þk0  6 C i¼Lþ1 ðkph  p k0 þ kp  q p k0 Þ 6

Pn

2

C½h þ h ðn  LÞ 6 Chðn  LÞ (L 6 n 6 N). Thus, from (4.10) and (4.25), we obtain (4.11), which completes the proof of Theorem 15. h Combining Theorem 10 and 15 yields the following result. Theorem 16. Under hypotheses of Theorem 15, the error estimates between the solutions for Problem II and the solutions for the reduced-order Problem IX are

kunh  und k0 þ k

n X 

 2 krðuih  uid Þk0 þ kpih  pid k0 6 Cðk þ h Þ þ CkL

i¼1

kunh  und k0 þ k

n X

l X

!1=2 kj

;

n ¼ 1; 2; . . . ; L;

ð4:26Þ

j¼dþ1

2

½krðuih  uid Þk0 þ kpih  pid k0  6 Cðk þ h Þ þ Ck

i¼Lþ1

1=2

hðn  LÞ þ CkL

l X

!1=2 kj

; L þ 1 6 n 6 N:

ð4:27Þ

j¼dþ1

Remark 3. Theorems 15 and 16 have presented the error estimates between the solution of ROEA Problem IX based on POD method and the solution of classical SFVE Problem VI and Problem II, respectively. The error estimates in (4.26) show that, in

Z. Luo / Applied Mathematics and Computation 247 (2014) 976–995

991

P 1=2 l order to get better accuracy, L can not be too large and choosing number d of POD basis should satisfy kL 6 j¼dþ1 kj 2

1=2

minfk; h g. While k hðn  LÞ (L þ 1 6 n 6 N) is errors caused by extrapolation, it is a criterion to guide the POD bases renewed duly with the newest solutions of Problem IX. 4.2. Algorithm implementation for solving ROEA Problem IX In the following, we give the implementation of algorithm for solving Problem IX which consists of seven steps. Step 1. Generate the snapshots ensemble

U i ðx; yÞ ¼ ðuih ; pih Þ;

i ¼ 1; 2; . . . ; L  N;

which may be the solutions for Problem VI, physical system trajectories by drawing samples from experiments with interpolation (or data assimilation), or previous results; Step 2. Form the correlation matrix A ¼ ðAik ÞLL , where Aik ¼ 1L ðU i ; U k ÞU~ ; ðU i ; U k ÞU~ ¼ ðruih ; rukh Þ þðpih ; pkh Þ, and ð; Þ is L2 -inner product; Step 3. Let v ¼ ða1 ; a2 ; . . . ; aL ÞT . Solving the eigenvalue problem Av ¼ kv obtains eigenvalues k1 P k2 P . . . P kl > 0 T ðl ¼ dimfU 1 ; U 2 ; . . . ; U L gÞ and corresponding eigenvectors v k ¼ ðak1 ; ak2 ; . . . ; akL Þ ðk ¼ 1; 2; . . . ; lÞ; Step 4. For the triangulation parameter h, the time step increment k, and error d needed, decide on the amount d of POD P 1=2 2 l 6 d; basis such that k þ h þ kL j¼dþ1 kj Step 5. Generate POD basis wk ðx; yÞ ¼ ðwuk ðx; yÞ; wpk ðx; yÞÞ ðk ¼ 1; 2; . . . ; dÞ:

ðwuk ðx; yÞ; wpk ðx; yÞÞ ¼

! L L 1 X 1 X pffiffiffiffiffiffiffi aki uih ; pffiffiffiffiffiffiffi aki pih ; Lkk i¼1 Lkk i¼1

Step 6. Taking U d ¼ span fwu1 ðx; yÞ; wu2 ðx; yÞ; . . . ; wud ðx; yÞg and M d ¼span fwp1 ðx; yÞ; wp2 ðx; yÞ; . . . ; wpd ðx; yÞg and solving Problem IX which only includes 3d degrees of freedom yield the solutions ðund ; pnd Þðn ¼ 1; 2; . . . ; L; L þ 1; . . . ; NÞ; 1=2

Step 7. If k

hðn  LÞ 6 d, then ðund ; pnd Þðn ¼ 1; 2; . . . ; NÞ are the solutions for Problem IX satisfying accuracy needed. Else,

1=2

namely, If k hðn  LÞ > d, let U i ¼ ðuid ; pid Þði ¼ n  L; n  L  1; . . . ; n  1Þ, repeat Step 1 to Step 6 until the solutions ðund ; pnd Þðn ¼ 1; 2; . . . ; L; L þ 1; . . . ; NÞ are obtained. Remark 4. Step 7 could be that if kun1  und k0 P kund  unþ1 k0 and kpn1  pnd k0 P kpnd  pnþ1 k0 ðn ¼ L; L þ 1; . . . ; N  1Þ, then d d d d  und k0 < ðund ; pnd Þðn ¼ 1; 2; . . . ; NÞ are the solutions for Problem IX satisfying accuracy needed. Else, namely, if kun1 d nþ1 n1 n n i i k or kp  p k < kp  p k ðn ¼ L; L þ 1; . . . ; N  1Þ, let U ¼ ðu ; p Þði ¼ n  L  1; n  L  2; . . . ; nÞ, repeat Step kund  unþ1 i 0 0 d d 0 d d d d d 1 to Step 6 until the solutions ðund ; pnd Þðn ¼ 1; 2; . . . ; L; L þ 1; . . . ; NÞ are obtained. 5. Some numerical examples In this section, we present some numerical examples of the physical formulation of backward facing step flow and Reynolds number Re ¼ 103 to validate the feasibility and efficiency of ROEA Problem IX based on POD method. Let the computational field X of the physical formulation of the backward facing step flow be as follows Fig. 2. All initial and boundary value conditions and f ¼ ðf 1 ; f 2 Þ, except inflow of left boundary with a velocity of u ¼ ðu1 ; u2 Þ ¼ ðyð2  yÞ; 0Þ and outflow of right boundary with velocity of u ¼ ðu1 ; u2 Þ satisfying u2 ¼ 0 and m@u1 =@x ¼ 0, are taken as 0. We first divide  into small squares with side length Mx ¼ My ¼ 103 , and then link diagonal of square to divide each square into two the X triangles on same direction and refine the each triangular element p inffiffiffivicinity of corner point ð0; 0Þ into two triangular elements, which constitutes a triangularization Ih with diameter h ¼ 2  103 and nodes  56  106 . In order to satisfy the condition k ¼ OðhÞ in Theorem 15, we take the time step increment as k ¼ 103 . The dual decomposition Ih is taken as barycenter dual decomposition, namely, the barycenter of the right triangle K 2 Ih is taken as the node of the dual decomposition. We first find a set of numerical solution ðunh ; pnh Þ of classical SFVE Problem VI when n ¼ 1; 2; . . . ; 20, constructing 20 snapshots U i ¼ ðuih ; pih Þði ¼ 1; 2; . . . ; 20Þ. And then, by using Step 3 in Section 4.2 we find 20 eigenvalues which are arranged in a non-decreasing order and 20 eigenvectors corresponding to the 20 eigenvalues and construct l POD bases P 1=2 20 6 3  103 . Thus, we only wj ðj ¼ 1; 2; . . . ; l ¼ dimfU 1 ; U 2 ; . . . ; U 20 g). It is shown by computing that eigenvalues kL j¼7 kj need to take the most main 6 eigenvectors ðwuj ; wpj Þðj ¼ 1; 2; . . . ; 6Þ from l POD bases wj ¼ ðwuj ; wpj Þðj ¼ 1; 2; . . . ; lÞ to expand into subspaces U d  M d and compute a numerical solution at t ¼ 2 and Re ¼ 103 with ROEA according to seven steps of Section 4.2, without renewing POD basis we obtain the numerical solutions ðund ; pnd Þðn ¼ 2000, namely, at t ¼ 2), which are depicted graphically in Figs. 3–5 at the bottom chart, but the solutions obtained with the classical SFVE Problem VI are

992

Z. Luo / Applied Mathematics and Computation 247 (2014) 976–995

Fig. 2. The flow region and boundary conditions of the backward facing step flow.

Fig. 3. When Re ¼ 1000, comparison of the component u1 of the velocity u between classical SFVE formulation (top chart: at the time level t ¼ 2) and reduced-order SFVE formulation (bottom chart: at the time level t ¼ 2).

Fig. 4. When Re ¼ 1000, comparison of the component u2 of the velocity u field between the classical SFVE formulation (top chart: at the time level t ¼ 2) and reduced-order SFVE formulation (bottom chart: at the time level t ¼ 2).

depicted respectively graphically in Figs. 3–5 at the top chart at the time level t ¼ 2. Every two charts in Figs. 3–5 are exhibiting quasi-identical similarity, respectively. Fig. 6 shows the errors between solutions obtained ROEA Problem IX with different number of optimal POD bases and solutions obtained with the classical SFVE Problem VI when t ¼ 2 and Re ¼ 1000. Comparing the classical SFVE Problem

Z. Luo / Applied Mathematics and Computation 247 (2014) 976–995

993

Fig. 5. When Re ¼ 1000, comparison of the velocity p field between the classical SFVE formulation (top chart: at the time level t ¼ 2) and reduced-order SFVE formulation (bottom chart: at the time level t ¼ 2).

Fig. 6. Error (log10 ) for Re ¼ 1000 when POD basic is different and at the time level t ¼ 2.

VI with ROEA Problem IX based on POD method containing 6 optimal bases implementing the numerical simulation computations when t ¼ 2 and Re ¼ 1000, we find that for the classical SFVE Problem VI the required computing time is 8 min, while for ROEA Problem IX with 6 optimal bases the corresponding time is only 6 s, namely, the computing time of the classical SFVE Problem VI is 80 times that of ROEA Problem IX with 6 optimal bases, while the error between their solutions does not exceed 3  103 . Our examples use only L solutions the classical SFVE Problem VI at the first few L steps as snapshots, when we compute actual problems, we may structure the snapshots and POD basis with interpolation or data assimilation by drawing samples from experiments, then solve directly ROEA Problem IX, while it is unnecessary to solve the classical SFVE Problem VI. Thus, the time-consuming calculations and resource demands in the computational process will be greatly saved. It is also shown that finding the approximate solutions for non-stationary Stokes equations with ROEA Problem IX is computationally very effective. And the results for numerical examples are consistent with those obtained for the theoretical case. 6. Conclusions and perspective In this paper, we have employed POD technique and FVE method to establish an ROEA for non-stationary Stokes equations. We first compile ensembles of data from the first few LðL  NÞ transient solutions computing an equation system derived with the classical SFVE formulation for non-stationary Stokes equations, while in actual applications, one may obtain the ensemble of snapshots from physical system trajectories by drawing samples from experiments with interpolation (or data assimilation). Next, we employ the POD method to deal with ensembles of data obtaining the POD basis. And then, the trial function spaces of the classical SFVE formulation are replaced with the subspaces spanning with the most main POD basis to derive ROEA for non-stationary Stokes equations. Finally, we provide the error estimates between the classical SFVE solutions and the ROEA solutions based on POD method and the algorithm implementation for solving ROEA of nonstationary Stokes equations. Comparing with the theoretical errors, the error estimates have been verified to provide quite good results, namely both the theoretical errors and the computing errors coincide within plot accuracy, thus validating both the feasibility and efficiency of our ROEA. Though some Galerkin methods based on POD method for parabolic equations (see [31,32]) and a reduced-order FVE formulation based on POD method for parabolic equations (see [48]) has been presented, they are completely different from

994

Z. Luo / Applied Mathematics and Computation 247 (2014) 976–995

here ROEA for non-stationary Stokes equations. Especially, ROEA for non-stationary Stokes equations is different from their reduced-order FD schemes and FE formulations and is of more advantage than their reduced-order FD schemes and FE formulations. The present ROEA based on SFVE method and POD technique has at least three features. Firstly, the treating with non-stationary Stokes equations here is different from parabolic equations. Since non-stationary Stokes equations do not only include the velocity field of fluid but also do the pressure field, establishing ROEA based on SFVE method and POD technique and estimating the errors of ROEA solutions for non-stationary Stokes equations have more great challenge and more difficult than those for parabolic problems in [31,32,48]. Secondly, the method here especially is different from existing other methods (see [35–47]). Almost all existing reduced-order computational methods based on POD technique are actually to test and verify only the comparison of solutions on the same time span ½0; T and belong to repeating computations. In this paper, we do only use the first fewer given numerical solutions on very short time span ½0; T 0  (T 0  T) as snapshots to construct POD basis and establish ROEA based on SFVE method and POD technique for finding the numerical solutions on total time span ½0; T. Thus, we sufficiently adopt the advantage of POD method, namely, sufficiently utilize the given data (on very short time span ½0; T 0  and T 0  T) to forecast future physic phenomenons (on time span ½T 0 ; T), whose idea has very important guiding role if it is extended to weather forecast and nature evolution forecast. Thirdly, in this paper, we first establish a semi-discrete formulation with respect to time for non-stationary Stokes equations and derive the existence and error estimates of its solutions, and then start directly from semi-discrete formulation with respect to time to establish a classical fully discrete SFVE formulation based on two local Gauss integrals and parameter-free and to do error analysis. Thus, we could circumvent the semi-discrete SFVE formulation with respect to space variable such that our method is simpler and more convenient than existing methods (see, e.g., [17,13]) and methodologically different from those in [17,13]. Therefore, the methodology establishing ROEA based on fully discrete SFVE formulation and POD method is the improvement and innovation for the existing methods (see, e.g., [31–48] or others). Future work in this area will aim to extend the reduced-order FVE formulation, implementing it for a realistic atmosphere quality forecast system and more complicated PDEs. Acknowledgement The author would sincerely like to thank the Editors and Reviewers for their valuable suggestions. Especially, the author would like to thank Professor Navon of Florida State University in USA for providing the help and support in English polishing. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28]

F. Brezzi, J. Douglas Jr, Stabilized mixed method for the Stokes problem, Numer. Math. 53 (1988) 225–235. F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer-Verlag, New York, 1991. J. Douglas Jr, J.P. Wang, An absolutely stabilized finite element method for the stokes problem, Math. Comput. 52 (1989) 495–508. V. Girault, P.A. Raviart, Finite Element Methods for Navier–Stokes Equations: Theory and Algorithms, Springer-Verlag, Berlin Heidelberg, 1986. J.G. Heywood, R. Rannacher, Finite element approximation of the non-stationary Navier–Stokes problem, I. Regularity of solutions and second order estimates for spatial discretization, SIAM J. Numer. Anal. 19 (1982) 275–311. T.J. Hughes, L.P. France, A new finite element formulation for computational fluid dynamics, VII. The Stokes problem with various well posed boundary conditions, symmetric formulations that converge for all velocity pressure space, Comput. Methods Appl. Mech. Eng. 65 (1987) 85–96. Z. Cai, S. McCormick, On the accuracy of the finite volume element method for diffusion equations on composite grid, SIAM J. Numer. Anal. 27 (3) (1990) 636–655. W.P. Jones, K.P. Menziest, Analysis of the cell-centred finite volume method for the diffusion equation, J. Comput. Phys. 165 (2000) 45–68. E. Süli, Convergence of finite volume schemes for Poisson’s equation on nonuniform meshes, SIAM J. Numer. Anal. 28 (5) (1991) 1419–1430. R.E. Bank, D.J. Rose, Some error estimates for the box methods, SIAM J. Numer. Anal. 24 (4) (1987) 777–787. R.H. Li, Z.Y. Chen, W. Wu, Generalized Difference Methods for Differential Equations-Numerical Analysis of Finite Volume Methods, Monographs and Textbooks in Pure and Applied Mathematics, 226, Marcel Dekker Inc., New York, 2000. Y. Li, R.H. Li, Generalized difference methods on arbitrary quadrilateral networks, J. Comput. Math. 17 (1999) 653–672. P. Blanc, R. Eymerd, R. Herbin, A error estimate for finite volume methods for the Stokes equations on equilateral triangular meshes, Numer. Methods PDEs 20 (2004) 907–918. P. Chatzipantelidis, R.D. Lazarrov, V. Thomée, Error estimates for a finite volume element method for parabolic equations in convex in polygonal domains, Numer. Meth. PDEs 20 (2004) 650–674. S.H. Chou, D.Y. Kwak, A covolume method based on rotated bilinears for the generalized Stokes problem, SIAM J. Numer. Anal. 35 (1998) 494–507. W.P. Jones, K.P. Menziest, Analysis of the cell-centred finite volume method for the diffusion equation, J. Comput. Phys. 165 (2000) 45–68. M. Yang, H.L. Song, A postprocessing finite volume method for time-dependent Stokes equations, Appl. Numer. Math. 59 (2009) 1922–1932. X. Ye, On the relation between finite volume and finite element methods applied to the Stokes equations, Numer. Meth. PDEs 17 (2001) 440–453. L.H. Shen, J. Li, Z.X. Chen, Analysis of a stabilized finite volume method for the transient stationary Stokes equations, Int. J. Numer. Anal. Model. 6 (3) (2009) 505–519. P. Holmes, J.L. Lumley, G. Berkooz, Turbulence, Coherent Structures, Dynamical Systems and Symmetry, Cambridge University Press, Cambridge, 1996. K. Fukunaga, Introduction to Statistical Recognition, Academic Press, New York, 1990. I.T. Jolliffe, Principal Component Analysis, Springer-Verlag, Berlin, 2002. N. Aubry, P. Holmes, J.L. Lumley, E. Stone, The dynamics of coherent structures in the wall region of a turbulent boundary layer, J. Fluid Dyn. 192 (1988) 115–173. J.L. Lumley, Coherent structures in turbulence, in: R.E. Meyer (Ed.), Transition and Turbulence, Academic Press, New York, 1981, pp. 215–242. H.V. Ly, H.T. Tran, Proper orthogonal decomposition for flow calculations and optimal control in a horizontal CVD reactor, Quart. Appl. Math. 60 (2002) 631–656. A.J. Majda, I. Timofeyev, E. Vanden-Eijnden, Systematic strategies for stochastic mode reduction in climate, J. Atmos. Sci. 60 (2003) 1705–1723. P. Moin, R.D. Moser, Characteristic-eddy decomposition of turbulence in channel, J. Fluid Mech. 200 (1989) 417–509. M. Rajaee, S.K.F. Karlsson, L. Sirovich, Low dimensional description of free shear flow coherent structures and their dynamical behavior, J. Fluid Mech. 258 (1994) 1401–1402.

Z. Luo / Applied Mathematics and Computation 247 (2014) 976–995

995

[29] R.D. Roslin, M.D. Gunzburger, R.A. Nicolaides, G. Erlebacher, M.Y. Hussaini, A self-contained automated methodology for optimal flow control validated for transition delay, AIAA J. 35 (1997) 816–824. [30] L. Sirovich, Turbulence and the dynamics of coherent structures: part I–III, Quart. Appl. Math. 45 (3) (1987) 561–590. [31] K. Kunisch, S. Volkwein, Galerkin proper orthogonal decomposition methods for parabolic problems, Numer. Math. 90 (2001) 117–148. [32] K. Kunisch, S. Volkwein, Galerkin proper orthogonal decomposition methods for a general equation in fluid dynamics, SIAM J. Numer. Anal. 40 (2002) 492–515. [33] K. Kunisch, S. Volkwein, Control of Burgers’ equation by a reduced order approach using proper orthogonal decomposition, J. Optim. Theory Appl. 102 (1999) 345–371. [34] D. Ahlman, F. Södelund, J. Jackson, A. Kurdila, W. Shyy, Proper orthogonal decomposition for time-dependent lid-driven cavity flows, Numer. Heat Trans. Part B–Fund. 42 (2002) 285–306. [35] Y.H. Cao, J. Zhu, Z.D. Luo, I.M. Navon, Reduced order modeling of the upper tropical pacific ocean model using proper orthogonal decomposition, Comput. Math. Appl. 52 (2006) 1373–1386. [36] Y.H. Cao, J. Zhu, I.M. Navon, Z.D. Luo, A reduced order approach to four-dimensional variational data assimilation using proper orthogonal decomposition, Int. J. Numer. Methods Fluids 53 (2007) 1571–1583. [37] J. Du, J. Zhu, Z.D. Luo, I.M. Navon, An optimizing finite difference scheme based on proper orthogonal decomposition for CVD equations, Int. J. Numer. Meth. Biomed. Eng. 27 (2011) 78–94. [38] Z.D. Luo, J. Chen, I.M. Navon, X.Z. Yang, Mixed finite element formulation and error estimates based on proper orthogonal decomposition for the nonstationary Navier–Stokes equations, SIAM J. Numer. Anal. 47 (1) (2008) 1–19. [39] Z.D. Luo, J. Chen, I.M. Navon, J. Zhu, An optimizing reduced PLSMFE formulation for non-stationary conduction–convection problems, Int. J. Numer. Meth. Fluids 60 (4) (2009) 409–436. [40] Z.D. Luo, J. Chen, P. Sun, X.Z. Yang, Finite element formulation based on proper orthogonal decomposition for parabolic equations, Sci. China Ser. A Math. 52 (3) (2009) 585–596. [41] Z.D. Luo, J. Chen, J. Zhu, R.W. Wang, I.M. Navon, An optimizing reduced order FDS for the tropical Pacific Ocean reduced gravity model, Int. J. Numer. Meth. Fluids 55 (2) (2007) 143–161. [42] Z.D. Luo, R.W. Wang, J. Zhu, Finite difference scheme based on proper orthogonal decomposition for the non-stationary Navier–Stokes equations, Sci. China Ser. A: Math. 50 (8) (2007) 1186–1196. [43] Z.D. Luo, X.Z. Yang, Y.J. Zhou, A reduced finite difference scheme based on singular value decomposition and proper orthogonal decomposition for Burgers equation, J. Comput. Appl. Math. 229 (1) (2009) 97–107. [44] Z.D. Luo, Y.J. Zhou, X.Z. Yang, A reduced finite element formulation based on proper orthogonal decomposition for Burgers equation, Appl. Numer. Math. 59 (8) (2009) 1933–1946. [45] Z.D. Luo, J. Zhu, R.W. Wang, I.M. Navon, Proper orthogonal decomposition approach and error estimation of mixed finite element methods for the tropical Pacific Ocean reduced gravity model, Comput. Meth. Appl. Mech. Engin. 196 (41–44) (2007) 4184–4195. [46] P. Sun, Z.D. Luo, Y.J. Zhou, A reduced finite difference scheme based on POD for the non-stationary conduction–convection equations, Math. Numer. Sin. 31 (3) (2009) (323–234). [47] P. Sun, Z.D. Luo, Y.J. Zhou, Some reduced finite difference schemes based on a proper orthogonal decomposition technique for parabolic equations, Appl. Numer. Math. 60 (2010) 154–164. [48] Z.D. Luo, Z.H. Xie, Y.Q. Shang, J. Chen, A reduced finite volume element formulation based on POD for parabolic equations, J. Comput. Appl. Math. 235 (2011) 2098–2111. [49] R.A. Adams, Sobolev Spaces, Academic Press, New York, 1975. [50] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978. [51] Z.D. Luo, Mixed Finite Element Methods and Applications, Chinese Science Press, Beijing, 2006. [52] M. Zlámal, Finite element methods for nonlinear parabolic equations, RAIRO-Anal. Numer. 11 (1) (1977) 93–107. [53] J. Li, Z.X. Chen, A new stabilized finite volume method for the stationary Stokes equations, Adv. Comput. Math. 30 (2009) 141–152. [54] W. Rudin, Functional and Analysis, second ed., McGraw-Hill Companies, Inc., 1973.