A New Reduced-Order fve Algorithm Based on POD Method for Viscoelastic Equations

Acta Mathematica Scientia 2013,33B(4):1076–1098 http://actams.wipm.ac.cn

A NEW REDUCED-ORDER FVE ALGORITHM BASED ON POD METHOD FOR VISCOELASTIC EQUATIONS∗

)

Hong LI (

School of Mathematical Sciences, Inner Mongolia University, Huhhot 010021, China E-mail: [email protected]



Zhendong LUO (

)†


)

Junqiang GAO (

School of Mathematics and Physics, North China Electric Power University, Beijing 102206, China E-mail: [email protected]; [email protected]

Abstract A proper orthogonal decomposition (POD) technique is used to reduce the finite volume element (FVE) method for two-dimensional (2D) viscoelastic equations. A reduced-order fully discrete FVE algorithm with fewer degrees of freedom and sufficiently high accuracy based on POD method is established. The error estimates of the reducedorder fully discrete FVE solutions and the implementation for solving the reduced-order fully discrete FVE algorithm are provided. Some numerical examples are used to illustrate that the results of numerical computation are consistent with theoretical conclusions. Moreover, it is shown that the reduced-order fully discrete FVE algorithm is one of the most effective numerical methods by comparing with corresponding numerical results of finite element formulation and finite difference scheme and that the reduced-order fully discrete FVE algorithm based on POD method is feasible and efficient for solving 2D viscoelastic equations. Key words proper orthogonal decomposition; finite volume element method; viscoelastic equations; error estimate 2010 MR Subject Classification

1

65N30; 65M10; 76M10

Introduction

Let Ω ⊂ R2 be a bounded domain with piecewise smooth boundary ∂Ω and consider the following initial boundary value problems for two-dimensional (2D) viscoelastic equations on ∗ Received March 26, 2012; revised October 8, 2012. Research of this work was supported by the National Science Foundation of China (12271227, 11061009 and 11061021), Natural Science Foundation of Inner Mongolia (2012MS0106), Science Research Program of Guizhou (GJ[2011]2367), Science Research Program of Inner Mongolia Advanced Education (NJ10006), and Special Funds for Co-construction Project of Beijing and North China Electric Power University. † Corresponding author: Zhendong LUO.

No.4

H. Li et al: A NEW REDUCED-ORDER FVE ALGORITHM

1077

¯ × [0, T ]. Ω Problem I

Find u such that utt − ε∆ut − γ∆u = f, (x, y, t)∈Ω × (0, T ]

(1)

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

(2)

and initial value conditions u(x, y, 0) = ϕ0 (x, y), ut (x, y, 0) = ϕ1 (x, y),

(x, y) ∈ Ω,

(3)

where ut = ∂u/∂t, utt = ∂ 2 u/∂t2 , ε and γ are two positive constants, source term f (x, y, t) and initial value functions ϕ0 (x, y) and ϕ1 (x, y) are all smooth enough to ensure the following analysis validity, and T is the total time. For the sake of convenience and without loss of generality, we may as well suppose that ϕ0 (x, y) and ϕ1 (x, y) are two zero functions in the following theoretical analysis. The viscoelastic equations, i.e., Problem I is an important modeling equation. It can be used to describe many practical problems, for example, the heat conduction with memory materials, nuclear reaction kinetics, viscous elastic mechanics, loose medium pressure, and so on (see [1–7]). However, it usually includes complex computing domain, initial value functions, and source term which are dependent on real practical applied background. Thus, it isn’t usually easy to find their exact solutions for the practical viscoelastic equations; on the contrary, an efficient approach is to find their numerical solutions. The finite volume element (FVE) method (see [8–10]) is considered as one of the most effective numerical methods because of its following advantages. First, it preserves the integral invariants of conservation of mass and that of total energy. Second, it has higher accuracy and is more suitable for computational field with complicated boundary than the finite difference (FD) method. Third, 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 box method (see [11]) or generalized difference method in China (see [12, 13]). The FVE method was widely applied to finding numerical solutions of different types of partial differential equations, for example, second order elliptic equations, parabolic equations, Stokes equations, and Sobolev equation (see [8–15]). Though a semi-discrete FVE formulation [15] and a fully discrete FVE formulation [16] for 2D viscoelastic equations were presented, they include many degrees of freedom too. Thus, an important problem is how to reduce their degrees of freedom and alleviate the computational load as well as save time of calculations and resource demands in the computational process. This is done in a way that guarantees some sufficiently accurate numerical solutions. A proper orthogonal decomposition (POD) was introduced in the context of analysis of turbulent flow by Holmes and Lumley [17]. In other disciplines, the same procedure goes by the names of Karhunen-Loeve decomposition. The POD method was widely and successfully applied to numerous fields, including signal analysis and pattern recognition (see [18]), statistics (see [19]), geophysical fluid dynamics and meteorology (also see [19]). The POD method

1078

ACTA MATHEMATICA SCIENTIA

Vol.33 Ser.B

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 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. In early times, the POD method was mainly used to perform principal component analysis in computations of statistics and search the main behavior of a dynamic system (see [20–25]), until the method of snapshots was introduced by Sirovich (see [26]) 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 were not presented (see [27, 28]), and the singular value decomposition approach combined with POD technique was not used to treat the Burgers equation (see [29]) and the cavity flow problem (see [30]). More recently, some reduced-order FD models and FE (or mixed FE or least-square mixed FE) formulations based on POD method were presented by our research group and are used to treat different types of partial differential equations (see [31–42]), for example, the upper tropical Pacific Ocean model, parabolic problems, Burgers equations, the non-stationary Navier-Stokes equations, the non-stationary conduction-convection problems, and CVD equations. Though a reduced-order FVE formulation for parabolic equations were presented (see for instance, mentioned in [43]) and a reduced-order semi-discrete FVE formulation based on the POD method for 2D viscoelastic equations was posed in [44], to the best of our knowledge, there are no published results addressing the case that the POD method is used to deal with any fully discrete FVE formulation for 2D viscoelastic equations or providing the error estimates between the classical fully discrete FVE solutions and the reduced-order fully discrete FVE solutions. In this paper, we extend the developments in [43] and [44], i.e., combine the POD method with the classical FVE method to establish a reduced-order fully discrete FVE algorithm with lower dimensions and sufficiently high accuracy for 2D viscoelastic equations, derive the error estimates between the reduced-order FVE solutions and the classical FVE solutions, and provide the implementation for solving the reduced-order fully discrete FVE algorithm. Some numerical examples illustrate that the results of numerical computation are consistent with theoretical conclusions and that the reduced-order fully discrete FVE algorithm is one of the most effective numerical methods by comparing with corresponding numerical results of FE formulation and FD scheme. Moreover, it has shown that the reduced-order fully discrete FVE algorithm based on POD method is feasible and efficient for solving 2D viscoelastic equations. The plan of the paper is organized as follows. Section 2 is to derive the classical fully discrete FVE formulation for 2D viscoelastic equations and to generate snapshots from transient solutions computed from the equation system derived by the classical fully discrete FVE formulation. In Section 3, the optimal orthonormal POD bases are reconstructed from the elements of the snapshots with POD method and a reduced-order fully discrete FVE algorithm with lower dimensions and sufficiently high accuracy based on POD method for 2D viscoelastic equations is developed. In Section 4, the error estimates between the classical fully discrete FVE solutions and the reduced-order fully discrete FVE solutions for 2D viscoelastic equations and the implementation for solving the reduced-order fully discrete FVE algorithm are provided. In Section 5, some numerical examples are presented to illustrate that the numerical results

No.4

H. Li et al: A NEW REDUCED-ORDER FVE ALGORITHM

1079

are consistent with previously obtained theoretical conclusions. Moreover, it is shown that the reduced-order fully discrete FVE algorithm is one of the most effective numerical methods by comparing with corresponding numerical results of FE formulation and FD scheme. Section 6 provides main conclusions and future tentative ideas.

2 A Classical Fully Discrete FVE Formulation for Problem I and Generation of Snapshots The Sobolev spaces and their norms used in this context are standard (see [45]). Let U = H01 (Ω). Then, the variational formulation for Problem I can be written as follows. Problem II Find u(t) : [0, T ] → U such that (utt , w) + εa(ut , w) + γa(u, w) = (f, w), ∀w ∈ U

(4)

u(x, y, 0) = ut (x, y, 0) = 0, (x, y) ∈ Ω,

(5)

subject to where a(u, w) = (∇u, ∇w) and (·, ·) denotes the inner product in L2 (Ω) or L2 (Ω)2 . The existence and uniqueness of solution for Problem II are known and there is the following estimate (see [1–3]): kut k0 + k∇ut kL2 (L2 ) + k∇uk0 6 C(kϕ1 k0 + k∇ϕ0 k0 + kf kL2(L2 ) ), p wherek · kH m (H l ) is the normal in space H m (0, T ; H l(Ω)), C = 3 exp(T )/ min{1, 2ε, γ}. Let N be the positive integer, k = T /N denote the time step increment. For any function g(x, y, t), let g n be the semi-discrete approximation with respect to time of g(x, y, t) at time tn = nk (0 6 n 6 N ). Put g n,1/2 =

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

Then, the semi-discrete formulation with respect to time t for Problem II is written as follows. Problem III Find un+1 ∈ U (n = 1, 2, · · · , N − 1) such that ε (∂¯t ∂¯t un , v) + a(∂¯t un + ∂¯t un−1 , v) + kγa(un,1/2 , v) = k(f n,1/2 , v), ∀v ∈ U 2

(6)

subject to u0 (x, y) = u1 (x, y) = 0, (x, y) ∈ Ω,

(7)

where f n,1/2 = [f (x, y, tn+1 ) + f (x, y, tn−1 )]/2. If ϕ0 (x, y) and ϕ1 (x, y) are non-zero functions, it is necessary to define u0 = ϕ0 (x, y) and u1 = ϕ1 (x, y). It has proved in Reference [16] that Problem III has a unique series of solutions un ∈ U (n = 2, · · · , N ) and, if f ∈ W 2,∞ (0, T ; L2(Ω)), there are the following error estimates: ku(tn ) − un k1 6 Ck 2 [kukW 4,∞ (L2 ) + kukW 3,∞ (H 1 ) ],

(8)

where C = 5 max{1, ε}/(12γ). In order to find the fully discrete FVE solutions for Problem II, it is necessary to introduce a FVE approximation for the spatial variables of Problem III.

1080

ACTA MATHEMATICA SCIENTIA

Vol.33 Ser.B

First, let ℑh = {K} be a quasi-uniform triangulation of Ω with h = max hK , where hK is the diameter of the triangle K ∈ ℑh (see [8–13, 46–48]). In order to describe the FVE method we shall introduce a dual partition ℑ∗h based on ℑh whose elements are called the control volumes. We construct the control volume in the same way as in [8,12]. Let zK be the barycenter of K ∈ ℑh . We connect zK with line segments to the midpoints of the edges of K, thus partitioning K into three quadrilaterals Kz (z = (xz , yz ) ∈ Zh (K), where Zh (K) are the S vertices of K). Then with each vertex z ∈ Zh = K∈ℑh Zh (K) we associate a control volume Vz , which consists of the union of the sub-regions Kz , sharing the vertex z. Finally, we obtain a group of control volumes covering the domain Ω, which is called the dual partition ℑ∗h of the triangulation ℑh (see Figure 1). We denote the set of interior vertices of Zh by Zh◦ .

Fig.1

Left chart: A triangle K partitioned into three sub-regions Kz . Right chart: A sample region with dotted lines indicating the corresponding control volume Vz

We call the partition ℑ∗h regular or quasi-uniform, if there exist two positive constants C1 and C2 , being independent of the spatial mesh size h and temporal mesh size, such that C1 h2 6 mes(Vz ) 6 C2 h2 , ∀Vz ∈ ℑ∗h ,

(9)

where mes(Vz ) denotes the measure of dual element Vz . The barycenter-type dual partition can be introduced for any FE triangulation ℑh and leads to relatively simple calculations. Besides, if the FE triangulation ℑh is quasi-uniform, then the dual partition ℑ∗h is also quasi-uniform (see [8, 12, 46]). The trial function space Uh chosen as the linear element space related to ℑh is the set of all the functions wh satisfying the following conditions: (i) wh ∈ C(Ω), wh |∂Ω = 0; (ii) wh |K ∈ P1 , namely wh is a linear function of x and y on each triangular element K ∈ ℑh , determined only by its values on the three vertexes. It is obvious that Uh ⊂ U . For w ∈ U , let Πh w be the interpolation projection of w onto the trial function space Uh . If w ∈ H 2 (Ω), by the interpolation theory of Sobolev spaces (see [8, 12, 46–48]), we have that |w − Πh w|m 6 Ch2−m |w|2 , m = 0, 1,

(10)

where C in this context indicates a positive constant which is possibly different inequality at different occurrences, being independent of the spatial mesh size h and temporal mesh size k. ˜h is chosen as the piecewise constant function space with respect to ℑ∗ The test space U h but is zero on any boundary dual element Vz ∈ ℑ∗h , i.e.,  ˜h = wh ∈ L2 (Ω); wh |Vz ∈ P0 (Vz ) ∀Vz ∈ ℑ∗h , wh |Vz = 0 (Vz ∩ ∂Ω 6= ∅) U

No.4

H. Li et al: A NEW REDUCED-ORDER FVE ALGORITHM

spanned by the following basis functions: for any point z ∈ Zh◦ ,   1, (x, y) ∈ V , z φz (x, y) =  0, elsewhere. ˜h , For any wh ∈ U

wh =

X

wh (z)φz .

1081

(11)

(12)

◦ z∈Zh

˜h , i.e., For w ∈ U , let Π∗h w be the interpolation projection of w onto the test space U Π∗h w =

X

(13)

w(z)φz .

◦ z∈Zh

By the interpolation theory (see [12, 14, 46, 47]) we have that kw − Π∗h wk0 6 Ch|w|1 .

(14)

Moreover, the interpolation projection Π∗h satisfies the following properties (see [12, 14]). Lemma 1 If vh ∈ Uh , then Z (vh − Π∗h vh )dxdy = 0, K ∈ ℑh ; kvh − Π∗h vh kLr (Ω) 6 ChK kvh kW 1 ,r(Ω) , 1 6 r 6 ∞. (15) K

˜h Though the trial function space Uh satisfies Uh ⊂ U like FE methods, the test space U 6⊂ Uh . As in the case of nonconforming FE methods, this is due to the loss of continuity of ˜h on the boundary of two neighboring elements. So the bilinear form a(u, w) the functions in U must be revised accordingly. Resulting from the piecewise integrations by parts on the dual elements Vz , we have that Z X Z w△udxdy = w△udxdy Ω

=

Vz ∈ℑ∗ h

Vz

Vz ∈ℑ∗ h

∂Vz

X Z



 X Z ∂u ∂u wdy − wdx − ∇u∇wdxdy, ∂x ∂y Vz ∗

(16)

Vz ∈ℑh

R where ∂Vz denotes the line integral, with the counter clockwise direction, on the boundary ∂Vz of the dual element. So we have   X Z X Z ∂u ∂u a(u, w) = ∇u∇wdxdy − wdy − wdx , ∀u, w ∈ U. (17) ∂x ∂y Vz ∂Vz ∗ ∗ Vz ∈ℑh

Vz ∈ℑh

˜h is the piecewise constant function space with the characteristic functions of the dual Since U elements Vz as the basis functions, a(u, w) can be revised into   X Z ∂uh ∂uh ah (uh , wh ) = − wh dy − dx ∂x ∂y ∂Vz Vz ∈ℑ∗ h X ˜h , =− wh (z)˜ a(uh , φz ), ∀uh ∈ Uh , ∀wh ∈ U (18) Vz ∈ℑ∗ h

1082

ACTA MATHEMATICA SCIENTIA

Vol.33 Ser.B

R
∂uh h where a ˜(uh , φz ) = ∂Vz ∂u ∂x dy − ∂y dx , z = (xz , yz ). Then, the classical fully discrete FVE formulation for Problem II reads as follows. Problem IV Find un+1 ∈ Uh (n = 1, 2, · · · , N − 1) such that h ε n,1/2 (∂¯t ∂¯t unh , Π∗h vh ) + ah (∂¯t unh + ∂¯t un−1 , Π∗h vh ) + γah (uh , Π∗h vh ) = (f n,1/2 , Π∗h vh ), ∀vh ∈ Uh h 2 (19) subject to u0h (x, y) = u1h (x, y) = 0, (x, y) ∈ Ω. (20) It has proved that Problem IV has a unique set of solutions {unh }N n=2 ⊂ Uh and there are the following results (see [16]). Theorem 2 If the solution u ∈ W 2,3 (Ω) of Problem II, there are the following stabilization of solutions for Problem IV kun+1 k0 + k 2 h

n+1 X i=1

˜ kuih k1 6 Ck

n+1 X

kf (ti )k0 ,

(21)

i=0

where C˜ = nk/ min{C1 , γC2 } 6 T / min{C1 , γC2 }, and the following error estimates between the solution u of Problem II and the solutions unh of Problem IV ku(tn ) − unh k0 + k

n X

ku(ti ) − uih k1 6 Ck (h + k) , n = 1, 2, · · · , N,

(22)

i=0

where C is a constant independent of h but dependent ε, γ, and f of Problem II. Remark 1 If ϕ0 (x, y) and ϕ1 (x, y) are non-zero functions, then it is necessary to let ϕj (x, y) ∈ W 2,3 (Ω) and ujh (x, y) = Πh ϕj (x, y) (j = 0, 1), one may get the same results as Theorem 2. Thus, if only f (x, y, t), ε, γ, ϕ0 (x, y), and ϕ1 (x, y), the triangulation parameter h, and trail function space Uh are given, we can obtain a set of fully discrete solutions unh ∈ Uh (n = 1, 2, · · · , N ) by solving Problem IV. We chose the first L (in general, L ≪ N , for example, L = 20, N = 200) solutions unh (1 6 n 6 L) as snapshots from N solutions unh (1 6 n 6 N ). Remark 2 When one computes actual problems, he may obtain the ensemble of snapshots from physical system trajectories by drawing samples from experiments and interpolation (or data assimilation). For example, when one computes real practical viscoelastic equations, one can use the previous given results to construct the ensemble of snapshots, then reconstruct the POD optimal basis for the ensemble of snapshots by using the following POD method, and finally the trail function space Uh is replaced with the subspace generated with POD basis in order to derive a reduced-order dynamical system with lower dimensions. Thus, the future change of viscoelastic system can be quickly simulated, which is a result of major importance for real-life applications.

3 Generation of POD Basis and Reduced-order FVE Formulation Based on POD Technique for Problem IV For uih (x, y) (i = 1, 2, · · · , L) in Section 2, let Wi (x, y) = uih (x, y) (1 6 i 6 L) and V = span{W1 , W2 , · · · , WL },

(23)

No.4

H. Li et al: A NEW REDUCED-ORDER FVE ALGORITHM

1083

and refer to V as the space generating by the snapshots {Wi }L i=1 at least one of which is assumed to be non-zero vector. Let {ψj }lj=1 denote an orthonormal basis of V with l = dim V. Then each member of the ensemble can be expressed as Wi =

l X

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

(24)

j=1

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

subject to

2 L d X

1 X

min Wi − (Wi , ψj )U ψj

d L {ψj }j=1 U i=1 j=1

(25)

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

(26)

where kWi k2U = k∇uih k20 . A solution {ψj }dj=1 of (25) and (26) is known as a POD basis of rank d. By (24) and orthonormality of ψj , we can rewrite (25) as follows

2

2 L d L X l X

1 X 1 X

Wi −

(W , ψ ) ψ = (W , ψ ) ψ i j U j i j U j

L i=1 L i=1 U U j=1 j=d+1

=

l X

j=d+1



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

(27)

Thus, in order to assure (25) minimum, it is equivalent to find orthonormal basis ψj (j = 1, 2, · · · , l) such that  L d  X 1 X (Wi , ψj )2U max (28) L i=1 {ψj }d j=1 j=1 subject to (ψr , ψj )U = δij , 1 6 r 6 d, 1 6 j 6 r.

(29)

In other words, (25) and (26) are equivalent to look for a vector function ψ, or the so-called POD basis vector, such that most resembles {Wi (x)}L i=1 in mean that it maximizes L

1X 2 |(Wi , ψ)U | subject to (ψ, ψ)U = k∇ψk20 = 1. L i=1

(30)

We cite the idea of snapshots introduced by Sirovich in [26] and choose a special class of trail function for ψ to be of the form: L X ai Wi , (31) ψ= i=1

1084

ACTA MATHEMATICA SCIENTIA

Vol.33 Ser.B

where the coefficients ai (i = 1, 2, · · · , L) are to be determined so that ψ given by the expression (31) provides a maximum for (30). Thus, (30) is equivalent to the eigenvalue problem (see [26, 27, 35]) Av = λv,

(32)

where v = (a1 , a2 , · · · , aL )T , A = (Aik )L×L , and Z 1 Aik = ∇Wi (x, y) · ∇Wk (x, y)dxdy L Ω

(33)

and λ depends on h and τ due to V depending on them. Since the matrix A is a nonnegative Hermitian matrix which has rank l, it has a complete set of orthonormal eigenvectors v 1 = (a11 , a12 , · · · , a1L )T , · · · , a2L )T , , v l = (al1 , al2 , · · · , alL )T

