A reduced-order LSMFE formulation based on POD method and implementation of algorithm for parabolic equations

Finite Elements in Analysis and Design 60 (2012) 1–12

Contents lists available at SciVerse ScienceDirect

Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel

A reduced-order LSMFE formulation based on POD method and implementation of algorithm for parabolic equations$ Zhendong Luo a,n, Hong Li b,n, Yueqiang Shang c, Zhichao Fang b a

School of Mathematics and Physics, North China Electric Power University, Beijing 102206, China School of Mathematical Sciences, Inner Mongolia University, Huhhot 010021, China c School of Mathematics and Computer Science, Guizhou Normal University, Guiyang 550001, China b

a r t i c l e i n f o

abstract

Article history: Received 12 September 2011 Received in revised form 10 May 2012 Accepted 12 May 2012 Available online 15 June 2012

Proper orthogonal decomposition (POD) technology has been successfully used in the reduced-order modeling of complex systems. In this paper, we extend the applications of POD method to a usual leastsquares mixed finite element (LSMFE) formulation for two-dimensional parabolic equations with real practical applied background. POD bases here are constructed using the method of snapshots. Karhunen–Loeve projection is used to develop a reduced-order LSMFE formulation obtained by projecting the usual LSMFE formulation onto the most important POD bases. The reduced-order LSMFE formulation has lower dimensions and sufficiently high accuracy, is suitable for finding approximate solutions for two-dimensional parabolic equations, and can circumvent the constraint of the convergence stability what is called Brezzi–Babuˇska condition. We derive the error estimates between the reduced-order LSMFE solutions and the usual LSMFE solutions for parabolic equations and provide the implementation of algorithm for solving the reduced-order LSMFE formulation so as to supply scientific theoretic basis and computing way for practical applications. Two numerical examples are used to validate that the results of numerical computations are consistent with theoretical conclusions. Moreover, it is shown that the reduced-order LSMFE formulation based on POD method is feasible and efficient for finding numerical solutions of parabolic equations. & 2012 Elsevier B.V. All rights reserved.

Keywords: Proper orthogonal decomposition Least-squares method Mixed finite element formulation Error estimate Parabolic equation

1. Introduction Many physical phenomena of the natural environment, engineering equipment, and living organisms, such as the proliferation of gas flow, the infiltration of liquid, the conduction of heat, and the spread of impurities in semiconductor materials, are described with parabolic equations. It is not easy to find their analytic solutions for practical engineering problems which usually include complex computing domains, initial value functions, and source terms which are dependent on real practical background; on the contrary, an efficient approach is to find their numerical solutions. Since the mixed finite element (MFE) method can simultaneously find unknown function (displacement) and gradient (stress), its theory is well established, and it is widely applied in scientific computations,

$ This work is jointly supported by National Science Foundation of China (11061009 and 11061021), Natural Science Foundation of Hebei Province (A2010001663), Science Research Program of Guizhou (GJ[2011]2367), and Science Research Program of Inner Mongolia Advanced Education (NJ10006). n Corresponding authors. Tel./fax: þ86 10 61772167. E-mail addresses: [email protected], [email protected] (Z. Luo), [email protected] (H. Li), [email protected] (Y. Shang), [email protected] (Z. Fang).

0168-874X/$ - see front matter & 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.finel.2012.05.002

so it is regarded as one of the most effective methods finding the numerical solutions for parabolic equations (see [1,2]). However, some MFE formulations for parabolic equations include too many degrees of freedom and the combination of finite element subspaces need to satisfy an important convergence stability condition what is called Brezzi–Babuˇska (BB) inequality (see [2,3]) so as to bring about many difficulties for practical applications. Thus, an important problem is how to circumvent the constraint of the BB inequality, alleviate the computational load, and save time-consuming calculations and resource demands in the computational process in a way that guarantees a sufficiently accurate numerical solution. The least-squares mixed finite element (LSMFE) methods can circumvent the constraint of the convergence stability BB condition, which are widely applied in MFE methods for second order elliptic equations, parabolic equations, Stokes equations, the stationary Navier–Stokes equations, the non-stationary conduction– convection problems, and nonlinear non-stationary convection– diffusion problems (see [4–13]). The proper orthogonal decomposition (POD) is a technique offering adequate approximation for representing fluid flow with fewer degrees of freedom, i.e., with lower dimensional models. It was introduced in the context of analysis of turbulent flow by Holmes and Lumley [14]. In other disciplines, the same procedure

2

Z. Luo et al. / Finite Elements in Analysis and Design 60 (2012) 1–12

goes by the names of Karhunen–Loeve decomposition. The POD method has been widely and successfully applied to numerous fields, including signal analysis and pattern recognition (see [15]), statistics (see [16]), geophysical fluid dynamics and meteorology (also see [16]). 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 of the given data. Combining the POD method with some numerical methods for partial differential equations can provide an efficient means of generating some reduced-order models and alleviating the computational load and memory requirements. Though POD method has been widely applied in computations of statistics and fluid dynamics (see [14–25]), it is mainly used to perform principal component analysis and search the main behavior of a dynamic system. More recently, some Galerkin POD methods for parabolic problems and a general equation in fluid dynamics are presented (see [26,27]), the singular value decomposition approach combined with POD technique is used to treat the Burgers equation (see [28]) and the cavity flow problem (see [29]), and some reduced order finite difference models and finite element (or MFE) formulations and error estimates for the upper tropical pacific ocean model, parabolic problems, Burgers equations, the non-stationary Navier–Stokes equations, the nonstationary conduction–convection problems, and CDV equations based on POD method are presented (see [30–43]). However, to the best of our knowledge, there are no published results addressing the case that a combination of POD method with LSMFE methods is used to deal with parabolic equations or providing the error estimates between the usual LSMFE solutions and the reduced-order LSMFE solutions. Therefore, in this paper, the POD technique is applied to a usual LSMFE formulation for parabolic equations such that it is simplified into a reduced-order LSMFE formulation with fewer dimensions and sufficiently high accuracy. The error estimates between the reduced-order LSMFE solutions and the usual LSMFE solutions and, especially, the implementation of algorithm for solving reduced-order LSMFE formulation are provided so as to supply scientific theoretic basis and computing way for practical applications. Finally, it is shown by a numerical example that the results of numerical computations are consistent with theoretical conclusions and validated that the reduced-order LSMFE formulation based on POD method is feasible and efficient for solving two-dimensional parabolic equations. Though Kunisch and Volkwein [26,27] have combined the POD method with Galerkin method to study parabolic equations, and a finite element formulation based on POD method for parabolic has been presented in [35], these methods are completely different from the reduced-order LSMFE formulations for parabolic equations. The present reduced-order LSMFE formulation has three features. One is that the MFE spaces of usual LSMFE formulation are substituted by the subspaces generated with POD basis such that the usual LSMFE formulation with higher dimensions is changed into the reduced-order LSMFE formulation with lower dimensions. The other is that the unknown function (displacement) and gradient (stress) can be simultaneously found. And the last is that the constraint of the convergence stability BB condition can be circumvented. While the reduced-order finite element formulation based on POD method in [35] can only get the unknown function (displacement). Furthermore, the sample values (known as snapshots in POD method) in [26,27] are taken as from the solutions of the physical system at all time instances, while our reduced-order LSMFE formulation needs only to take the solutions of the physical system at the first few time steps as snapshots so that the computational load finding POD bases could be reduced. It is shown that the present method has improved and innovated the existing methods (see, e.g., [26,27,35]).

The paper is organized as follows. Section 2 is to derive the usual LSMFE formulation for parabolic equations and to generate snapshots from the first fewer transient solutions computed from the equation system derived by the usual LSMFE method. In Section 3, the optimal orthonormal bases are reconstructed from the elements of the snapshots with POD method and a reducedorder LSMFE formulation with lower dimensions based on POD method for parabolic equations is established. In Section 4, the error estimates between the usual LSMFE solutions and the reduced-order LSMFE solutions and the implementation of algorithm for solving reduced-order LSMFE formulation are provided. In Section 5, two numerical examples are presented illustrating that the errors between the reduced-order LSMFE approximate solutions and the usual LSMFE solutions are consistent with previously obtained theoretical results, thus validating the feasibility and efficiency of the reduced-order formulation. Section 6 provides main conclusions and future tentative ideas.

2. Usual LSMFE formulation and generation of snapshot Let O  R2 be a bounded, connected and polygonal domain. Consider the following parabolic equations. Problem I. Find u such that 8 > < ut Du ¼ f ,ðx,y,tÞ A O  ð0,TÞ, uðx,y,tÞ ¼ jðx,y,tÞ,ðx,y,tÞ A @O  ½0,TÞ, > : uðx,y,0Þ ¼ u ðx,yÞ,ðx,yÞ A O, 0