(34)

with the corresponding eigenvalues λ1 > λ2 > · · · > λl > 0. Thus, the solution to the optimization for (30) is given by ψ1 = √

L 1 X 1 aj Wj , Lλ1 j=1

(35)

where a1j (j = 1, 2, · · · , L) are the elements of the eigenvector v 1 corresponding to the largest eigenvalue λ1 . The remaining POD basis elements ψi (i = 2, 3, · · · , l) are obtained by using the elements of other eigenvectors v i (i = 2, 3, · · · , l), i.e., L 1 X i ψi = √ a Wj , i = 2, 3, · · · , l. Lλi j=1 j

(36)

Moreover, the POD basis {ψ1 , ψ2 , · · · , ψl } forms an orthonormal set and there hold the following results (see [26, 27, 35]). Proposition 3 Let λ1 > λ2 > · · · > λl > 0 denote the positive eigenvalues of A and v 1 , v 2 , · · ·, v l the associated orthonormal eigenvectors. Then a POD basis of rank d 6 l is given by 1 ψi = p (W1 , W2 , · · · , WL ) · v i , Lλi

1 6 i 6 d 6 l.

(37)

Furthermore, the following error formula holds

2 L d l X X

1 X

Wi −

= (W , ψ ) ψ λj . i j U j

L i=1 U j=1