where u represents unknown function, T the total time, f the given source term, u0 ðx,yÞ the given initial function, and jðx,y,tÞ the given boundary value function. For the sake of convenience and without loss of generality, we may as well suppose that jðx,y,tÞ is a zero function in the following theoretical analysis. Let p ¼ ru, then Problem I can be written as follows. Problem II. Find (u,p) such that 8 u þdiv p ¼ f , ðx,y,tÞ A O  ð0,TÞ, > > > t > < p þ ru ¼ 0, ðx,y,tÞ A O  ð0,TÞ, uðx,y,tÞ ¼ 0, ðx,y,tÞ A @O  ½0,TÞ, > > > > : uðx,y,0Þ ¼ u ðx,yÞ, ðx,yÞ A O: 0

ð1Þ

The Sobolev spaces used in this context are standard (see [44]). Put X ¼ H10 ðOÞ and M ¼ fq A L2 ðOÞ2 ; div q A L2 ðOÞg equipped with the norm JqJM ¼ ½JqJ20 þ Jdiv qJ20 1=2 . Then a mixed variational formulation for Problem II is written as follows: Problem III. Find ðu,pÞ A X  M such that, for any t A ð0,TÞ, 8 > < ðut þ div p,v þ k div qÞ þðp þ ru,q þ rvÞ ¼ ðf ,v þ k div qÞ, 8ðv,qÞ A X  M, > : uðx,y,0Þ ¼ u ðx,yÞ, ðx,yÞ A O,

ð2Þ

0

where ð,Þ represents L2-inner product. In order to prove the existence and uniqueness of solution of the variational formulation (III), we need the following Gronwall Lemma (see [2,45]). Lemma 1 (Gronwall Lemma). Let g(t) be a positive and integrable function on ½0,T, c be a nonnegative constant. If cðtÞ A C 0 ð½0,TÞ satisfies Z t 0 r cðtÞ rc þ gðsÞcðsÞ ds, 8t A ½0,T, 0

Z. Luo et al. / Finite Elements in Analysis and Design 60 (2012) 1–12

"

then, cðtÞ also satisfies Z t  gðsÞ ds , 0 r cðtÞ r c  exp

2 Ju0 J20 þ ðk þkÞ

r

3

#

n X

i Jf J20

expðnkÞ,

ð7Þ

i¼1

t A ½0,T:

0

Especially, if c ¼0, then cðtÞ  0.

and if the solution u of Problem II satisfies utt A L1 ðOÞ, the following error estimates hold

Theorem 2. If f A L2 ðOÞ and u0 ðx,yÞ A H1 ðOÞ, then there is a unique solution of Problem III.

Juðt n Þun J0 þ

Proof. Since Problem I has a generalized solution u A H10 ðOÞ \ H2 ðOÞ ðt Z0Þ (see [1,2]), therefore, ðu,pÞ is a solution of Problem III, i.e., Problem III has at least one solution. It is necessary to prove uniqueness, i.e., to do that, if f¼0 and u0 ðx,yÞ ¼ 0, Problem III has only zero solution. If f ¼0, u0 ðx,yÞ ¼ 0, and taking v¼0 and q ¼ p ¼ ru, from Problem III we have that 1d JruJ20 þ Jdiv pJ20 þJpJ20 JruJ20 ¼ 0: 2 dt

ð3Þ

Z

t 0

2JruJ20 ds,

8t A ½0,T:

r Ck,

where C in this context indicates a positive constant which is possibly different at different occurrences, being independent of the spatial and temporal mesh sizes.

ð5Þ

Problem IV. Find ðun ,pn Þ A X  M ð1 rn rNÞ such that 8 n n n n > < ðu þ k div p ,v þ k div qÞ þ kðp þ ru ,q þ rvÞ n n1 ¼ ðkf þ u ,v þ k div qÞ, 8ðv,qÞ A X  M, > : u0 ¼ u ðx,yÞ, ðx,yÞ A O:

2

ð6Þ

l 4 0,a0 þ b0 r c0 ,

ai ,

then an þbn r cn expðnl Þ,

2

Jun J20 þ k

Jdiv pi J20 þ k

i¼1 n X

þk

n X

Jpi J20

i¼1

Jr

ui J20 þ k

i¼1

n X

Jpi þ rui J20

i¼1 n X

i

Jf J20 þ k

n1 X

Jui J20 þ Ju0 J20 :

ð11Þ

i¼0

Using discrete Gronwall Lemma for (11) yields (7). Subtracting Problem IV from Problem III at t ¼ t n , taking v ¼ uðt n Þun and q ¼ pðt n Þpn , and using Taylor series yield that ðuðt n Þun þ k divðpðt n Þpn Þ,uðt n Þun þk divðpðtn Þpn ÞÞ þkðpðt n Þpn þ rðuðt n Þun Þ, pðt n Þpn þ rðuðt n Þun ÞÞ

¼ ðuðt n1 Þun1 ,uðt n Þun þ k divðpðt n Þpn ÞÞ Z tn ðt n tÞðutt ,uðt n Þun þ k divðpðt n Þpn ÞÞ dt: þ

ð12Þ

¨ inequality and Cauchy If Jutt J0,1 is a bounded constant, by Holder inequality we obtain that

þ kJpðt n Þpn þ rðuðt n Þun ÞJ20

Theorem 4. Under the assumptions of Theorem 2, Problem IV has a unique solution ðun ,pn Þ A X  M such that Jrui J20 þk

i¼1

n X

Juðt n Þun þ k divðpðt n Þpn ÞJ20

n Z0:

For Problem IV, we have the following results.

n X

ð10Þ

Summing (10) from 1 to n yields that

r

i¼1

n

tn1

i¼0

n X

n

r k Jf J20 þJun1 J20 þkJf J20 þ kJun1 J20 :

i¼1

Lemma 3 (Discrete Gronwall Lemma). If fan g, fbn g, and fcn g are three positive sequences, and fcn g is monotone that satisfy an þbn r cn þ l

n

þ kJrun J20 þkJpn þ run J20 r Jkf þun1 J20

2

The following discrete Gronwall Lemma is well-known and very useful in the next analysis (see [2,45]).

n1 X

2

Jun J20 þ k Jdiv pn J20 þ kJpn J20

r ðk þkÞ

0

2

ð9Þ

Therefore, by Green’s formula we obtain that

For any positive integer N, let k¼T/N denote time-step, and write tn ¼ nk ð0 rn r NÞ. un and pn denote the approximation of uðt n Þ  uðx,t n Þ and pðt n Þ  pðx,t n Þ, respectively. We use a one-step Euler backward finite difference scheme to discretize the time derivative. Then the semi-discrete formulation for Problem III with respect to time t is written as follows.

þk

ð8Þ

Jun þk div pn J20 þ kJpn þ run J20 1 n r ½Jkf þ un1 J20 þJun þ k div pn J20 : 2

By Gronwall Lemma and (5) we get that JruJ0 ¼ 0. Moreover, u ¼0 since Jr  J0 is equivalent to J  J1 in H10 ðOÞ. From (3) we obtain that Jdiv pJ20 þ JpJ20 ¼ 0, which implies that p ¼ 0. Thus, Problem III has a unique solution, which completes the proof of Theorem 2.

Jun J20 þk

1 r n r N,

ð4Þ

namely, JruJ20 r

N N 1X 1X Jpðt i Þpi J0 þ Jrðuðt i Þui ÞJ0 Ni¼1 Ni¼1

þ

Proof. Since the left side of Problem IV is a coercive bilinear functional on X  M, Problem IV has a set of unique solution ðun ,pn Þ A X  M ð1 rn rNÞ. Taking q ¼ pn and v ¼ un in Problem ¨ IV, by Holder inequality and Cauchy inequality we have that

Therefore, we obtain that 1d JruJ20 r JruJ20 , 2 dt

1=2 X N k Jdivðpðt i Þpi ÞJ0 N i¼1

n X

1 ½Juðt n1 Þun1 J20 2 þ Juðt n Þun þk divðpðt n Þpn ÞJ20 4

3

þ Ck þ Ck þ kJuðt n1 Þun1 J20 :

ð13Þ

By using Green’s formula, from (13) we obtain that Jpi J20

i¼1

Jdiv pi J20 þk

n X i¼1

Jpi þ rui J20

2

Juðt n Þun J20 þ k Jdivðpðt n Þpn ÞJ20 þ kJpðt n Þpn J20 þ kJrðuðt n Þun ÞJ20 þkJpðt n Þpn þ rðuðt n Þun ÞJ20 3

r Juðt n1 Þun1 J20 þ Ck þkJuðt n1 Þun1 J20 :

ð14Þ

4

Z. Luo et al. / Finite Elements in Analysis and Design 60 (2012) 1–12