(38)

j=d+1

Let U d = span {ψ1 , ψ2 , · · · , ψd }. For any uh ∈ Uh , define a Ritz-projection P d : Uh → U d as follows: (∇P d uh , ∇wd ) = (∇uh , ∇wd ), ∀wd ∈ U d . Then there is an extension P h : U → Uh of P d such that P h |Uh = P d : Uh → U d denoted by (see [49]) (∇P h u, ∇wh ) = (∇u, ∇wh ),

∀wh ∈ Uh ,

(39)

No.4

H. Li et al: A NEW REDUCED-ORDER FVE ALGORITHM

1085

where u ∈ U . Due to (39), the operator P h is well-defined and bounded (see [35, 44]) k∇(P h u)k0 6 k∇uk0 , ∀u ∈ U,

(40)

and the following inequality (see [35, 44, 46]) ku − P h uk0 6 Chk∇(u − P h u)k0 ,

∀u ∈ U.

(41)

And there holds the following result (see [35, 44, 46]). Lemma 4 For every d (1 6 d 6 l), the projection operator P d satisfies L l X 1X k∇(uih − P d uih )k20 6 λj , L i=1

(42)

j=d+1

where uih ∈ V (i = 1, 2, · · · , L) are the solutions of Problem IV. Thus, by using U d , we can obtain a reduced-order fully discrete FVE algorithm based on POD method for Problem IV as follows. Problem V Find un+1 ∈ U d (n = 1, 2, · · · , N − 1) such that d ε n,1/2 (∂¯t ∂¯t und , Π∗h vd ) + ah (∂¯t und + ∂¯t un−1 , Π∗h vd ) + γah (ud , Π∗h vd ) = (f n,1/2 , Π∗h vd ), ∀vd ∈ U d d 2 (43) subject to u0d (x, y) = u1d (x, y) = 0, (x, y) ∈ Ω. (44) Let und = αn1 ψ1 + αn2 ψ2 + · · · + αnd ψd (n = 0, 1, · · · , N ). By using the definition of operator Π∗h , Problem V can be rewritten as follows. Problem VI Find (αn+1 , αn+1 , · · · , αn+1 )T ∈ Rd (n = 1, 2, · · · , N − 1) such that 1 2 d d X

αn+1 [2Bij + k(ε + γk)Cij ] = Gn,n−1 , i = 1, 2, · · · , d, j i

(45)

j=1

where Bij =

P

ψi (xz , yz )

Vz ∈ℑ∗ h

R

Vz

ψj (x, y)dxdy, Gn,n−1 = Fin,n−1 + i

d P

j=1

{4Bij αnj + [(kε − k 2 γ)Cij


R P ψi (xz ,yz ) P ∂ψ (x,y) ∂ψj (x,y) dx − j∂x dy , Fin,n−1 = −2Bij ]αn−1 }, Cij = ψi (xz , yz ) ∂Vz j ∂y 2 ∗ Vz ∈ℑ∗ Vz ∈ℑh h R 0 1 Vz [f (x, y, tn+1 ) + f (x, y, tn−1 )]dxdy, and αi = αi = 0 (i = 1, 2, · · · , d). Remark 3 If ϕ0 (x, y) and ϕ1 (x, y) are non-zero functions, then αrj (j = 1, 2, · · · , d, r = 0, 1) satisfy d X j=1

αrj Bij =

X

Vz ∈ℑ∗ h

ψi (xz , yz )

Z

ϕr (x, y)dxdy, i = 1, 2, · · · , d, r = 0, 1.

(46)

Vz

Remark 4 If ℑh is a uniformly regular triangulation and Uh is the space of piecewise linear function, Problem IV on each time level includes Nh (where Nh is the number of vertices of triangles in ℑh , see [12, 46]) degrees of freedom for , i.e., unknown quantities, while Problem V on each time level does only contain d (d ≪ l 6 L) degrees of freedom. For actual scientific engineering problems, the number of vertices of triangles in ℑh is more than ten of thousands

1086

ACTA MATHEMATICA SCIENTIA

Vol.33 Ser.B

or 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, in Section 5, d = 6, while Nh = 2000 × 2000 = 4 × 106 ). Therefore, Problem V is a reduced-order fully discrete FVE algorithm based on POD method for Problem IV. Moreover, since the development and change of future physical phenomena are closely related to previous given results, one may truly capture the laws of the development and change of future physical phenomena by using existing results as snapshots. Therefore, the POD method provides an useful and important application.

4 Error Estimates of Solutions for Problem V and Implementation of Its Algorithm 4.1

Error Estimates of Solutions for Problem V In this subsection, we recur to the classical FVE method to derive the error estimates of solutions for Problem V. To this end, it is necessary to introduce some preparative lemmas. We first introduce the following discrete Gronwall lemma (see [46, 50]). Lemma 5 (Discrete Gronwall Lemma) If {an }, {bn }, and {cn } are three positive sequences, and {cn } is monotone, that satisfy ¯ a n + b n 6 cn + λ

n−1 X

¯ > 0; ai , λ

a 0 + b 0 6 c0 ,

i=0

¯ then an + bn 6 cn exp(nλ), n > 0. By using the same as approaches of proof in [12] or [14], we yield the following two lemmas. Lemma 6 The bilinear form a(uh , Π∗h vh ) is symmetric, bounded, and positive definite, i.e., a(uh , Π∗h vh ) = a(vh , Π∗h uh ) = a(uh , vh ), ∀uh , vh ∈ Uh , and there are two positive constants C˜1 and C˜2 such that a(uh , vh ) 6 C˜2 k∇uh k0 k∇vh k0 , a(uh , uh ) > C˜1 k∇uh k20 , ∀uh , vh ∈ Uh . Lemma 7

The following statement holds: (uh , Π∗h u ¯h ) = (¯ uh , Π∗h uh ), ∀uh , u¯h ∈ Uh .

For any u ∈ H m (Ω)2 (m = 0, 1) and vh ∈ Uh , |(u, vh ) − (u, Π∗h vh )| 6 Chm+s kukmkvh ks , s = 0, 1. Set |kuh k|0 = (uh , Π∗h uh )1/2 , then |k·k|0 is equivalent to k·k0 on Uh , i.e., there exist two positive constants C3 and C4 such that C3 kuh k0 6 |kuh k|0 6 C4 kuh k0 ,

∀uh ∈ Uh .

We have the following main results for Problem V. Theorem 8 Under hypotheses of Theorem 2, Problem V has a unique set of solution und ∈ U d such that kun+1 k0 d

+k

2

n+1 X i=1

kuid k1

6 Ck

n+1 X i=0

kf (ti )k0 ,

1 6 n 6 N,

(47)

No.4

H. Li et al: A NEW REDUCED-ORDER FVE ALGORITHM

1087

where C = nk/ min{C˜1 , αγC3 } 6 T / min{C˜1 , αγC3 } and α = min{2C3 , kε + k 2 γ}. And if k = O(h) and L2 = O(N ), then the following error estimates with respect to the norms in H01 (Ω) and L2 (Ω) hold kunh

−

und k0

+k

n X

k∇(uih

−

i=0

uid )k0

 1/2 l X 1/2 6C k λj , n = 1, 2, · · · , L.

(48)

j=d+1

Proof Noting that U d ⊂ Uh . For fixed k, let Ad (un+1 , vd ) = 2(un+1 , Π∗h vd ) + kεa(un+1 , d d d n+1 n−1 n−1 ∗ 2 ∗ n ∗ 2 n,1/2 ∗ ∗ Πh vd )+k γa(ud , Πh vd ), Fd (vd ) = (2ud −ud , Πh vd )+2k (f , Πh vd )+kεa(ud , Πh vd )− ∗ k 2 γa(un−1 , Π v ). Thus, Problem V is rewritten as follows. h d d Problem VII Find un+1 ∈ U d (n = 1, 2, · · · , N − 1) such that d Ad (un+1 , vd ) = Fd (vd ), d

∀vd ∈ U d

(49)

u0d (x, y) = u1d (x, y) = 0, (x, y) ∈ Ω.

(50)

subject to We get by Lemmas 6 and 7 that Ad (vd , vd ) = 2k|vd |k20 + (kε + k 2 γ)k∇vd k20 > αkund k21 ,

∀vd ∈ U d ,

(51)

where α = min{2C3 , kε+k 2 γ}, i.e., the left hand of (49) is a symmetric positive definite bilinear function on U d , and Ad (·, ·) and Fd (·, ·) are bounded obviously. Thus, by Lax-Milgram (see [46–48]), Problem VII, i.e., Problem V has a unique set of solution und ∈ U d (n = 2, 3, · · · , N ). Let Qnd = un+1 − und . Taking vd = un+1 − un−1 = un+1 − und + und − un−1 = Qnd + Qn−1 in d d d d d d Problem V yields that (Qnd − Qn−1 , Π∗h (Qnd + Qn−1 )) + d d

kε a(Qnd + Qn−1 , Π∗h (Qnd + Qn−1 )) d d 2

k2 γ a(un+1 + un−1 , Π∗h (un+1 − un−1 )) = k 2 (f n,1/2 , Π∗h (Qnd + Qn−1 )), 1 6 n 6 N − 1, d d d d d 2 Q1h = Q0h = 0. (52) +