Noting that uðx,y,t 0 Þ ¼ uðx,y,0Þ ¼ u0 ðx,yÞ ¼ u0 . Summing (14) from 1 to n yields that 2

Juðt n Þun J20 þ k

n X

Theorem 7. Under the assumptions of Theorem 4, Problem V has a unique solution sequence ðunh ,pnh Þ A X h  M h ð1 rn rNÞ such that n X

Junh J20 þ k

Jdivðpðt i Þpi ÞJ20

i¼1

i¼1 n X

þk

n X

Jpðt i Þpi J20 þk

i¼1

2

þk

Jrðuðt i Þui ÞJ20

3

r Cnk þ k

Juðt i Þui J20 :

ð15Þ

Jpih J20

i¼1

n X

Jdiv pih J20 þ k

n X

Jpih þ ruih J20

i¼1

Ju0 J20 þ kJ

r

n X

i¼1

"

i¼1 n1 X

Jruih J20 þ k

r

2 u0 J20 þ ðk þ kÞ

2

Juðt n Þun J20 þ k

n X

1=2

1=2

i¼1 n X

þk

n X

Jpðt i Þpi J20 þk

i¼1

Jrðuðt i Þui ÞJ20

2

r Ck expðknÞ:

ð16Þ

PN 2 Noting that for any ai Z 0 ði ¼ 1; 2, . . . ,NÞ, i ¼ 1 ai Zð1=NÞ PN 2 ð i ¼ 1 ai Þ . From (16) we obtain (8), which completes the proof of Theorem 4. Let fIh g be a uniformly regular family of triangulation of O (see [2,3,46]). MFE spaces are taken as follows:

ð18Þ

Jrðun unh ÞJ0

N 1X Jpi pih J0 Ni¼0 i¼0 ! N N X X mþ1 kþ2 i i rC h Ju Jm þ 1 þh Jp Jk þ 1 ,

þ

i¼1

expðnkÞ

and if the solution sequence ðun ,pn Þ of Problem IV satisfy ðun ,pn Þ A Hm þ 1 ðOÞ  Hk þ 1 ðOÞ2 ð0 r n rNÞ and k ¼ OðhÞ, the following error estimates hold Jun unh J0 þ k

Jdivðpðt i Þpi ÞJ20

# i Jf J20

i¼1

i¼0

Using discrete Gronwall Lemma for (15) and noting that nk r T yield that

n X

k N

N X

Jdivðpi pih ÞJ0 þ

i¼0

0 r n r N:

i¼1

ð19Þ

X h ¼ fvh A X \ C 0 ðOÞ; vh 9K A Pm ðKÞ,8K A Ih g,

Proof. Using the same approaches as proving Theorem 4 and Lemma 5, we could prove that Problem V has a unique group of solution ðunh ,pnh Þ A X h  Mh ð1r n r NÞ satisfying (18). Subtracting Problem V from Problem IV taking v ¼ vh and q ¼ qh , we obtain that u0 u0h ¼ u0 P h u0 and

M h ¼ fqh A M; qh 9K A Pk ðKÞ2 , 8K A Ih g,

ðun unh þ k divðpn pnh Þ,vh þk div qh Þ þ kðpn pnh þ rðun unh Þ,qh þ rvn Þ

where m Z1, k Z 1, and Pm(K) is the space of polynomials of degree r m on K. By using dual principle, we easily obtain the following two Lemmas (see [2,3,46]). Lemma 5. There is an operator P h : X-X h such that ðvP h v,vh Þ þkðrðvP h vÞ, rvh Þ ¼ 0,

8vh A X h ,

8v A X:

JP h vJ20 þ kJrP h vJ20 r JvJ20 þ kJrvJ20 : If v A Hr ðOÞ, then rs

9v9r ,s ¼ 1; 0,1,s r r rm þ 1:

Lemma 6. There is an operator rh : M-M h such that

¨ By Lemmas 5 and 6, Green’s formula (20), Holder inequality and Cauchy inequality, we have that

þ ðun1 un1 ,P h un un þk divðrh pn pn ÞÞ h

ðqrh q,qh Þ þ kðdivðqrh qÞ,div qh Þ ¼ 0,

þ ðun1 un1 ,un unh þ k divðpn pnh ÞÞ h

8qh A M h , 8q A M:

r C½Jun P h un J20 þ kJrðun Ph un ÞJ20 þCkJpn rh pn J20

Moreover,

2

þ k Jdivðpn rh pn ÞJ20 þ k Jun P h un J21  k 1 ÞJ20 þ Jun1 un1 J20 þ Jrðun1 un1 h h 2 2 1 þ Jun unh þ k divðpn pnh ÞJ20 : 2

Jrh qJ20 þ kJdiv rh qJ20 rJqJ20 þ kJdiv qJ20 : If q A Hr ðOÞ, then Jqrh qJs r Ch

rs

9q9r ,s ¼ 0; 1,

ð20Þ

Jun unh þk divðpn pnh ÞJ20 þkJpn pnh þ rðun unh ÞJ20 ¼ ðun unh þ k divðpn pnh Þ,un unh þ k divðpn pnh ÞÞ þ kðpn pnh þ rðun unh Þ,pn pnh þ rðun unh ÞÞ ¼ ðun unh þ k divðpn pnh Þ,un Ph un þ k divðpn rh pn ÞÞ þ ðun unh þ k divðpn pnh Þ,Ph un unh þ k divðqh pn pnh ÞÞ þ kðpn pnh þ rðun unh Þ,pn rh pn þ rðun P h un ÞÞ þ kðpn pnh þ rðun unh Þ, rh pn pnh þ rðP h un unh ÞÞ ¼ ðun unh þ k divðpn pnh Þ,un Ph un þ k divðpn rh pn ÞÞ þ kðpn pnh þ rðun unh Þ,pn rh pn þ rðun P h un ÞÞ

Moreover,

JvP h vJs rCh

,vh þk div qh Þ, ¼ ðun1 un1 h 8ðvh ,qh Þ A X h  M h , 1 rn r N:

s rr r k þ1:

1

ð21Þ

Therefore, we obtain that Let ðunh ,pnh Þ be the fully discrete approximation of ðu,pÞ. Then, the fully discrete LSMFE formulation for Problem II may be written as follows. Problem V. Find ðunh ,pnh Þ A X h  M h ð1 r n rNÞ such that 8 n ðu þk div pnh ,vh þk div qh Þ þkðpnh þ runh ,qh þ rvh Þ > > < h n ¼ ðkf þ un1 ,vh þk div qh Þ,8ðvh ,qh Þ A X h  Mh , h > > : u0 ðx,yÞ ¼ P h u0 ðx,yÞ ¼ P h u0 ðx,yÞ,ðx,yÞ A O: h

Jun unh þk divðpn pnh ÞJ20 þ 2kJpn pnh þ rðun unh ÞJ20 r C½Jun P h un J20 2

þ kJrðun Ph un ÞJ20 þk Jdivðpn rh pn ÞJ20 1

þk

Jun Ph un J21 þkJpn rh pn J20 

þ kJrðun1 un1 ÞJ20 þJun1 un1 J20 : h h ð17Þ

Thus, by using Lemmas 5 and 6, we get that 2

Jun unh J20 þ k Jdivðpn pnh ÞJ20

ð22Þ

Z. Luo et al. / Finite Elements in Analysis and Design 60 (2012) 1–12

where JU i J2XM ¼ Jruih J20 þ Jpih J20 þJdivðpih ÞJ20 . A set of solution fwj gdj¼ 1 of (27) and (28) is known as a POD basis of rank d.

þkJpn pnh J20 þ kJrðun unh ÞJ20 2m þ 2

rCðh

2m

þ kh

2 2k

þCðk h

1 2m þ 4

þk

2k þ 2

þ kh

h

ÞJun J2m þ 1

n 2

ÞJp Jk þ 1

þkJrðun1 un1 ÞJ20 þ Jun1 un1 J20 : h h Noting that Ju0 u0h J0 ¼ Ju0 P h u0 J0 r Ch summing (23) from 1 to n yields that

mþ1

ð23Þ 2

Ju0 Jm þ 1 . If k ¼ Oðh Þ,

rC h

i¼1

n X 2m þ 2

Jui J2m þ 1 þh

n X 2k þ 4

i¼0

We introduce the correlation matrix A ¼ ðAij ÞLL A RLL corresponding to the snapshots fU i gLi ¼ 1 by Aij ¼

1 ðU ,U Þ : L i j XM

Proposition 9. Let l1 Z l2 Z    Z ll 40 denote the positive eigenvalues of A and v1 ,v2 , . . . ,vl the associated orthonormal eigenvectors. Then a POD basis of rank d r l is given by

! Jpi J2k þ 1 ,

i¼1

1 rn rN:

ð24Þ PN

2 Noting that for ai Z 0 ði ¼ 1; 2, . . . ,NÞ, we have i ¼ 1 ai Z ð1= P 2 a Þ , from (24) we obtain (19), which completes the NÞð N i¼1 i proof of Theorem 8.

Remark 1. The inequality (18) of Theorem 8 has shown that the solution of Problem (V) is stable and continuously dependent of source term f ðx,y,tÞ and initial conditions u0 ðx,yÞ. If source term f ðx,y,tÞ, initial conditions u0 ðx,yÞ, triangulation parameter h, the time-step k, finite elements Xh and Mh are given, we can obtain a solution ensemble fðunh ðx,yÞ,pnh ðx,yÞÞgN n ¼ 1 by solving Problem V. And then we choose the first L instantaneous solutions ðuih ðx,yÞ,pih ðx,yÞÞ ð1 r ir LÞ (in general, L 5 N, for example, L¼20, N ¼200, which are useful and of interest for us) from the set of N solutions fðunh ðx,yÞ,pnh ðx,yÞÞgN n ¼ 1 for Problem V, which are referred to as snapshots. Remark 2. When one computes actual problems, one may obtain the ensemble of snapshots from physical system trajectories by drawing samples from experiments and interpolation (or data assimilation). 3. Generation of POD basis and a reduced-order LSMFE formulation based on POD method for parabolic equations

L 1 X

wi ¼ pffiffiffiffi ðvi Þj U j , 1r i rd, li j ¼ 1

L d X 1X JU  ðU , w Þ w J2 ¼ L i ¼ 1 i j ¼ 1 i j XM j XM

ðU i , wj ÞXM wj ,

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

ð26Þ

j¼1

where ðU i , wj ÞXM ¼ ðruih , rcuj Þ þ ðpih , wpj Þ þðdiv pih ,div wpj Þ, cuj and wpj are the orthonormal bases corresponding to u and p, respectively. Definition 1. The method of POD consists in finding an orthonormal basis sequence wi ði ¼ 1; 2, . . . ,lÞ such that for every d ð1r d r lÞ the mean square error between the elements U i ð1r i rLÞ and corresponding d-th partial sum of (26) is minimized on average: min

fwj gdj ¼ 1

L d X 1X JU i  ðU i , wj ÞXM wj J2XM Li¼1 j¼1

ð27Þ

ðP h u,vh Þ þkðrPh u, rvh Þ ¼ ðu,vh Þ þ kðru, rvh Þ, d

ðwi , wj ÞXM ¼ dij ,

1 r ir d,1 r j ri,

ð28Þ

8vh A X h :

ð32Þ

d

For any p A M, define div-projection r : M-M denoted by ðrd p,qd Þ þ kðdiv rd p,div qd Þ ¼ ðp,qd Þ þ kðdiv p,div qd Þ,

8qd A M d ,

ð33Þ

where u A X, p A M. Due to (32) and (33) the operators Ph and rd are well-defined and bounded

Jrd pJ0 þ kJdivðrd pÞJ0 r CJpJM ,

8u A X

8p A M:

ð34Þ

ð35Þ

Lemma 10. For every d ð1 rd rlÞ the projection operators Pd and rd satisfy L l X 1X ½Juih Pd uih J20 þ kJrðuih P d uih ÞJ20  rCk lj , Li¼1 j ¼ dþ1

ð36Þ

L l X 1X ½Jpih rd pih J20 þkJdivðpih rd pih ÞJ20  r lj , Li¼1 j ¼ dþ1

ð37Þ

where ðuih ,pih Þ A V ði ¼ 1; 2, . . . ,LÞ are the solution sequence of Problem V. Proof. For any u A X, since uPh u A X, there is a unique wA H2 ðOÞ \ H10 ðOÞ such that kðrw, rvÞ þðw,vÞ ¼ ðv,uP h uÞ, JwJ2 rcJuPh uJ0 :

such that

ð31Þ

Then there is an extension Ph: X-X h of Pd such that Ph 9X ¼ P d : h X h -X d denoted by (see [47])

and

l X

lj :

j ¼ dþ1

ðP d uh ,vd Þ þkðrP d uh , rvd Þ ¼ ðuh ,vd Þ þ kðruh , rvd Þ:

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

and refer to V as the space generated by the snapshots fU i gLi ¼ 1 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 ðl r LÞ. Then each member of the ensemble can be expressed as

l X

Let X d ¼ spanfcu1 , cu2 , . . . , cud g and M d ¼ spanfwp1 , wp2 , . . . , wpd g. For any uh A X h , define the Ritz-projection P d : X h -X d as follows: 8vd A X d ,

JP h uJ0 þ kJrðPh uÞJ0 r CJuJ1 ,

ð25Þ

ð30Þ

where ðvi Þj ðj ¼ 1; 2, . . . ,LÞ denote the j-th component of the eigenvector vi . Furthermore, the following error formula holds

For ðuih ðx,yÞ,pih ðx,yÞÞ ð1 ri rLÞ in Section 2, let U i ðx,yÞ ¼ ðuih ðx,yÞ, pih ðx,yÞÞT ð1 r ir LÞ and

Ui ¼

ð29Þ

Since the matrix A is positive semi-definite and has rank l, the solution of (27) and (28) can be found, in addition, the following result holds (for example, see [14,25–27,35]).

Jun unh J20 þ kJrðun unh ÞJ20 n n X X 2 þk Jdivðpi pih ÞJ20 þk Jpi pih J20 i¼1

5

8v A X, ð38Þ

Let ph w A X h be a interpolation function of w in Xh. Then, by ¨ Holder inequality and Cauchy inequality, we deduce from (32)

6

Z. Luo et al. / Finite Elements in Analysis and Design 60 (2012) 1–12

Theorem 4 and (18) of Theorem 7, we could obtain that

and theory of interpolation (see [2,3,46]) that JuP

h

uJ20

h

h

¼ kðrw, ruP uÞ þ ðw,uP uÞ

Jund J20 þ k

¼ kðrðwph wÞ, rðuP h uÞÞ þ ðwph w,uPh uÞ 2

2

þk

2

rCkhJuP h uJ0 JrðuP h uÞJ0 þCh JuPh uJ20

"

2 2

rCk h JrðuPh uÞJ20 1 2 JuP h uJ20 þ Ch JuP h uJ20 : 4

2

ð40Þ

¨ Thus, by Holder inequality and Cauchy inequality, we deduce from (32) and (40) that

¼ ðuPh u,uPh uÞ þ kðrðuPh uÞ, rðuPh uÞÞ ¼ ðuPh u,uvd Þ þkðrðuP h uÞ, rðuvd ÞÞ r CJuP h uJ21 Jrðuvd ÞJ20 k JrðuPh uÞJ20 þ kJrðuvd ÞJ20 4

k JrðuP h uÞJ20 þ CkJrðuvd ÞJ20 , 2

8vd A X d  X h :

ð41Þ

Therefore, we obtain that JuP h uJ20 þ kJrðuP h uÞJ20 r CkJrðuvd ÞJ20 ,8vd A X d :

ð42Þ

If u ¼ uih , and Ph is restricted to Ritz-projection from Xh to Xd, i.e., P P h uih ¼ P d uih A X d , taking vd ¼ dj¼ 1 ðuih , cj ÞX cj A X d in (42), we can obtain (36) from (31). ¨ Using Holder inequality, Cauchy inequality, and (33) yields Jpih  d pih J20 þ kJdivðpih  d pih ÞJ20 ¼ ðpih  d pih ,pih  d pih Þ

r

r

þkðdivðpih rd pih Þ,divðpih rd pih ÞÞ ¼ ðpih rd pih ,pih qd Þ þ kðdivðpih rd pih Þ,divðpih qd ÞÞ r

n X

1 i 1 Jp rd pih J20 þ Jpih qd J20 2 h 2 k k þ Jdivðpih rd pih ÞJ20 þ Jdivðpih qd ÞJ20 , 2 2

Jpid þ ruid J20

i¼1 n X

# i

Jf J20 expðnkÞ,

8qd A M d ,

consequently, Jpih rd pih J20 þ kJdivðpih rd pih ÞJ20 r Jpih qd J20 þkJdivðpih qd ÞJ20 ,8qd A M d : ð43Þ Pd i Taking qd ¼ j ¼ 1 ðph , wpj Þ0 wpj in (43) and using (31) yield (37), which completes the proof of Lemma 10. Thus, substituting Xd and Md for Xh and Mh, respectively, we can obtain a reduced-order formulation for Problem V as follows. Problem VI. Find ðund ,pnd Þ A X d  M d ð1 r n r NÞ s such that 8 n ðu þk div pnd ,vd þk div qd Þ þ kðpnd þ rund ,qd þ rvd Þ > > < d n ¼ ðkf þ un1 ,vd þ k div qd Þ, 8ðvd ,qd Þ A X d  Md , d > > : u0 ðx,yÞ ¼ P h u0 ðx,yÞ, ðx,yÞ A O:

ð45Þ

which has shown that the solution of Problem VI is stable and continuously dependent of source term f ðx,y,tÞ and initial conditions u0 ðx,yÞ. n n n n Let pnd ¼ ðpn1d ,pn2d Þ ¼ ðb11 cp1 1 þ b12 cp1 2 þ    þ b1d cp1 d , b21 cp2 1 þ n n n n n b22 cp2 2 þ    þ b2d cp2 d Þ and ud ¼ a1 cu1 þ a2 cu2 þ    þ and cud . Then, by Green’s formula, Problem VI can be rewritten as follows. n

JuP h uJ20 þ kJrðuP h uÞJ20

r

Jdiv pid J20 þ k

i¼1

ð39Þ

JuP h uJ0 rCk JrðuP h uÞJ0 :

r

i¼1

Jpid J20

i¼1

2

2

r

n X

n X

r Ju0 J20 þ kJru0 J20 þ ðk þ kÞ

If h is sufficiently small such that Ch r 1=4 and k ¼ OðhÞ, from (39) we get that

þ

Jruid J20 þ k

i¼1

rCkhJwJ2 JrðuP h uÞJ0 þCh JwJ2 JuPh uJ0

þ

n X

n

n

n

n

Problem VII. Find ðan1 , an2 , . . . , and , b11 , b12 , . . . , b1d , b21 , b22 , . . . , bn2d ÞT A R3d ðn ¼ 1; 2, . . . ,NÞ such that 8 d d X X > n > n > a B ¼ kðf , c Þ þ an1 ðcui , cuj Þ, i ¼ 1; 2, . . . ,d; > ij ui j j > > > j¼1 j¼1 > > > > d > > >X n 2 n n > ½b1j C ij þ b2j Dij  ¼ k ðf ,@cp1 i =@xÞ > > > > j ¼ 1 > > > > > d X < þk an1 ðcuj ,@cp1 i =@xÞ,i ¼ 1; 2, . . . ,d; j > > j ¼ 1 > > > > d > >X n ~ n ~ 2 n > > ½b D > ij þ b2j C ij  ¼ k ðf ,@cp2 i =@yÞ > > j ¼ 1 1j > > > > > d > X > > > þk an1 ðcuj ,@cp2 i =@yÞ,i ¼ 1; 2, . . . ,d, > j > : j¼1

Pd 0 where a0j ðj ¼ 1; 2, . . . ,dÞ satisfy ðu0 , cui Þ j ¼ 1 aj ðcui , cuj Þ ¼ 2 ði ¼ 1; 2, . . . ,dÞ, Bij ¼ ðcui , cuj Þ þ kðrcui , rcuj Þ, C ij ¼ k ð@cp1 j =@x, 2 ~ ij ¼ k2 ð@c = @cp1 i =@xÞ þkðcp1 j , cp1 i Þ, Dij ¼ k ð@cp2 j =@y,@cp1 i =@xÞ, D p1 j 2 @x,@cp2 i =@yÞ, C~ ij ¼ k ð@cp2 j =@y,@cp2 i =@yÞ þkðcp2 j , cp2 i Þ, and ð,Þ is L2-inner product. Remark 3. If Ih is a uniformly regular triangulation and Xh and Mh, respectively, are the spaces of piecewise linear polynomials and polynomial vectors, the number of total degrees of freedom for Problem V, i.e., the number of unknown quantities is 3N h (where Nh is the number of vertices of triangles in Ih ), while the number of total degrees of freedom for Problem VI, i.e., Problem VII is 3d ðd 5 l r L5 NÞ. For scientific engineering problem, the number Nh of vertices of triangles in Ih is more than ten of thousands, even more than a hundred million, while d is only the number of the first few maximal eigenvalues so that it is very small (for example, the first example in Section 5, d ¼6, while Nh ¼ 2  103  2  103 ¼ 4  106 ). Therefore, Problem VI is a reduced-order LSMFE formulation based on POD method for Problem V. 4. Error estimates of solution for problem VI 4.1. Error estimates of solutions for Problem VI

ð44Þ

In this subsection, we adopt the notation from Section 3 and recur to usual LSMFE method to derived the error estimates for Problem VI. We have the following main result.

Since ðund þ k div pnd ,vd þ k div qd Þ þkðpnd þ rund ,qd þ rvd Þ in X d  M d is positive definite, Problem VI has a unique solution sequence ðund ,pnd Þ A X d  Md . Using the same approaches as proving (7) of

Theorem 11. Under hypotheses of Theorem10, if ðunh ,pnh Þ A X h  Mh , ðn ¼ 1; 2, . . . ,NÞ are the solutions of Problem V, ðund ,pnd Þ A X d  Md ðn ¼ 1; 2, . . . ,NÞ are the solutions of Problem VI, k ¼ OðhÞ, and

d

Z. Luo et al. / Finite Elements in Analysis and Design 60 (2012) 1–12

L2 ¼ OðNÞ, then for n ¼ 1; 2, . . . ,L, the following error estimates hold 1=2

Junh und J0 þk Jrðunh und ÞJ0 n n X 1 X þk2 Jpjh pjd J0 þ k Jdivðpjh pjd ÞJ0 j¼1

0 1=2

rC @k

11=2

lj A

Summing (49) from 1; 2, . . . ,n, noting that u0h ¼ u0d , if k ¼ OðhÞ, 1 k ¼ OðNÞ ¼ OðL2 Þ, from Lemma 10 we obtain that Junh und J20 þkJrðunh und ÞJ20 n n X X 2 þk Jpjh pjd J20 þ k Jdivðpjh pjd ÞJ20

j¼1 l X

7

j¼1

:

ð46Þ

j ¼ dþ1

rC

n X

j¼1

ðJujh P d ujh J20 þ kJrðujh Pd ujh ÞJ20

j¼1 2

Proof. Since X d  X h and M d  M h , subtracting Problem VI from Problem V taking vh ¼ vd A X d and qh ¼ qd A M d yields that u0h u0d ¼ 0 and ðunh und þ k divðpnh pnd Þ,vd þ k div qd Þ un1 ,vd þ k div qd Þ, ¼ ðun1 h d ð47Þ

¨ By (32), (33), (47), Green’s formula, Holder inequality, Cauchy inequality, and (40), we have Junh und J20 þ kJrðunh und ÞJ20

l X

lj :

ð50Þ

j ¼ dþ1

Combining Theorems 4, 7, and 11 yields the following result. Pn i Theorem 12. Under hypotheses of Theorem 11, if i ¼ 1 Ju Jm þ 1 Pn i and i ¼ 1 Jp Jk þ 1 are all bounded, the error estimates between the solutions for Problem II and the solutions for the reduced-order Problem VI are, for n ¼ 1; 2, . . . ,L,

2

n X j¼1

Jpðt j Þpjd J0 þ k

n X

Jdivðpðt j Þpjd ÞJ0

j¼1

0 11=2 3 l X 6 m k @ 1=2 r C 4k þh þ h þ k lj A 7 5:

þkðpnh pnd þ rðunh und Þ, pnh pnd þ rðunh und ÞÞ ¼ ðunh und þk divðpnh pnd Þ, unh Pd unh þ k divðpnh rd pnh ÞÞ

ð51Þ

j ¼ dþ1

þkðpnh pnd þ rðunh und Þ, pnh rd pnh þ rðunh P d unh ÞÞ þðunh und þk divðpnh pnd Þ,P d unh und þ k divðrd pnh pnd ÞÞ þ kðpnh pnd þ rðunh und Þ, qd pnh pnd þ rðP d unh und ÞÞ ¼ Junh P d unh J20 þ kJrðunh Pd unh ÞJ20 2

þkJpnh rd pnh J20 þk Jdivðpnh rd pnh ÞJ20 þðun1 un1 , P d unh und þk divðrd pnh pnd ÞÞ h d ¼ Junh P d unh J20 þ kJrðunh Pd unh ÞJ20 2

þkJpnh rd pnh J20 þk Jdivðpnh rd pnh ÞJ20 þðun1 un1 , P d unh unh þk divðrd pnh pnh ÞÞ h d un1 , unh und þ k divðpnh pnd ÞÞ þðun1 h d rC½Junh P d unh J20 þkJrðunh P d unh ÞJ20 2

þkJpnh rd pnh J20 þk Jdivðpnh rd pnh ÞJ20  þCJrðun1 un1 ÞJ0 Junh P d unh J1 h d 1 n1 n1 2 Ju ud J0 2 h 1 þ Junh und þ k divðpnh pnd ÞJ20 2 þ