By using Lemmas 6 and 7, H¨ older inequality, and Cauchy inequality, we get that k2 γ kε k∇(Qnd + Qn−1 [k∇un+1 k20 − k∇un−1 )k20 + k20 ] d d d 2 2   k2 [kf (tn+1 )k0 + kf (tn−1 )k0 ] k|Qnd − Qn−1 |k0 . )) 6 = k 2 (f n,1/2 , Π∗h (Qnd + Qn−1 d d 2 k|Qnd |k20 − k|Qn−1 |k20 + d

(53)

Noting that k|Qnd − Qn−1 |k0 6 k|Qnd |k0 + k|Qn−1 |k0 or k|Qnh − Qn−1 |k0 6 k|un+1 |k0 + k|un−1 |k0 . d d h h h We obtain from (53) that k|Qnd |k0 − k|Qn−1 |k0 + d

 k2 αk 2 γ  k∇un+1 k0 − k∇un−1 k0 6 [kf (tn+1 )k0 + kf (tn−1 )k0 ] . (54) d h 2 2

Noting again that Q0h = 0 and u1h = u0h = 0. Summing (54) from 1 to n yields that k|Qnh |k0

n+1 X  αk 2 γ  n+1 n 2 + kuh k1,h + kuh k1,h 6 k kf (ti )k0 . 2 i=0

(55)

1088

ACTA MATHEMATICA SCIENTIA

Vol.33 Ser.B

Since k|Qnh |k0 = k|un+1 − unh |k0 > k|un+1 |k0 − k|unh |k0 , summing (55) from 1 to n yields that h h k|un+1 |k0 + k 2 αγ h

n+1 X

kuih k1,h 6 k 2

i=1

j+1 n+1 XX

kf (ti )k0 6 nk 2

j=0 i=0

n+1 X

kf (ti )k0 ,

(56)

i=0

which yields (47). Since U d ⊂ Uh , subtracting Problem V from Problem IV taking vh = vd ∈ U d yields that ε (∂¯t ∂¯t (unh − und ), Π∗h vd ) + a(∂¯t (unh − und ) + ∂¯t (un−1 − un−1 ), Π∗h vd ) h d 2 n,1/2 n,1/2 +kγa(uh − ud , Π∗h vd ) = 0, ∀vd ∈ U d . (57) Let θn = unh − und = (unh − P d unh ) + (P d unh − und ) = ρn + ̺n . On the one hand, by using Lemmas 6 and 7, we have that ε (∂¯t ∂¯t θn , Π∗h (θn+1 − θn−1 )) + a(∂¯t θn + ∂¯t θn−1 , Π∗h (θn+1 − θn−1 )) 2 γ + a(θn+1 + θn−1 , Π∗h (θn+1 − θn−1 )) 2 ε ¯ = (∂t θn − ∂¯t θn−1 , ∂¯t θn + ∂¯t θn−1 ) + a(θn+1 − θn−1 , θn+1 − θn−1 ) 2k γ + a(θn+1 + θn−1 , θn+1 − θn−1 ) 2 ε γ = k|∂¯t θn |k20 − k|∂¯t θn−1 |k20 + k∇(θn+1 − θn−1 )k20 + [k∇θn+1 k20 − k∇θn−1 k20 ] 2k 2 γ n+1 2 n 2 n−1 2 n−1 2 ¯ ¯ k0 − k∇θ > k|∂t θ |k0 − k|∂t θ k0 ]. |k0 + [k∇θ (58) 2 On the other hand, from the definition of P d , Lemmas 6 and 7, H¨older inequality, error equation (57), and (41), we obtain if h = O(k) that ε (∂¯t ∂¯t θn , Π∗h (θn+1 − θn−1 )) + a(∂¯t θn + ∂¯t θn−1 , Π∗h (θn+1 − θn−1 )) 2 γ + a(θn−1 + θn−1 , Π∗h (θn+1 − θn−1 )) 2 ε = (∂¯t ∂¯t θn , Π∗h (ρn+1 − ρn−1 )) + a(∂¯t θn + ∂¯t θn−1 , Π∗h (ρn+1 − ρn−1 )) 2 γ + a(θn+1 + θn−1 , Π∗h (ρn+1 − ρn−1 )) 2 ε a(ρn + ̺n − ρn−1 − ̺n−1 , Π∗h (ρn+1 − ρn−1 )) = (∂¯t ∂¯t θn , Π∗h (ρn+1 − ρn−1 )) + 2k γ + a(ρn+1 + ρn−1 + ̺n+1 + ̺n−1 , Π∗h (ρn+1 − ρn−1 )) 2 ε a(ρn − ρn−1 , ρn+1 − ρn−1 ) = (∂¯t ∂¯t θn , Π∗h (ρn+1 − ρn−1 )) + 2k γ + a(ρn+1 + ρn−1 , ρn+1 − ρn−1 ) 2 6 Ck[k|∂¯t θn |k20 + k|∂¯t θn−1 |k20 ] + C[k −2 k|ρn+1 |k20 + k −2 k|ρn−1 |k20 + k∇ρn k20 +k∇ρn+1 k20 + k∇ρn−1 k20 ] 6 Ck[k|∂¯t θn |k2 + k|∂¯t θn−1 |k2 + C(k∇ρn+1 k2 + k∇ρn−1 k2 ). 0

0

0

0

(59)

Combining (58) with (59) yields that γ k|∂¯t θn |k20 − k|∂¯t θn−1 |k20 + [k∇θn+1 k20 − k∇θn−1 k20 ] 2 6 Ck[k|∂¯t θn |k20 + k|∂¯t θn−1 |k20 + C(k∇ρn+1 k20 + k∇ρn−1 k20 ).

(60)

No.4

1089

H. Li et al: A NEW REDUCED-ORDER FVE ALGORITHM

Summing (60) from 1 to n and using Lemma 4 yield that n l X X γ n 2 n+1 2 n 2 i 2 ¯ ¯ k|∂t θ |k0 + [k∇θ k0 + k∇θ k0 ] 6 Ck k|∂t θ |k0 + CL λj . 2 i=0

(61)

j=d+1

If k is sufficiently small such that Ck 6 1/2, we obtain from (61) that k|∂¯t θn |k20 + γ[k∇θn+1 k20 + k∇θn k20 ] 6 Ck

n−1 X

k|∂¯t θi |k20 + CL

i=0

l X

λj .

(62)

j=d+1

Applying Lemma 5 (Discrete Gronwall Lemma) to (62) yields k|∂¯t θn |k20 + [k∇θn+1 k20 + k∇θn k20 ] 6 CL

l X

λj exp(Ckn) 6 CL

j=d+1

i.e.,

l X

λj ,

(63)

j=d+1

 X 1/2 l k|∂¯t θn |k0 + k∇θn+1 k0 6 C L λj .

(64)

j=d+1

Since k|∂¯t θn |k0 > [k|θn+1 |k0 − k|θn |k0 ]/k, we obtain from (64) that  X 1/2 l k|θn+1 |k0 + kk∇θn+1 k0 6 k|θn |k0 + Ck L λj .

(65)

j=d+1

Summing (65) from 0, 1, to n − 1 yields that kθn k0 + k

n X i=0

 X 1/2 l k∇θi k0 6 Cnk L λj ,

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

(66)

j=d+1

If k = O(N −1 ) = O(L−2 ), from (66) we obtain (48), which completes the proof of Theorem 8. 2 Combining Theorem 2 with 8 yields the following result. Theorem 9 Under hypotheses of Theorem 8, the error estimates between the solutions for Problem II and the solutions for the reduced-order fully discrete Problem V are ku(tn ) −

und k0

+k

n X i=0

k∇(u(ti ) −

uid )k0

 1/2 l X 1/2 6 Ck(h + k) + C k λj , n = 1, 2, · · · , L. j=d+1

Remark 5 Inequality (47) in Theorem 8 shows that the solutions und (n = 1, 2, · · · , N ) to Problem V are stabilized and continuously dependent on source term f (x, y, t) and initial value function ϕ0 (x, y) and ϕ1 (x, y) if ϕ0 (x, y) and ϕ1 (x, y) are two nonzero functions. The condition L2 = O(N ) in Theorem 8 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, for instance in references [27, 28]. Theorems 8 and 9 have respectively provided the error estimates between the solutions of the reducedorder fully discrete FVE algorithm Problem V and the solutions of classical fully discrete FVE formulation Problem IV and Problem II, which can guide choosing the number d of POD

1090

ACTA MATHEMATICA SCIENTIA

Vol.33 Ser.B

1/2 1/2   l l P P λj . Our λj = O(kh) = O(k 2 ) and hk + k 2 + k 1/2 basis such that k 1/2 j=d+1

j=d+1

method here employs first few classical fully discrete FVE solutions unh (n = 1, 2, · · · , L) for Problem IV as assistant analysis. However, when one computes actual problems, one may obtain the ensemble of snapshots from physical system trajectories by drawing samples from experiments with interpolation (or data assimilation) or previous results. Thus, the assistant unh (n = 1, 2, · · · , L) could be replaced with the interpolation functions of experimental or previous results, thus rendering it unnecessary to solve Problem IV, and requiring only to solve directly reduced-order Problem V such that Theorem 8 is satisfied. And then, time instances are continuously extrapolated forward and POD basis is ceaselessly renewed, the rules of development and change of future physical phenomena would be very well simulated. 4.2 Implementation for Solving Problem V In the following, we give the implementation for solving Problem V which consists of the following seven steps. Step 1 Generate the snapshots ensemble Wi (x, y) = uih , i = 1, 2, · · · , L ≪ N, which may be the solutions for Problem IV, 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 )L×L , where Aij = (Wi , Wj )U /L = (∇uih , ∇ujh )/L, and (·, ·) is L2 -inner product. Step 3 Solving the eigenvalue problem Av = λv, v = (a1 , a2 , · · · , aL )T

(67)

obtains eigenvalues λ1 > λ2 > · · · > λl > 0 (l = dim{W1 , W2 , · · · , WL }) and corresponding eigenvectors v j = (aj1 , aj2 , · · · , ajL ) (j = 1, 2, · · · , l). Step 4 For the triangulation parameter h, the time step increment k, and the error δ  1/2 l P needed, decide on the number d of POD basis such that hk + k 2 + k 1/2 λj 6 δ. j=d+1

L P

p aji uih / Lλj (j = 1, 2, · · · , d).

Step 5

Generate POD basis ψj =

Step 6

Solving the following system of equations which only includes d degrees of freedom

i=1

Z

d X

αrj Bij

d X

αn+1 [2Bij + k(ε + γk)Cij ] = Gn,n−1 , i = 1, 2, · · · , d, j i

=

X

ψi (xz , yz )

ϕr (x, y)dxdy, i = 1, 2, · · · , d, r = 0, 1,

Vz

Vz ∈ℑ∗ h

j=1

j=1

where Gn,n−1 = Fin,n−1 + i

d X  4Bij αnj + [(kε − k 2 γ)Cij − 2Bij ]αn−1 , j j=1

Fin,n−1

X ψi (xz , yz ) Z = [f (x, y, tn+1 ) + f (x, y, tn−1 )]dxdy, 2 Vz ∗ Vz ∈ℑh

No.4

1091

H. Li et al: A NEW REDUCED-ORDER FVE ALGORITHM

Bij =

X

ψi (xz , yz )

X

ψi (xz , yz )

Vz ∈ℑ∗ h

ψj (x, y)dxdy,

Vz

Vz ∈ℑ∗ h

Cij =

Z Z

∂Vz



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



obtains (αn1 , αn2 , · · · , αnd )T ∈ Rd (n = 1, 2, · · · , N ) and further yields the solutions und =

d P

i=1

αni ψi

(x, y) (n = 1, 2, · · · , L, L + 1, · · · , N ); Step 7 If kun−1 −und k0 > kund −un+1 k0 (n = L, L+1, · · · , N −1), then und (n = 1, 2, · · · , N ) d d  1/2 l P are the solutions for Problem V, whose the errors are not more than hk+k 2 + k 1/2 λj . j=d+1

Else, i.e., if kun−1 − und k0 < kund − un+1 k0 (n = L, L + 1, · · · , N − 1), let Wi = uid (i = d d n − L, n − L − 1, · · · , n − 1) and repeat Step 1 to Step 6.

5

Some Numerical Experiments

In this section, we present some simple numerical examples of 2D viscoelastic equations to validate the feasibility and efficiency of the reduced-order fully discrete FVE formulation based on POD method. Moreover, it is shown that the reduced-order FVE algorithm is one of the most effective numerical methods by comparing with the corresponding numerical results of FE formulation and FD scheme. 5.1 The Comparison of Reduced-order Fully Discrete FVE Solutions with Classical Fully Discrete FVE Solutions For the sake of convenience and without loss of generality, herein we consider the simple 2D viscoelastic equations as an example, whose computational domain consists of a square with ¯ = [0, 2] × [0, 2]. We take ε = γ = 1/(2π 2 ), the initial function dimensions 2cm ×2cm, i.e., Ω ϕ0 (x, y) = ϕ1 (x, y) = sin πx sin πy, and source term f (x, y, t) = 3 sin πx sin πy exp(t). Though this example is simpler, the ideas and approaches here can be easily applied to numerical computations for other 2D viscoelastic equations with real practical background. ¯ into small squares with side length △x = △y = 0.001cm, and We first divide the Ω then link diagonal of square to divide each square into two triangles on same direction, which √ constitutes a triangularization ℑh with diameter h = 2 × 10−3 and nodes Nh = 4 × 106 . The dual decomposition ℑ∗h is taken as barycenter dual decomposition, i.e., the barycenter of the right triangle K ∈ ℑh is taken as the node of the dual decomposition. In order to satisfy the condition k = O(h) in Theorem 8, we take the time step increment as k = 10−3 . We chose the first L = 20 solutions unh (n = 1, 2, · · · , 20, i.e., at time t = 0.001, 0.002, · · · , 0.02) as snapshots from 200 solutions unh (n = 1, 2, · · · , 200, i.e., at time t = 0.001, 0.002, · · · , 0.2) for the classical fully discrete FVE formulation, i.e., Problem IV. It is shown by computing  1/2 20 P that when d = 6 and k = 10−3 , the error k 1/2 λj 6 4 × 10−6 in Theorem 8, which j=7

guides one needs only to choose 6 POD bases. When we solve the reduced-order fully discrete FVE algorithm, i.e., Problem V based on POD method with 6 optimal POD bases, according seven steps of implementation of algorithm in Subsection 4.2, we find that the reduced-order fully discrete FVE algorithm at t = 2 is still convergent, without renewing POD basis, and

1092

ACTA MATHEMATICA SCIENTIA

Vol.33 Ser.B

the reduced-order fully discrete FVE solution obtained with the reduced-order Problem V is depicted graphically in Fig. 2. The solution obtained with classical fully discrete FVE formulation, i.e., Problem IV is depicted graphically in Fig. 3 at the time level t = 2. Figures 2 and 3 are exhibiting quasi-identical similarity, but the reduced-order fully discrete FVE solution is better than classical fully discrete FVE solution (that is due to use more initial data, i.e., six POD bases which is a open problem).

Fig. 2

Reduced-order fully discrete FVE solution figure when t = 2

Fig. 3

Classical fully discrete FVE solution figure when t = 2

Fig. 4 shows the errors between the POD solutions obtained by the reduced-order fully discrete FVE formulation Problem V with different number of optimal POD bases and the classical solutions obtained by classical fully discrete FVE formulation Problem IV at time t = 2. Comparing the classical fully discrete FVE solutions of Problem IV with the reducedorder fully discrete FVE solutions of Problem V based on POD method containing 6 optimal POD bases implementing numerical simulation for t = 2, we find that for the classical fully discrete FVE formulation Problem IV (including 4 × 106 unknowns) the required computing time is 6 minutes, while for the reduced-order fully discrete FVE formulation Problem V (only including 6 unknowns) with 6 optimal POD bases the corresponding time is only 3 seconds,

No.4

1093

H. Li et al: A NEW REDUCED-ORDER FVE ALGORITHM

i.e., the computing time of the classical fully discrete FVE formulation Problem IV is 120 times that of the reduced-order fully discrete FVE formulation Problem V with 6 optimal bases, while the error between their solutions does not exceed 4 × 10−6 . It is shown that finding the approximate solutions for 2D viscoelastic equations with the reduced-order fully discrete FVE formulation Problem V is computationally very effective. And the results for numerical examples are consistent with those obtained for the theoretical case.

Fig. 4

The errors between the reduced-order fully discrete FVE solutions and the classical fully discrete FVE solutions with different number of POD bases at t = 2

5.2

The Comparisons of FVE Solutions with FE Solutions and FD Solutions Using the same time step increment and spatial step increment, data, and initial and boundary conditions as those in Subsection 5.1, we respectively obtain the solution from the following classical FE formulation with FE space Uh of piecewise linear polynomial. Problem VIII Find un+1 ∈ Uh (n = 1, 2, · · · , N − 1) such that h   (∂¯ ∂¯ un , v ) + εa(∂¯ un + ∂¯ un−1 , v ) + γa(un,1/2 , v ) = (f n,1/2 , v ), t t h t t h h h h

h

h

h

∀vh ∈ Uh ,

 u0 = u1 = Πh sin πx sin πy h h

at the time level t = 2 and the solution from the following reduced-order FE formulation based on POD method with 6 POD basis: Problem IX Find un+1 ∈ U d (n = 1, 2, · · · , N − 1) such that d   (∂¯ ∂¯ un , v ) + εa(∂¯ un + ∂¯ un−1 , v ) + γa(un,1/2 , v ) = (f n,1/2 , v ), t d d d t t d t d

d

d

 u0 = u1 = P d Πh sin πx sin πy d d

d

∀vd ∈ U d ,

at the time level t = 2, which is still convergent without renewing POD basis with the similar seven steps of implementation of algorithm in Subsection 4.2, but the snapshots of Step 1 are substituted with solutions of Problem VIII and Step 6 is substituted with Problem IX, and are respectively depicted graphically in Figures 6 and 5, the required computing time for the reduced-order FE formulation is 6 seconds but the classical FE formulation 9 minutes. Though Figures 5 and 6 are also exhibiting quasi-identical similarity and the reduced-order FE solution (see Fig. 5) is better than the classical FE solution (see Fig. 6), they have more evident

1094

ACTA MATHEMATICA SCIENTIA

Vol.33 Ser.B

dissipation than Figures 2 and 3 and haven’t far better than the reduced-order fully discrete FVE solution and the classical fully discrete FVE solution, respectively.

Fig. 5

Reduced-order FE solution figure when t = 2

Fig. 6

Classical FE solution figure when t = 2

Also using the same time step increment and spatial step increment, data, and initial and boundary conditions as those in Subsection 5.1, we respectively obtain the solution from the following classical implicit FD scheme: Problem X Find un+1 i,j (n = 1, 2, · · · , N − 1) such that " n+1 # n+1 un+1 − 2unij + un−1 ui+1,j − 2un+1 uni,j+1 − 2uni,j + uni,j−1 ij ij i,j + ui−1,j +ε + k2 2∆x2 k 2∆y 2 k " n,1/2 # n,1/2 n,1/2 n,1/2 n,1/2 n,1/2 ui+1,j − 2ui,j + ui−1,j ui,j+1 − 2ui,j + ui,j−1 +γ + ∆x2 ∆y 2 " n−1 # n−1 n−1 n−1 ui+1,j − 2un−1 un−1 i,j + ui−1,j i,j+1 − 2ui,j + ui,j−1 =ε + + f n,1/2 , i = j = 0, 1, 2, · · · , 2000, ∆x2 k ∆y 2 k u0i,j = ϕ0 (xi , yj ), u1i,j = ϕ1 (xi , yj )

No.4

H. Li et al: A NEW REDUCED-ORDER FVE ALGORITHM

1095

at the time level t = 2 and the solution from the reduced-order FD scheme (which is derived with similar to [37]) based on POD method with 6 POD bases at the time level t = 2, which is still convergent without renewing POD basis with the similar seven steps of implementation of algorithm in Subsection 4.2, and are depicted graphically in Figures 8 and 7, respectively. The required computing time for the reduced-order FD formulation is 2.8 seconds but the classical FD formulation 5.6 minutes. Though their computing time is slight less than that of formulations corresponding to FVE method, Figures 7 and 8 are also exhibiting quasi-identical similarity, and the reduced-order FD solution (see Fig. 7) is better than the classical FD solution (see Fig. 8), they have also more evident dissipation than Figures 2 and 3 and haven’t also far better than the reduced-order fully discrete FVE solution and the classical FVE solution, respectively.

Fig. 7

Reduced-order FD solution figure when t = 2

Fig. 8

Classical FD solution figure when t = 2

By comparing the above results of numerical simulations of the reduced-order fully discrete FVE algorithm with classical fully discrete FVE formulation, reduced-order FE formulation, classical FE formulation, reduced-order FD scheme, and classical FD scheme for the 2D vis-

1096

ACTA MATHEMATICA SCIENTIA

Vol.33 Ser.B

coelastic equations, it has shown that the reduced-order fully discrete FVE algorithm is one of the most effective numerical methods.

6

Conclusions and Perspective

In this paper, by starting directly from semi-discrete formulation with respect to time and circumventing semi-discrete FVE formulation with respect to space variable, we have established a reduced-order fully discrete FVE algorithm based on POD technology for 2D viscoelastic equations, discussed theoretically the relation of the number of snapshots and the number of solutions at all time instances, analyzed the errors between the solutions of their classical fully discrete FVE formulation and solutions of the reduced-order fully discrete FVE algorithm, and provided the implementation for solving reduced-order fully discrete FVE algorithm. It has shown that our present method has improved and innovated the existing methods (see, e.g., [27, 28, 44]). Comparing with the theoretical errors, the error estimates have been verified to provide quite good results, namely both the theoretical error and the computing errors coincide within plot accuracy, thus validating both the feasibility and efficiency of our reduced-order fully discrete FVE algorithm. Though snapshots and POD basis of our numerical examples are constructed with the first few solutions of the classical fully discrete FVE formulation, when one computes actual viscoelastic equations, 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 solves Problem V according seven steps of implementation of algorithm in Subsection 4.2, while it is unnecessary to solve Problem IV or it is only necessary to solve Problem V at the first few steps such that the computational load could be alleviated and time-consuming of calculations in the computational process saved. We have also compared the results of numerical simulation of the reduced-order fully discrete FVE formulation with classical fully discrete FVE formulation, reduced-order FE formulation, classical FE formulation, reduced-order FD scheme, and classical FD scheme for the viscoelastic equations and shown that the reduced-order fully discrete FVE algorithm is one of the most effective numerical methods. Therefore, the method in this paper holds a good prospect of extensive applications. Future research work in this area will aim at extending the reduced-order fully discrete FVE formulation, applying it to a realistic viscoelastic equations operational forecast system and to a set of more complicated PDEs such as the atmosphere quality forecast system, the ocean fluid forecast system. References [1] Gurtin M, Pipkin A. A general theory of heat conduction with finite wave speeds. Arch Ration Mech Anal, 1968, 31: 113–126 [2] Yuan Y R. Finite difference method and analysis for three-dimensional semiconductor device of heat conduction. Sci China Ser A, 1996, 11: 21–32 [3] Yuan Y R, Wang H. Error estimates for the finite element methods of nonlinear hyperbolic equations. J Systems Sci Math Sci, 1985, 5(3): 161–171 [4] Lin Y P. A mixed boundary problem describing the propagation of disturbances in viscous media solution for quasi-linear equations. J Math Anal Appl, 1988, 135: 644–653 [5] Suveika K. Mixed problems for an equation describing the propagation of disturbances in viscous media. J Differ Equ, 1982, 19: 337–347

No.4

H. Li et al: A NEW REDUCED-ORDER FVE ALGORITHM

1097

[6] Raynal M. On some nonlinear problems of diffusion in volterra equations//London S, Staffans O. Lecture Notes in Math, 737. Berlin, New York: Springer-Verlag, 1979 [7] Cannon J R, Lin Y. A priori L2 error estimates for finite-element methods for nonlinear diffusion equations with memory. SIAM J Numer Anal, 1990, 27(3): 595–607 [8] Cai Z, McCormick S. On the accuracy of the finite volume element method for diffusion equations on composite grid. SIAM J Numer Anal, 1990, 27(3): 636–655 [9] S¨ uli E. Convergence of finite volume schemes for Poisson’s equation on nonuniform meshes. SIAM J Numer Anal, 1991, 28(5): 1419–1430 [10] Jones W P, Menziest K R. Analysis of the cell-centred finite volume method for the diffusion equation. J Comput Phy, 2000, 165(1): 45–68 [11] Bank R E, Rose D J. Some error estimates for the box methods. SIAM J Numer Anal, 1987, 24(4): 777–787 [12] Li R H, Chen Z Y, Wu W. Generalized Difference Methods for Differential Equations-Numerical Analysis of Finite Volume Methods. Monographs and Textbooks in Pure and Applied Mathematics 226. New York: Marcel Dekker Inc, 2000 [13] Blanc P, Eymerd R, Herbin R. A error estimate for finite volume methods for the Stokes equations on equilateral triangular meshes. Numer Methods Partial Differ Equ, 2004, 20(6): 907–918 [14] Li J, Chen Z X. A new stabilized finite volume method for the stationary Stokes equations. Adv Comput Math, 2009, 30(2): 141–152 [15] Li H R, Luo Z D, Li Q. Generalized difference methods for two-dimensional viscoelastic problems. Chinese J Numer Math Appl, 2007, 29(3): 251–262 [16] Li H, Sun P, Shang Y Q, Luo Z D. A fully discrete finite volume element formulation and numerical simulations for viscoelastic equations. Math Numer Sin, 2012, 34(4): 413–424 [17] Holmes P, Lumley J L, Berkooz G. Turbulence, Coherent Structures, Dynamical Systems and Symmetry. Cambridge: Cambridge University Press, 1996 [18] Fukunaga K. Introduction to Statistical Recognition. New York: Academic Press, 1990 [19] Jolliffe I T. Principal Component Analysis. Berlin: Springer-Verlag, 2002 [20] Moin P, Moser R. Characteristic-eddy decomposition of turbulence in channel. J Fluid Mech, 1989, 200: 471–509 [21] Rajaee M, Karlsson S, Sirovich L. Low dimensional description of free shear flow coherent structures and their dynamical behavior. J Fluid Mech, 1994, 258: 1–29 [22] Joslin R, Gunzburger M, Nicolaides R, Erlebacher G, Hussaini M Y. A self-contained automated methodology for optimal flow control validated for transition delay. AIAA J, 1997, 35: 816–824 [23] Ly H, Tran H. Proper orthogonal decomposition for flow calculations and optimal control in a horizontal CVD reactor. Quart Appl Math, 2002, 60: 631–656 [24] Lumley J. Coherent structures in turbulence//Meyer R E. Transition and Turbulence. New York: Academic Press, 1981 [25] Aubry Y N, Holmes P, Lumley J L, Stone E. The dynamics of coherent structures in the wall region of a turbulent boundary layer. J Fluid Dyn, 1998, 192: 115–173 [26] Sirovich L. Turbulence and the dynamics of coherent structures, Part I-III. Quart Appl Math, 1987, 45: 561–590 [27] Kunisch K, Volkwein S. Galerkin proper orthogonal decomposition methods for parabolic problems. Numerische Mathematik, 2001, 90: 117–148 [28] Kunisch K, Volkwein S. Galerkin proper orthogonal decomposition methods for a general equation in fluid dynamics. SIAM J Numer Anal, 2002, 40: 492–515 [29] Kunisch K, Volkwein S. Control of Burgers’ equation by a reduced order approach using proper orthogonal decomposition. J Optim Theory Appl, 1999, 102: 345–371 [30] Ahlman D, S¨ odelund F, Jackson J, Kurdila A, Shyy W. Proper orthogonal decomposition for timedependent lid-driven cavity flows. Numer Heat Transfer Part B, 2002, 42(4): 285–306 [31] Cao Y H, Zhu J, Luo Z D, Navon I M. Reduced order modeling of the upper tropical pacific ocean model using proper orthogonal decomposition. Comput Math Appl, 2006, 52: 1373–1386 [32] Luo Z D, Ou Q L, Wu J R, Xie Z H. A reduced finite element formulation based on POD for two-dimensional hyperbolic equation. Acta Math Sci, 2012, 32B(5): 1997–2009

1098

ACTA MATHEMATICA SCIENTIA

Vol.33 Ser.B

[33] Luo Z D, Xie Z H, Chen J. A reduced MFE formulation based on POD for the non-stationary conductionconvection problems. Acta Math Sci, 2011, 31B(5) : 1765–1785 [34] Du J, Zhu J, Luo Z D, Navon I M. An optimizing finite difference scheme based on proper orthogonal decomposition for CVD equations. Int J Numer Methods Biom Eng, 2011, 27(1): 78–94 [35] Luo Z D, Chen J, Sun P, Yang X. Finite element formulation based on proper orthogonal decomposition for parabolic equations. Sci China Ser A: Math, 2009, 52(3): 587–596 [36] Luo Z D, Chen J, Zhu J, Navon I M. An optimizing reduced order FDS for the tropical Pacific Ocean reduced gravity model. Int J Numer Meth Fluids, 2007, 55(2): 143–161 [37] Luo Z D, Wang R W, Zhu J. Finite difference scheme based on proper orthogonal decomposition for the non-stationary Navier-Stokes equations. Sci China Ser A: Math, 2007, 50(8): 1186–1196 [38] Luo Z D, Yang X Z, Zhou Y J. A reduced finite difference scheme based on singular value decomposition and proper orthogonal decomposition for Burgers equation. J Comput Appl Math, 2009, 229(1): 97–107 [39] Luo Z D, Zhou Y J, Yang X Z. A reduced finite element formulation based on proper orthogonal decomposition for Burgers equation. Appl Numer Math, 2009, 59(8): 1933–1946 [40] Luo Z D, Zhu J, Wang R W, Navon I M. Proper orthogonal decomposition approach and error estimation of mixed finite element methods for the tropical Pacific Ocean reduced gravity model. Comput Meth Appl Mech Eng, 2007, 196(41-44): 4184–4195 [41] Sun P, Luo Z D, Zhou YJ. Some reduced finite difference schemes based on a proper orthogonal decomposition technique for parabolic equations. Appl Numer Math, 2010, 60: 154–164 [42] Wang R W, Zhu J, Luo Z D, Navon I M. An Equation-Free Reduced Order Modeling Approach to Tropic Pacific Simulation. Advances in Geosciences Book Series. World Scientific Publishing, 2007 [43] Luo Z D, Xie Z H, Shang Y Q, Chen J. A reduced finite volume element formulation and numerical simulations based on POD for parabolic equations. J Comput Appl Math, 2011, 235(8): 2098–2111 [44] Luo Z D, Li H, Zhou Y J, Xie Z H. A reduced FVE formulation based on POD method and error analysis for two-dimensional viscoelastic problem. J Math Anal Appl, 2012, 385(1): 371–383 [45] Adams R A. Sobolev Spaces. New York: Academic Press, 1975 [46] Luo Z D. Mixed Finite Element Methods and Applications. Beijing: Science Press, 2006 [47] Brezzi F, Fortin M. Mixed and Hybrid Finite Element Methods. New York: Springer-Verlag, 1991 [48] Ciarlet P G. The Finite Element Method for Elliptic Problems. Amsterdam, New York: North-Holland, 1978 [49] Rudin W. Functional and Analysis. 2nd ed. McGraw-Hill Companies Inc, 1973 [50] Temam R. Navier-Stokes Equations. 3rd ed. Amsterdam, New York: North-Holland, 1984