rC½Junh P d unh J20 þkJrðunh P d unh ÞJ20 2

þkJpnh rd pnh J20 þk Jdivðpnh rd pnh ÞJ20  þ

ð48Þ

Thus, from (48) we have Junh und J20 þ kJrðunh und ÞJ20 þ kJpnh pnd J20 2

þk Jdivðpnh pnd ÞJ20 þ kJpnh pnd þ rðunh und ÞJ20 rCðJunh P d unh J20 þkJrðunh P d unh ÞJ20 2

þkJpnh rd pnh J20 þk Jdivðpnh rd pnh ÞJ20 Þ þkJr

lj r k1=2

j ¼ dþ1

1

¼ ðunh und þk divðpnh pnd Þ, unh und þ k divðpnh pnd ÞÞ

ðun1 un1 ÞJ20 þ Jun1 un1 J20 : h d h d

r kL

Juðt i Þuid J0 þ k2

2

þkJpnh pnd J20 þ k Jdivðpnh pnd ÞJ20

k 1 Jrðun1 un1 ÞJ20 þ Jun1 un1 J20 h d d 2 2 h 1 þ Junh und þ k divðpnh pnd ÞJ20 : 2

l X

By (50) we get the second inequality of (46), which completes Theorem 11.

þkðpnh pnd þ rðunh und Þ,qd þ rvd Þ 8ðvd ,qd Þ A X d  M d ,1 rn r N:

þ kJpjh rd pjh J20 þ k Jdivðpjh rd pjh ÞJ20 Þ

ð49Þ

Remark 4. The condition L2 ¼ OðNÞ in Theorem 11 shows the relation between the number L of snapshots and the number N of all time instances. Therefore, it is unnecessary to take total transient solutions at all time instances tn as snapshots (see [26,27]). Theorems 11 and 12 have, respectively, presented the error estimates between the solution of the reduced-order LSMFE formulation Problem VI and the solution of usual LSMFE formulation Problem V and Problem II, which can guide one to choose the number d of POD basis and the number L of snapshots such that 1=2 Pl minfm, kg minfm, kg 1=2 1=2 ðk ¼ Oðh Þ ¼ OðkÞ and h þ k þðk j ¼ d þ 1 lj Þ Pl 1=2 satisfies the accuracy needed. Our method here j ¼ d þ 1 lj Þ employs the first L solutions ðunh ,pnh Þ ðn ¼ 1, 2, . . . ,LÞ for Problem V as the snapshots and to construct the POD bases. However, when one computes actual problems, one may obtain the ensemble of snapshots from physical system trajectories by drawing samples from experiments and interpolation (or data assimilation). Therefore, the snapshots ðunh ,pnh Þ ðn ¼ 1; 2, . . . ,LÞ could be replaced with the interpolation functions of experimental and previous results, thus rendering it unnecessary to solve Problem V, and requiring only to solve directly Problem VI such that Theorem 11 is satisfied. Moreover, since the development and change of a large number of future nature phenomena are closely related to previous given results, for example, weather change, biology anagenesis, by using existing results (which may be from physical system trajectories by drawing samples from experiments) as snapshots to construct POD bases and substituting the corresponding MFE spaces by the subspace generated with POD basis, the high order or infinity dimensional dynamical system, i.e., PDEs which can describe the laws of change of natural phenomena is reduced a lower order dynamical system with fewer dimensions. Thus, the future laws of change of natural phenomena can be quickly simulated and forecasted. And then, time instances are continuously extrapolated forward and POD basis is ceaselessly renewed, the rules of future development and change of natural phenomena would be very well simulated, which is a result of major importance for real-life applications.

8

Z. Luo et al. / Finite Elements in Analysis and Design 60 (2012) 1–12

4.2. Implementation of algorithm for solving Problem VI In the following, we give the implementation of algorithm for solving Problem VI which consists of seven steps. Step 1: Generate the snapshots ensemble U i ðx,yÞ ¼ ðuih ðx,yÞ,pih ðx,yÞÞT ,1 ri r L, which may be the solutions for Problem V, physical system trajectories by drawing samples from experiments with interpolation (or data assimilation), or previous given results. Step 2: Generate the correlation matrix A ¼ ðAij ÞLL , where Aij ¼ ðU i ,U j ÞXM =L ¼ ½ðruih , rujh Þ þ ðpih ,pjh Þ þ ðdiv pih ,div pjh Þ=L and ð,Þ is L2-inner product. Step 3: Solving the eigenvalue problem Av ¼ lv,

v ¼ ða1 ,a2 , . . . ,aL ÞT

ð52Þ

obtains eigenvalues l1 Z l2 Z    Z ll 40ðl ¼ dimfU 1 ,U 2 , . . . ,U L gÞ and corresponding eigenvectors vj ¼ ðaj1 ,aj2 , . . . ,ajL Þ ðj ¼ 1; 2, . . . ,lÞ. Step 4: For the triangulation parameter h, the time step increment k, and the error d needed, decide on the amount d of POD minfm, kg 1=2 Pl 1=2 basis, m, and k such that k þ h þðk r d. j ¼ d þ 1 lj Þ P Step 5: q Generate POD basis wj ¼ ðcuj , cp1 j , cp2 j Þ ¼ Li ¼ 1 aji ffiffiffiffiffiffiffi ðuih ,pi1h ,pi2h Þ= Llj ðj ¼ 1; 2, . . . ,dÞ. 2 Step 6: Let Bij ¼ ðcui , cuj Þ þkðrcui , rcuj Þ, C ij ¼ k ð@cp1 j =@x, 2 2 ~ @cp1 i =@xÞ þ kðcp1 j , cp1 i Þ, Dij ¼ k ð@cp2 j =@y,@cp1 i =@xÞ, D ij ¼ k ð@cp1 j = 2 @x,@cp2 i =@yÞ and, C~ ij ¼ k ð@cp2 j =@y,@cp2 i =@yÞ þ kðcp2 j , cp2 i Þ. Solving the following system of equations which only includes 3d degrees of freedom: 8 d X > > > a0j ðcui , cuj Þ ¼ ðu0 , cui Þ, i ¼ 1; 2, . . . ,d; > > > > j¼1 > > > > d d > X > >X n > aj Bij ¼ kðf n , cui Þ þ an1 ðcui , cuj Þ, i ¼ 1; 2, . . . ,d; > j > > > j¼1 j¼1 > > > > > d > >X n n 2 n > > ½b1j C ij þ b2j Dij  ¼ k ðf ,@cp1 i =@xÞ > > X > > > þk an1 ðcuj ,@cp1 i =@xÞ,i ¼ 1; 2, . . . ,d; > j > > > j ¼ 1 > > > > d > X > n ~ n ~ 2 n > > ½b1j D > ij þ b2j C ij  ¼ k ðf ,@cp2 i =@yÞ > > > >j¼1 > > > d > X > > > þk an1 ðcuj ,@cp2 i =@yÞ,i ¼ 1; 2, . . . ,d > j > : j¼1

obtains n n n n n n ðan1 , an2 , . . . , and , b11 , b12 , . . . , b1d , b21 , b22 , . . . , b2d ÞT A R3d ðn ¼ 1; 2, . . . ,NÞ n n n n and further yields the solutions pd ¼ ðp1d ,pn2d Þ ¼ ðb11 cp1 1 þ b12 cp1 2 n n n n þ    þ b1d cp1 d , b21 cp2 1 þ b22 cp2 2 þ    þ b2d cp2 d Þ and und ¼ an1 cu1 þ an2 cu2 þ    þ and cud ðn ¼ 1; 2, . . . ,NÞ of Problem VI;

Step 7: If Jðun1 ,pn1 Þðund ,pnd ÞJXM Z Jðund ,pnd Þðund þ 1 , pnd þ 1 Þ d d JXM ðn ¼ L,L þ1, . . . ,N1Þ, then ðund ,pnd Þ ðn ¼ 1; 2, . . . ,NÞ are the solutions for Problem VI, whose errors are not more than minfm, kg 1=2 Pl 1=2 kþ h þðk . Else, i.e., if Jðun1 ,pn1 Þ j ¼ d þ 1 lj Þ d d ðund ,pnd ÞJXM o Jðund ,pnd Þðund þ 1 ,pnd þ 1 ÞJXM ðn ¼ L,L þ1, . . . ,N1Þ, let U i ¼ ðuid ,pid Þ ði ¼ nL1,nL2, . . . ,nÞ and repeat Step 1–Step 6.

5. Two numerical experiments In this section, we present two numerical examples showing the advantage of the reduced-order POD LSMFE formulation, i.e., Problem VI for the two-dimensional parabolic equations.

5.1. The case of regular domain and discretization In this subsection, we consider the simple two-dimensional parabolic equations, whose computational domain is regular and consists of a square with dimensions 2 cm  2 cm, namely, O ¼ ½0; 2  ½0; 2, the source term is taken as f ðx,y,tÞ ¼ 2p2 expðp2 tÞ sin px sin py, the initial function is taken as u0 ðx,yÞ ¼ sin px sin py which satisfies uðx,y,tÞ ¼ 0 ððx,yÞ A @O  ½0,TÞ, and MFE spaces Xh and Mh are, respectively, taken as piecewise linear polynomials and piecewise linear polynomial vectors. Though it is a very simple example, the ideas and approaches here can be easily applied to numerical computations for other two-dimensional parabolic equations models with real practical background. We first divide the field O into 2000  2000 small squares with side length Dx ¼ Dy ¼ 103 , and then link the diagonal of the square to divide each square into two triangles in the pffiffiffisame direction which consists of triangularization Ih . Thus h ¼ 2  103 . In order to make k ¼ OðhÞ satisfied, we take time-step as k ¼ 103 . We find the numerical solutions unh and pnh  ðpn1h ,pn2h Þ with usual LSMFE method, i.e., Problem V when t¼ 1, depicted graphically on right hand side in Figs. 1–3, respectively. We output a set of 1000 values at time t ¼ 0:001,0:002, . . . ,1 when solving the usual LSMFE formulation, i.e., Problem V. And then, we choose the first 20 values from 1000 values to constitute a set of snapshots. It is shown by computing that when d ¼6 and k ¼ 103 , the error 1=2 P20 1=2 ðk r 2  103 in Theorem 11, which shows that one j ¼ 7 lj Þ only needs to choose six POD bases. When we solve the reducedorder Problem VI based on POD method with six optimal POD bases, according seven steps of implementation of algorithm in Subsection 4.2, we find that the algorithm of reduced-order LSMFE formulation at t ¼1 are still convergent, without renewing POD basis. The reduced-order LSMFE solutions obtained with the reduced-order Problem VI are depicted graphically on left hand side in Figs. 1–3, respectively. Every two images in Figs. 1–3 are exhibiting quasi-identical similarity, respectively, but the reduced-order solution is better than usual LSMFE solution (it is

Fig. 1. Figure on left hand side is a figure of solution und of reduced-order LSMFE formulation when d ¼ 6, t¼ 1; figure on right hand side is a figure of solution unh of usual LSMFE formulation when t ¼1.

Z. Luo et al. / Finite Elements in Analysis and Design 60 (2012) 1–12

9

Fig. 2. Figure on left hand side is a figure of pn1d solution of the reduced-order LSMFE formulation when d ¼ 6, t¼ 1; figure on right hand side is a figure of pn1h solution of usual LSMFE formulation when t¼ 1.

Fig. 3. Figure on left hand side is a figure of pn2d solution of the reduced-order LSMFE formulation when d ¼ 6, t¼ 1; figure on right hand side is a figure of pn2h solution of usual LSMFE formulation when t¼ 1.

only 8 s, namely, the required computing time to solve the usual LSMFE formulation Problem V is as 120 times as that to do the reduced-order LSMFE formulation Problem VI with six POD bases, while the errors between their respective solutions do not exceed 2  103 . 5.2. The case of irregular domain and discretization In this subsection, we consider the two-dimensional parabolic equations, whose computational domain is irregular and consists of a set O ¼ ð½0; 2  ½0; 2Þ [ ð½0:65,1:3  ½2,2:03Þ cm2 , the source term is taken as f ðx,y,tÞ ¼ 0, the initial and boundary valve functions are taken as follows, for 0 r t rT,

jðx,y,tÞ ¼ u0 ðx,yÞ Fig. 4. When t¼ 1, the errors ðlog 10Þ between solutions of Problem VI with different number of POD bases for a group of 20 snapshots and the usual LSMFE formulation Problem V with piecewise 1st degree polynomials and piecewise 1st degree polynomial vectors.

due to use more initial numerical data, i.e., six POD bases which is our opinion). Fig. 4 shows the errors between 20 solutions ðund ,pnd Þ  ðund ,pn1d ,pn2d Þ of the reduced-order Problem (VI) with 20 different numbers of POD bases and the solutions ðunh ,pnh Þ  ðunh ,pn1h ,pn2h Þ of usual LSMFE formulation, namely, Problem V at t ¼1, respectively. Comparing the usual LSMFE formulation Problem (V) with the reduced-order LSMFE formulation Problem (VI) containing six POD bases implementing the numerical simulation computations when total time T ¼1, we find that for usual LSMFE formulation Problem (V) with piecewise linear polynomials for unh and piecewise linear polynomial vector for pnh  ðpn1h ,pn2h Þ, which includes 3  2  103  2  103 ¼ 12  106 degrees of freedom, the required computing time is approximately 16 min, while for the reducedorder LSMFE formulation Problem VI with six POD bases, which includes only 18 degrees of freedom, the corresponding time is

8 > < 1:5 ¼ 1:0 > : 0:0

if ðx,yÞ A ½0:65,1:3  ½2,2:03, if ðx,yÞ A ð½0,0:5 [ ð½1:5,2Þ  ½2; 2, others:

The MFE spaces Xh and Mh are still taken as piecewise linear polynomials and piecewise linear polynomial vectors, respectively. The triangularization Ih for computational domain O adopts local refining meshes such that the scale of meshes on ½0:65,1:3  ½2,2:03 and nearby ðx,2Þ ð0 r xpr2Þ are one-third of meshes ffiffiffi nearby ðx,0Þ ð0 r x r2Þ. But h ¼ 2  102 . The time-step is still k ¼ 103 . We also find the numerical solutions unh and pnh  ðpn1h ,pn2h Þ with usual LSMFE method, i.e., Problem (V) when t ¼1, depicted graphically on right hand side in Figs. 5–7, respectively. We output a set of 1000 values at time t ¼ 0:001,0:002, . . . ,1 when solving the usual LSMFE formulation, i.e., Problem V. And then, we still choose the first 20 values from 1000 values to constitute a set of snapshots. It is shown by computing that when d ¼6 and 1=2 P20 1=2 k ¼ 103 , the error ðk r4  104 in Theorem 11, j ¼ 7 lj Þ which shows that one only needs to choose six POD bases. When we solve the reduced-order Problem VI based on POD method

10

Z. Luo et al. / Finite Elements in Analysis and Design 60 (2012) 1–12

Fig. 5. Figure on left hand side is a figure of solution ud formulation when t¼ 1.

n

of reduced-order LSMFE formulation when d ¼6, t¼ 1; figure at bottom is a figure of solution unh of usual LSMFE

Fig. 6. Figure on left hand side is a figure of pn1d solution of the reduced-order LSMFE formulation when d ¼ 6, t¼ 1; figure on right hand side is a figure of pn1h solution of usual LSMFE formulation when t ¼ 1.

Fig. 7. Figure on left hand side is a figure of pn2d solution of the reduced-order LSMFE formulation when d ¼ 6, t¼ 1; figure on right hand side is a figure of pn2h solution of usual LSMFE formulation when t ¼ 1.

with six optimal POD bases, according seven steps of implementation of algorithm in Subsection 4.2, we find that the algorithm of reduced-order LSMFE formulation at t ¼1 are still convergent, without renewing POD basis too. The reduced-order LSMFE solutions obtained with the reduced-order Problem VI are depicted graphically on left hand side in Figs. 5–7, respectively. Every two images in Figs. 5–7 at this case are also exhibiting quasi-identical similarity, respectively, but the reduced-order solution is far better than usual LSMFE solution (it is due to use more initial numerical data, namely, six POD bases which is our viewpoint). Fig. 8 shows the absolute errors between 20 solutions ðund ,pnd Þ  ðund ,pn1d ,pn2d Þ of the reduced-order Problem VI with 20 different numbers of POD bases and the solutions ðunh ,pnh Þ  ðunh ,pn1h ,pn2h Þ of usual LSMFE formulation, namely, Problem V at t ¼1, respectively. And again, comparing the usual LSMFE formulation Problem V with the reduced-order LSMFE formulation Problem VI containing six POD bases implementing the numerical simulation computations when total time T¼ 1, we find that for

usual LSMFE formulation Problem V with piecewise linear polynomials for unh and piecewise linear polynomial vector for pnh  ðpn1h ,pn2h Þ, the required computing time is approximately 10 min, while for the reduced-order LSMFE formulation Problem VI with six POD bases, the corresponding time is only 8 seconds too, namely, the required computing time to solve the usual LSMFE formulation Problem V is as 100 times as that to do the reduced-order LSMFE formulation Problem VI with six POD bases, while the errors between their respective solutions do not exceed 4  104 . Due to adopting local refining meshes, numerical simulation errors are improved greatly. Though our two examples are in sense recomputing what we have already computed the first L ðL5 NÞ steps by usual LSMFE formulation, when we compute actual problems, we may construct the snapshots and POD basis with interpolation or data assimilation by drawing samples from experiments, then solve directly the reduced-order LSMFE formulation, while it is unnecessary to solve usual LSMFE formulation, thus, the time-consuming calculations and resource demands in the computational

Z. Luo et al. / Finite Elements in Analysis and Design 60 (2012) 1–12

11

Fig. 8. When t ¼1, the absolute errors between solutions of Problem VI with different number of POD bases for a group of 20 snapshots and the usual LSMFE formulation Problem V with piecewise 1st degree polynomial vectors (on left hand side) and piecewise 1st degree polynomials (on right hand side).

process will be greatly saved, which has also shown that finding the approximate solutions for parabolic equations with the reduced-order LSMFE formulation Problem VI is computationally very effective. And the results above two numerical examples are consistent with those obtained for the theoretical cases.

6. Conclusions and perspective In this paper, we have employed the POD method to establish a reduced-order LSMFE formulation for parabolic equations, analyzed the errors between the solutions of their usual LSMFE formulation and solutions of the POD reduced-order LSMFE formulation, provided the implementation of algorithm for solving reduced-order LSMFE formulation, and discussed theoretically the relation between the number of snapshots and the number of solutions at all time instances, which have shown that our present method has improved and innovated the existing methods (see, e.g., [26–43]) and has supplied scientific theoretic basis for service of practical applications. We have validated the correctness of our theoretical results with two numerical examples. Though snapshots and POD basis of our numerical examples are constructed with the first L ðL 5NÞ solutions of the usual LSMFE formulation, when one computes actual problems, this process can be omitted in actual applications and one may construct the snapshots and POD basis with interpolation or data assimilation by drawing samples from experiments, then solve Problem VI, while it is unnecessary to solve Problem V. Thus, the time-consuming calculations and resource demands in the computational process are greatly saved and the computational efficiency is vastly improved. Therefore, the method in this paper holds a good prospect of extensive applications. In this paper, we only use the forward LSMFE formulation based POD to deal with parabolic equations. However, to solve an inverse problem of parabolic equations, for example, to find the initial conditions, boundary conditions, source term, coefficients (if need), and discuss the POD basic sensitivity of initial condition by using existing data with the POD technique is a very interesting work and an important applied topic of POD, which is our future research work. Future another research work in this area will aim at extending the reduced-order LSMFE formulation, applying it to a realistic atmospheric operational forecast system and to a set of more complicated PDEs such as the atmosphere quality forecast system, the ocean fluid forecast system, and so on.

References [1] V. Thome´e, Galerkin Finite Element Methods for Parabolic Problems, Springer, Berlin, 1997.

[2] Z.D. Luo, Mixed Finite Element Methods and Applications, Chinese Science Press, Beijing, 2006. [3] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, SpringerVerlag, New York, 1991. [4] F. Brezzi, J. Douglas Jr., Stabilized mixed method for the Stokes problem, Numer. Math. 53 (1988) 225–235. [5] P.B. Bochev, M.D. Gunzburger, Least-Squares Finite Element Methods, Springer, New York, 2009. [6] J. Douglas Jr., J.P. Wang, An absolutely stabilized finite element method for the Stokes problem, Math. Comput. 52 (1989) 495–508. [7] P. Hansbo, A. Szepessy, A velocity–pressure streamline diffusion finite element method for the incompressible Navier–Stokes equations, Comput. Methods Appl. Mech. Eng. 84 (1990) 175–192. [8] 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. [9] T.J. Hughes, L.P. France, M. Balestra, A new finite element formulation for computational fluid dynamics, V. Circumventing the Bubuka–Brezzi condition: a stable Petrov–Galerkin formulation of the Stokes problem accommodating equal-order interpolation, Comput. Methods Appl. Mech. Eng. 59 (1986) 85–99. [10] T.J. Hughes, T.E. Tezduyar, Finite element methods for first-order hyperbolic systems with particular emphasis on the compressible Euler equations, Comput. Methods Appl. Mech. Eng. 45 (1984) 217–284. [11] C. Johnson, J. Saranen, Streamline diffusion methods for the incompressible Euler and Navier–Stokes equations, Math. Comput. 47 (1986) 1–18. [12] P. Sun, Z.D. Luo, J. Chen, Petrov least squares mixed finite element method for the non stationary conduction convection problems, Math. Numer. Sin. 31 (1) (2009) 87–98. [13] T.X. Zhou, M.F. Feng, A least squares Petrov–Galerkin finite element method for the stationary Navier–Stokes equations, Math. Comput. 60 (1993) 531–543. [14] P. Holmes, J.L. Lumley, G. Berkooz, Turbulence, Coherent Structures, Dynamical Systems and Symmetry, Cambridge University Press, Cambridge, UK, 1996. [15] K. Fukunaga, Introduction to Statistical Recognition, Academic Press, New York, 1990. [16] I.T. Jolliffe, Principal Component Analysis, Springer-Verlag, Berlin, 2002. [17] 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 Dynamics 192 (1988) 115–173. [18] J.L. Lumley, Coherent Structures in Turbulence, in: R.E. Meyer (Ed.), Transition and Turbulence, Academic Press, 1981, pp. 215–242. [19] H.V. Ly, H.T. Tran, Proper orthogonal decomposition for flow calculations and optimal control in a horizontal CVD reactor, Q. Appl. Math. 60 (2002) 631–656. [20] A.J. Majda, I. Timofeyev, E. Vanden-Eijnden, Systematic strategies for stochastic mode reduction in climate, J. Atmos. Sci. 60 (2003) 1705–1723. [21] P. Moin, R.D. Moser, Characteristic-Eddy decomposition of turbulence in channel, J. Fluid Mech. 200 (1989) 417–509. [22] 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. [23] 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. [24] F. Selten, Barophilic empirical orthogonal functions as basis functions in an atmospheric model, J. Atmos. Sci. 54 (1997) 2100–2114. [25] L. Sirovich, Turbulence and the dynamics of coherent structures: Part I–III, Q. Appl. Math. 45 (3) (1987) 561–590. [26] K. Kunisch, S. Volkwein, Galerkin proper orthogonal decomposition methods for parabolic problems, Numer. Math. 90 (2001) 117–148.

12

Z. Luo et al. / Finite Elements in Analysis and Design 60 (2012) 1–12

[27] K. Kunisch, S. Volkwein, Galerkin proper orthogonal decomposition methods for a general equation in fluid dynamics, SIAM J. Numer. Anal. 40 (2) (2002) 492–515. [28] 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. ¨ [29] D. Ahlman, F. Sodelund, J. Jackson, A. Kurdila, W. Shyy, Proper orthogonal decomposition for time-dependent lid-driven cavity flows, Numer. Heat Transfer Part B – Fundamentals 42 (2002) 285–306. [7,8,11,26–34,43,44]. [30] 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. [31] Y.H. Cao, J. Zhu, I.M. Navon, Z.D. Luo, A reduced order approach to fourdimensional variational data assimilation using proper orthogonal decomposition, Int. J. Numer. Methods Fluids 53 (2007) 1571–1583. [32] 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. Methods Biomed. Eng. 27 (2011) 78–94. [33] 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. [34] Z.D. Luo, J. Chen, I.M. Navon, J. Zhu, An optimizing reduced PLSMFE formulation for non-stationary conduction–convection problems, Int. J. Numer. Methods Fluids 60 (4) (2009) 409–436. [35] Z.D. Luo, J. Chen, P. Sun, X.Z. Yang, Finite element formulation based on proper orthogonal decomposition for parabolic problems, Sci. China Ser. A: Math. 52 (3) (2009) 585–596. [36] Z.D. Luo, Z.H. Xie, J. Chen, A reduced MFE formulation based on POD for the non-stationary conduction–convection problems, Acta Math. Sci. 31 (5) (2011) 1765–1785.

[37] 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. Methods Fluids 55 (2) (2007) 143–161. [38] 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. [39] 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. [40] 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. [41] 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. Methods Appl. Mech. Eng. 196 (41–44) (2007) 4184–4195. [42] 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. [43] 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–334. [44] R.A. Adams, Sobolev Spaces, Academic Press, New York, 1975. [45] V. Girault, P.A. Raviart, Finite Element Approximations of the Navier–Stokes Equations, Theorem and Algorithms, Springer-Verlag, New York, 1986. [46] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, Amsterdam, North-Holland, 1978. [47] W. Rudin, Functional and Analysis, second ed., McGraw-Hill Companies, Inc., 1973.