Reduced-Order Extrapolation Finite Difference Schemes Based on Proper Orthogonal Decomposition

Chapter 1

Reduced-Order Extrapolation Finite Difference Schemes Based on Proper Orthogonal Decomposition The key objective of this book is to develop the numerical treatments of proper orthogonal decomposition (POD) for partial differential equations (PDEs). With regard to numerical methods for PDEs, the finite difference (FD) method essentially constitutes the basis of all numerical methods for PDEs. In order to introduce how POD works, it is natural to start with the FD method. For the sake of proper self-containedness, in this chapter we first review the basic theory of classical FD schemes. We then introduce the construction, theoretical analysis, and implementations of algorithms for the POD-based reducedorder extrapolation FD (PODROEFD) schemes for the two-dimensional (2D) parabolic equation, 2D nonstationary Stokes equation, and 2D shallow water equation with sediment concentration. Finally, we provide some numerical examples to show what the PODROEFD schemes have over the classical FD schemes. Moreover, it is shown that the PODROEFD schemes are reliable and effective for solving above-mentioned PDEs. The numerical models treated here include both simple equations and coupled systems. The systematic approach we take in this chapter, namely, following the logical sequence of rudiments, modeling equations, error estimatesstability-convergence, POD methods, error estimates for POD solutions, and finally concrete numerical examples, will be the standard for all chapters.

1.1 REVIEW OF CLASSICAL BASIC FINITE DIFFERENCE THEORY 1.1.1 Approximation of Derivative The FD schemes use difference quotients to approximate derivatives. Denote uni,j = u(ix, j y, nt) = u(xi , yj , tn ). Then uni±1,j ±1 = u(xi ± x, yj ± y, tn ), un±1 i,j = u(xi , yj , tn ± t). Derivative approximations have usually the following four forms.

Proper Orthogonal Decomposition Methods for Partial Differential Equations https://doi.org/10.1016/B978-0-12-816798-4.00006-1 Copyright © 2019 Elsevier Inc. All rights reserved.

1

2 Proper Orthogonal Decomposition Methods for Partial Differential Equations

1. Approximation to first-order derivative by a forward difference We have 

∂u ∂x

n

u(xi + x, yj , tn ) − u(xi , yj , tn ) x→0 x

= lim i,j

u(xi + x, yj , tn ) − u(xi , yj , tn ) + O(x) x uni+1,j − uni,j + O(x) = x uni+1,j − uni,j . ≈ x =

(1.1.1)

2. Approximation to first-order derivative by a backward difference We have 

∂u ∂x

n = lim

x→0

i,j

u(xi , yj , tn ) − u(xi − x, yj , tn ) x

u(xi , yj , tn ) − u(xi − x, yj , tn ) + O(x) x uni,j − uni−1,j + O(x) = x n n ui,j − ui−1,j . ≈ x

=

(1.1.2)

3. Approximation to first-order derivative by a central difference We have 

∂u ∂x

n = i,j

≈

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

+ O(x 2 ) (1.1.3)

.

4. Approximation to second derivative by a second-order central difference We have 

∂ 2u ∂x 2

n = i,j

≈

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

+ O(x 2 ) .

(1.1.4)

Reduced-Order Extrapolation Finite Difference Schemes Chapter | 1

3

1.1.2 Difference Operators 1. The definitions of difference operators We denote operators I, Ex , Ex−1 , x , ∇x , δx , μx by the following: I unj = unj ; I is known as the unit operator; Ex unj = unj+1 ; Ex is known as the forward shift operator and denoted by Ex = Ex+1 ; Ex−1 unj = unj−1 ; Ex−1 is known as the backward shift operator and denoted by Ex− = Ex−1 ; x unj = unj+1 − unj ; x is known as the forward difference operator and satisfies x = Ex − I ; ∇x unj = unj − unj−1 ; ∇x is known as the backward difference operator and satisfies ∇x = I − Ex− ; δx unj = unj+ 1 − unj− 1 ; δx is known as the one step central difference and 2

1 2

2

−1

satisfies δx = Ex − Ex 2 ; 1 μx unj = (uj + 1 + uj − 1 ); μx is known as the average operator and satisfies 2 2 2 1 12 − 12 μx = (Ex + Ex ). 2 2. Composite difference operators We have 1 1 i. (μδ)x = (Ex − Ex−1 ) = (x + ∇x ); 2 2 1 − 12 2 2 2 ii. (δx ) = δx δx = (Ex − Ex ) = (Ex − 2I + Ex−1 ); iii. (δx )n = δx (δxn−1 ); nx = x n−1 = · · · ; ∇xn = ∇x · ∇xn−1 = · · · . x 3. Derivative relations with difference operators We have  n x unj ∂u i. + O(x)—forward difference = ∂x j x = =  ii.

∂ 2u ∂x 2

n = j

= =

∇x unj x δx unj x δx2 unj x 2

+ O(x)—backward difference + O(x 2 )—central difference; + O(x 2 )—the second-order central difference

2x unj x 2 ∇x2 unj x 2

+ O(x 2 )—the second-order forward difference + O(x 2 )—the second-order backward difference.

4 Proper Orthogonal Decomposition Methods for Partial Differential Equations

1.1.3 The Formation of Difference Equations 1. Explicit FD schemes The explicit FD scheme implies that the time farthest point values appear only once. For example, un+1 − unj j t is an explicit FD scheme for ficient”.

=

δ (un − 2unj + unj−1 ) x 2 j +1

∂u ∂ 2u = δ 2 , where δ is a positive “diffusion coef∂t ∂x

2. Implicit FD schemes The implicit FD schemes imply that the time farthest point values in the FD scheme appear more than once. For example, un+1 − unj j

δ (un+1 − 2un+1 + un+1 j j −1 ) x 2 j +1

=

t

∂ 2u ∂u =δ 2. ∂t ∂x 3. Semidiscretized difference schemes is an implicit FD scheme for

The semidiscretized difference scheme implies that only the spatial variable is discretized and the time variable is not discretized. For example, 

du dt

 = j

δ (uj +1 − 2uj + uj −1 ) x 2

is a semidiscretized difference scheme for

∂ 2u ∂u =δ 2. ∂t ∂x

1.1.4 The Effectiveness of Finite Difference Schemes 1. The error of an FD scheme For the equation ∂u ∂ 2u =δ 2, ∂t ∂x consider an FD scheme discretized by forward difference for time and central difference for space (FTCS) as follows: − unj un+1 j t

=

δ (un − 2unj + unj−1 ), x 2 j +1

Reduced-Order Extrapolation Finite Difference Schemes Chapter | 1

5

which is expanded at reference point (xj , tn ) by Taylor’s formula into 

∂u ∂ 2u −δ 2 ∂t ∂x

n j



t + 2



∂ 2u ∂t 2

n j

x 2 −δ 12



∂ 4u ∂x 4



n

+ · · · = 0. j

The above equation is known as the modified PDE, where the first parenthesis above is called the source equation (i.e., PDE) and the bracket above is called the remainder (R) or truncation error (TE), denoted by R = TE. The TE is equal to the source equation (PDE) subtracting the FD equation (FDE), i.e., TE = PDE − FDE. The discretization error (DE, i.e., the error of the FD scheme) is equal to the TE plus the boundary error (BE), i.e., DE = TE + BE. Thus, DE = BE + PDE − FDE. The round-off error (ROE) denotes the rounding error of the computing procedure by the computer and the total error is denoted by CE. Then CE = DE + ROE = BE + PDE − FDE + ROE. However, ROE and BE are usually omitted. The error of FD schemes mainly considers TE, which is directly obtained from the approximation of the derivative (1.1.1)–(1.1.4). 2. The consistency of an FD scheme Definition 1.1.1. (1) Let the PDE ut = Lu be discretized by an FD scheme as follows:   αμ un+1 βγ unj+γ . j +μ = μ

(1.1.5)

(1.1.6)

γ

When an FD grid is indefinitely refined and if R = T E satisfies the property that it tends to zero, i.e., lim R = lim T E = 0,

x→0

x→0

t→0

t→0

(1.1.7)

then the FD scheme is said to be compatible with the source equation. t (2) If t = o(x γ ) (γ > 0), i.e., lim = 0, we have R → 0, then the x→0 x γ t→0

FD scheme is said to be compatible with the source equation on the conditions t = o(x γ ) (γ > 0). 3. The stability of an FD scheme

6 Proper Orthogonal Decomposition Methods for Partial Differential Equations

Definition 1.1.2. i. If an error disturbance εjn = u˜ nj − unj is added at certain time − un+1 of solutions level t = tn , i.e., u˜ nj = unj + εjn , and the errors εjn+1 = u˜ n+1 j j and un+1 obtained from the FD scheme (1.1.6) do not produce extra-large u˜ n+1 j j overall growth, i.e., there is a constant K > 0 independent of n and j such that εjn+1   Kεjn ,

or omitting j,

εn+1   Kε n ,

where  ·  represents a norm, then when 0 < K < 1, the FD scheme (1.1.6) is said to be strongly stable, else (when K  1) the FD scheme (1.1.6) is said to be weakly stable. ii. If there is no restriction between the time step and the spatial step in the FD scheme (1.1.6), it is said to be absolutely stable or unconditionally stable. In this case, the time step t can take a larger size. iii. If the stability of the FD scheme (1.1.6) is restricted by some relationship of time step and spatial steps (usually constrained in the form of some inequalities, for example, t = o(x γ ) (γ > 0)), then it is said to be conditionally stable. We have the following criteria for the stability of the FD scheme (1.1.6) (see [192,194]). Theorem 1.1.1. The FD scheme (1.1.6) is stable if and only if there are two positive constants t0 and x0 as well as a nonnegative constant K > 0 independent of n and j such that the solutions of FD scheme (1.1.6) satisfy unj   K,

n = 1, 2, . . . .

Theorem 1.1.2. If the FD scheme (1.1.6) of PDE (1.1.5) satisfies · unj > 0, un+1 j

un+1 j  unj 

< 1, n = 1, 2, · · · ,

(1.1.8)

where  ·  is a discrete norm, then the FD scheme (1.1.6) is stable. 4. Equivalence between stability and convergency of FD schemes Definition 1.1.3. i. Let u∗ (x, t) be an exact solution for PDE (1.1.5) and unj the approximate solutions of the FD scheme (1.1.6) compatible with the PDE (1.1.5). If when t → 0, x → 0 (i.e., grid is infinitely refined), for any sequence (xj , tn ) → (x ∗ , t ∗ ) ∈ , we have lim

x,t→0

unj = u∗ (x ∗ , t ∗ ),

(1.1.9)

then the FD scheme (1.1.6) is said to be convergent. The stability and convergence of FD schemes are their intrinsic properties with the following important equivalence (see [192] or [194, Theorem 1.3.18]).

Reduced-Order Extrapolation Finite Difference Schemes Chapter | 1

7

Theorem 1.1.3. If the PDE (1.1.5) is well posed, i.e., it has a unique solution and is continuously dependent on initial boundary value conditions and compatible with its FD scheme (1.1.6), then the stability of the FD scheme (1.1.6) is equivalent to its convergence. Theorem 1.1.3 is briefly described as “if a PDE is well posed and compatible with its FD scheme, then the stability of the FD scheme (1.1.6) is equivalent to its convergence”. Remark 1.1.1. For linear PDEs, if their FD schemes are stable, then their convergence is ensured. Therefore, it is only necessary to discuss their stability. However, for nonlinear PDEs, because of the complexity, at present we can only investigate nearly similar convergence and local stability analysis instead of global analysis of convergence. 5. Von Neumann’s stability analysis of FD schemes Von Neumann’s stability analysis of FD schemes is also known as the Fourier analysis method. In order to analyze the stability of FD schemes, we first discuss the stability of exact solutions. i. Stability analysis of exact solution The initial value problem ⎧ ∂ ∂ ⎨ ut = L( , 2 , · · · )u, ∂x ∂x ⎩ u(x, 0) = φ(x)

x ∈ R,

(1.1.10)

is said to be well posed, if it has a unique and stable solution. A so-called stable solution means that the solution is continuously dependent on the initial value and can maintain the boundedness of small disturbances. Assume that the initial value φ(x) is periodic and expandable into the following Fourier series: φ(x) =



fk eikx ,

(1.1.11)

k

where the fk s are the Fourier coefficients. Assume that the exact solution u(x, t) of the initial value problem (1.1.10) is also expandable into the following Fourier series:  u(x, t) = Fk (t)eikx , (1.1.12) k

where the Fk (t)s are the Fourier coefficients of the series depending on t.

8 Proper Orthogonal Decomposition Methods for Partial Differential Equations

By substituting (1.1.12) into the first equation of (1.1.10), we obtain 

Fk (t)eikx = L(

k

=

 k

 ∂ ∂ Fk (t)eikx , 2,···) ∂x ∂x k

∂ ∂ L( , 2 , · · · )eikx · Fk (t). ∂x ∂x

(1.1.13)

If L is a linear differential operator, we can rewrite (1.1.13) as follows: 

Fk (t)eikx =



k

L(ik, (ik)2 , · · · )eikx · Fk (t).

(1.1.14)

k

By comparing the coefficients eikx of the LHS and those of the RHS in (1.1.14) and setting them equal, we have Fk (t) = L(ik, (ik)2 , · · · ) · Fk (t).

(1.1.15)

Eq. (1.1.15) has a general solution as follows: Fk (t) = ck eL(ik,(ik)

2 ,··· )t

(1.1.16)

.

condition u(x, 0) = φ(x) of (1.1.10), we obtain

By ikxusing the initial ikx fk e = Fk (0)e , which implies ck = Fk (0) = fk . Thus, (1.1.12) can k

k

be rewritten as follows: u(x, t) =



fk eL(ik,(ik)

2 ,··· )t

· eikx .

(1.1.17)

k

Note that the exact solution u(x, t) is said to be stable, if there is a nonnegative constant M, independent of u, t, and x, such that u(x, t)L2  Mu(x, 0)L2 .

(1.1.18)

Because {eikx } is a set of standard orthogonal bases in L2 (−π, π), we have u(x, t)2L2 = 





| fk eL(ik,(ik)

2 ,··· )t

|2

k

| fk | sup | eL(ik,(ik)

k

u(x, 0)2L2 =

2



| fk |2 ,

2 ,··· )t

|2 ,

(1.1.19) (1.1.20)

k

where u(x, t)L2 =

 π

−π

| u(x, t) |2 dx

1/2

is the norm of u(·, t) in L2 (−π, π).

Reduced-Order Extrapolation Finite Difference Schemes Chapter | 1

9

From (1.1.19)–(1.1.20), we know that | eL(ik,(ik)

2 ,··· )t

| M

(1.1.21)

holds if and only if u(x, t)2L2 = M

2





| fk eL(ik,(ik)

2 ,··· )t

|2

k

| fk |2  M 2 u(x, 0)2L2 .

(1.1.22)

k

Thus, the exact solution u(x, t) of (1.1.10) is stable if and only if there is a nonnegative constant M independent of u, t, and x, such that | eL(ik)t | M. ii. Stability analysis of an FD scheme In the following, we analyze the von Neumann stability conditions of an FD scheme. Let (1.1.10) have the following FD scheme:   αμ un+1 βγ unj+γ . (1.1.23) j +μ = μ

γ

Let the initial time level (n = 0) be denoted as   gk0 eikxj = gk0 eikj x . u0j = φ(xj ) = k

(1.1.24)

k

Let the nth time level solution unj be denoted by unj =

 k

gkn eikxj =



gkn eikj x .

(1.1.25)

k

By inserting (1.1.25) into (1.1.23), and by the standard orthogonality of {eikx }, we obtain gkn+1 = Ggkn , k ∈ Z,

(1.1.26)

where Z is the set of all integers and G is known as the growth factor, which is denoted by the following equation: ⎞ −1 ⎛    αμ eikμx ·⎝ βγ eikγ x ⎠ . G= (1.1.27) μ

γ

Let g = (· · · , g−k , · · · , g−2 , g−1 , g0 , g1 , g2 , · · · , gk , · · · ). We have g n = Gg n−1 , n = 1, 2, · · · .

(1.1.28)

10 Proper Orthogonal Decomposition Methods for Partial Differential Equations

By using Eq. (1.1.28), we have g n 2 = Gg n−1 2 = G2 g n−2 2 = · · · = Gn g 0 2 , n = 1, 2, · · · , (1.1.29) where  · 2 represents the norm in l 2 . Thus, by (1.1.24) and (1.1.25), from (1.1.29), we obtain unj 2  Gn u0j 2 = Gn φ(xj )2 , n = 1, 2, · · · .

(1.1.30)

Then, by Theorem 1.1.1, we obtain the following result. Theorem 1.1.4. The FD scheme (1.1.13) is stable if and only if there are two positive constants t0 and x0 as well as a nonnegative constant K > 0 independent of n and j , such that its growth factor G satisfies Gn   K, n = 1, 2, · · · .

(1.1.31)

Corollary 1.1.5. The FD scheme (1.1.13) is stable if and only if its growth factor G satisfies G  1 + O(t), n = 1, 2, · · · .

(1.1.32)

iii. Von Neumann’s stability analysis of FD schemes The condition of (1.1.32) in Corollary 1.1.5 is usually known as the von Neumann stability condition. It follows that in order to distinguish the stability of the FD scheme (1.1.13), it is necessary to compute its growth factor G by formula (1.1.27) and determine the t and x such that (1.1.32) or G  1 is satisfied. For specific FD schemes, it is easy to compute their growth factor G by (1.1.27). By the standard orthogonality of {eikx }, it is necessary to substitute unj = Gn eikxj , n = 1, 2, · · ·

(1.1.33)

into the FD scheme (1.1.13), and then, eliminating some common factors, one can obtain the growth factor G of (1.1.27). For an FD scheme for two-dimensional linear PDEs, we need only to substitute unj,m = Gn eikxj eikym (n = 1, 2, · · · ) into the FD scheme and then simplify, so we can also obtain the growth factor G. Some relative examples can be found in [192].

1.2 A POD-BASED REDUCED-ORDER EXTRAPOLATION FINITE DIFFERENCE SCHEME FOR THE 2D PARABOLIC EQUATION In this section, we introduce the PODROEFD scheme for the two-dimensional (2D) parabolic equation. The work is based on Luo et al. [111].

Reduced-Order Extrapolation Finite Difference Schemes Chapter | 1

11

For convenience and without loss of generality, let = (ax , bx ) × (cy , dy ) and consider the following 2D parabolic equation. Find u such that ⎧ ⎪ (x, y, t) ∈ × (0, T ), ⎨ ut − u = f, (1.2.1) u(x, y, t) = g(x, y, t), (x, y, t) ∈ ∂ × [0, T ), ⎪ ⎩ u(x, y, 0) = s(x, y), (x, y) ∈ , where f (x, y, t), g(x, y, t), and s(x, y) are given source item, boundary function, and initial value function, respectively, and T is the total time duration. The main motivation and physical background of the parabolic equation are the modeling of heat conduction and diffusion phenomena.

1.2.1 A Classical Finite Difference Scheme for the 2D Parabolic Equation Let x and y be the spatial steps in the x and y directions, respectively, t the time step, and unj,k the function value of u at points (xj , yk , tn ) (0  j  J = [(bx − ax )/x], 0  k  K = [(dy − cy )/y], 0  n  N = [T /t], where [(bx − ax )/x], [(dy − cy )/y], and [T /t] denote the integral parts of (bx − ax )/x, (dy − cy )/y, and T /t, respectively). Thus, the forward difference explicit scheme for the 2D parabolic equation (1.2.1) at reference point (xj , yk , tn ) is given by n un+1 j,k = uj,k +

+

t n (u − 2unj,k + unj−1,k ) x 2 j +1,k

t n n (u − 2unj,k + unj,k−1 ) + tfj,k . y 2 j,k+1

(1.2.2)

For the FD scheme (1.2.2), we have the following stability and convergence. Theorem 1.2.1. If 4t/x 2  1 and 4t/y 2  1, the FD scheme (1.2.2) is stabile. Further, we have the following error estimates: unj,k − u(xj , yk , tn ) = O(t, x 2 , y 2 ), 1  n  N.

(1.2.3)

Proof. If 4t/x 2  1 and 4t/y 2  1, we have   2t 2t t n+1 − | unj,k | + 2 (| unj+1,k | + | unj−1,k |) | uj,k |  1 − x 2 y 2 x t n + (| unj,k+1 | + | unj,k−1 |) + t | fj,k | y 2  un ∞ + tf ∞ , (1.2.4) where  · ∞ is the L∞ ( ) norm. Thus, from (1.2.4), we obtain un+1 ∞  un ∞ + tf ∞ , n = 0, 1, 2, · · · , N − 1.

(1.2.5)

12 Proper Orthogonal Decomposition Methods for Partial Differential Equations

By summing (1.2.4) from 0 to n − 1, we obtain un ∞  u0 ∞ + ntf ∞ , n = 1, 2, · · · , N.

(1.2.6)

Because nt  T , from (1.2.6), we obtain un ∞  s(x, y)∞ + T f ∞ , n = 1, 2, · · · , N,

(1.2.7)

which shows that solutions of the FD scheme (1.2.2) are bounded and continuously dependent on the initial value s(x, y) and source term f (x, y, t). Thus, by Theorem 1.1.1, we deduce that the FD scheme (1.2.2) is stable. Furthermore, the error estimates (1.2.3) are directly obtainable from approximating difference quotients for derivatives. Thus, as long as the time step t, the spatial steps x and y, f (x, y, t), g(x, y, t), and s(x, y) are given, we can obtain classical FD solutions unj,k (0  j  J, 0  k  K, 0  n  N ) for the 2D parabolic equation (1.2.1) by computing the FD scheme (1.2.2).

1.2.2 Formulation of the POD Basis n (1  i  m ≡ (J + 1)(K + 1), i = kJ + k + j + 1, Set uni = unj,k and fin = fj,k 0  j  J, 0  k  K, 0  n  N ). Then, the classical FD solutions for the FD scheme (1.2.2) can be denoted by {uni }N n=1 (1  i  m). We extract the first L (1  i  m, L N ) as snapshots. Further, we sequence of solutions {uni }L n=1 formulate an m × L snapshot matrix ⎛ ⎞ u11 u21 · · · uL 1 ⎜ 1 ⎟ 2 L ⎟ ⎜ u ⎜ 2 u2 · · · u2 ⎟ A=⎜ . . (1.2.8) .. .. ⎟ .. ⎜ . ⎟ . . . ⎠ ⎝ . u1m u2m · · · uL m

By the singular value decomposition, the snapshot matrix A has a factorization    O l×(L−l) l×l (1.2.9) VT, A=U O (m−l)×l O (m−l)×(L−l)  where l×l = diag{σ1 , σ2 , · · · , σl } is a diagonal matrix consisting of the singular values of A according to the decreasing order σ1  σ2  · · ·  σl > 0, U = (ϕ 1 , ϕ 2 , · · · , ϕ m ) is an m × m orthogonal matrix consisting of the orthogonal eigenvectors of AAT , whereas V = (φ 1 , φ 2 , · · · , φ L ) is an L×L orthogonal matrix consisting of the orthogonal eigenvectors of AT A, and O is a zero matrix.

Reduced-Order Extrapolation Finite Difference Schemes Chapter | 1

13

Because the number of mesh points m is much larger than that of snapshots L extracted, the height m of A for matrix AAT is much larger than the width L of A for matrix AT A, but the positive eigenvalues λj (j = 1, 2, · · · , l) of AT A and AAT are identical and satisfy λj = σj2 (j = 1, 2, · · · , l). Therefore, we may first find the eigenvalues λ1  λ2  · · ·  λl > 0 (l = rank A) for matrices AT A and corresponding eigenvectors ϕ j . Then, by the relationship  ϕ j = Aφ j / λj ,

j = 1, 2, . . . , l,

we obtain these eigenvectors ϕ j (j = 1, 2, . . . , l) corresponding to the nonzero eigenvalues for matrix AAT . Take  AM = U





O M×(L−M)

M×M

O (m−M)×M

VT,

O (m−M)×(L−M)

(1.2.10)

 where M×M = diag{σ1 , σ2 , · · · , σM } is the diagonal matrix consisting of the  first M positive singular values of the diagonal matrix l×l . Define the norm of matrix A as A2,2 = supu∈R L Au2 /u2 (where u2 is the l 2 norm for vector u). We have the following. Lemma 1.2.2. Let  = (ϕ 1 , ϕ 2 , · · · , ϕ M ) consist of the first M eigenvectors U = (ϕ 1 , ϕ 2 , · · · , ϕ m ). Then, we have AM =

M 

σi ϕ j φ Ti = T A.

(1.2.11)

i=1

Proof. According to ⎛

AM =

M 

σi ϕ i φ Ti =



ϕ1

· · · ϕM

i=1

⎛

A=

we have



ϕ1

σ1 ⎜ . · · · ϕl ⎜ ⎝ .. 0

··· .. . ···

0 .. . σl

σ1 ⎜ . ⎜ . ⎝ . 0 ⎞⎛ φ T1 ⎟⎜ ⎟ ⎜ .. ⎠⎝ . φ Tl

··· .. . ··· ⎞ ⎟ ⎟, ⎠

⎞⎛ φ T1 0 ⎟⎜ .. ⎟ ⎜ .. . ⎠⎝ . σM φ TM

⎞ ⎟ ⎟, ⎠

14 Proper Orthogonal Decomposition Methods for Partial Differential Equations

⎛ ⎞ σ1 · · · ϕ T1 ⎜ ⎟ ⎜ . T . ⎜ ⎟ ..  A =  ⎜ . ⎝ .. ⎠ ϕ 1 · · · ϕ l ⎝ .. T 0 ··· ϕM ⎛ ⎞⎛ ⎞ σ1 · · · 0 φ T1 ⎟⎜ ⎟  ⎜ . . ⎟⎜ . ⎟ .. =  IM O ⎜ . .. ⎠ ⎝ .. ⎠ ⎝ .. 0 · · · σl φ Tl ⎛ ⎞ ⎞⎛ σ1 · · · 0 · · · 0 φ T1 ⎜ ⎟ ⎟⎜ . .. .. . ⎟⎜ . ⎟ = ⎜ . ⎝ .. . · · · .. ⎠ ⎝ .. ⎠ 0 · · · σM · · · 0 φ Tl ⎛ ⎞⎛ σ1 · · · 0 φ T1 ⎟⎜  ⎜ . .. ⎟ ⎜ .. .. = ϕ1 · · · ϕM ⎜ . ⎝ .. . ⎠⎝ . 0 · · · σM φT ⎛

⎞⎛ φ T1 0 ⎟⎜ .. ⎟ ⎜ .. . ⎠⎝ . σl φ Tl

⎞ ⎟ ⎟ ⎠

⎞ ⎟ ⎟ = AM . ⎠

M

Thus, by the relationship between the matrix norm and its spectral radius, we have  min A − B2,2 = A − AM 2,2 = A − T A2,2 = λM+1 , rank(B)M

(1.2.12) where  = (ϕ 1 , ϕ 2 , · · · , ϕ M ) consists of the first M eigenvectors U = (ϕ 1 , ϕ 2 , · · · , ϕ m ). If the L column vectors of A are denoted by un = (un1 , un2 , · · · , unm )T (n = 1, 2, · · · , L), we have  un − unM 2 = (A − T A)ε n 2  A − T A2,2 ε n 2 = λM+1 , (1.2.13)

n n where unM = M j =1 (ϕ j , u )ϕ j represents the projection of u onto  = n n (ϕ 1 , ϕ 2 , · · · , ϕ M ), (ϕ j , u ) is the inner product of ϕ j and u , and ε n denotes the unit vector with the nth component being 1. The inequality (1.2.13) √ shows that unM is an optimal approximation of un whose error is no more than λM+1 . Thus,  is just an orthonormal optimal POD base of A.

1.2.3 Establishment of the POD-Based Reduced-Order Finite Difference Scheme for the 2D Parabolic Equation We still denote the classical FD solution vectors from the FD scheme (1.2.2) by un = (un1 , un2 , · · · , unm )T (n = 1, 2, · · · , L, L + 1, · · · , N ). Thus, we can write

Reduced-Order Extrapolation Finite Difference Schemes Chapter | 1

15

the FD scheme (1.2.2) in vector form as follows: un+1 = un +

t t Bun + Cun + tF nm , x 2 y 2

(1.2.14)

where F nm = (f1n , f2n , · · · , fmn )T , and ⎛

−1 1 0 0 .. . 0 0

1 0 −2 1 1 −2 0 1 .. .. . . 0 0 0 0

⎜ −1 ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ C=⎜ ⎜ 1 ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝

   ··· ···

⎜ ⎜ ⎜ ⎜ ⎜ ⎜ B =⎜ ⎜ ⎜ ⎜ ⎜ ⎝

⎛

0 0 1 −2 .. . 0 0

1

−2 −2 ..

.

. ..

0 0 0 0 .. . −2 1

0 0 0 0 .. . 1 −1

⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

,

m×m

⎞

J times 0

..

··· ··· ··· ··· .. . ··· ···

. 1

⎟ ⎟ ⎟ .. ⎟ . ⎟ ⎟ ⎟ .. ⎟ . ⎟ ⎟ ⎟ 1 ⎟ ⎟ ⎟ ⎟ −2 ⎟ ⎟ ⎟ ⎟ −2 ⎟ ·· ·  · ·· −1 ⎠ J times 0

.

m×m

In order to estimate a tridiagonal matrix, it is necessary to introduce the following lemma (see [192, Theorems 1.3.1 and 1.3.2]). Lemma 1.2.3. The tridiagonal matrix ⎛ d b 0 0 ⎜ ⎜ c a b 0 ⎜ ⎜ 0 c a b ⎜ ⎜ B˜ = ⎜ 0 0 c a ⎜ . . . . ⎜ . . . . ⎜ . . . . ⎜ ⎝ 0 0 0 0 0 0 0 0

··· ··· ··· ··· .. . ··· ···

0 0 0 0 .. . a c

0 0 0 0 .. . b d

⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ m×m

16 Proper Orthogonal Decomposition Methods for Partial Differential Equations

has the following eigenvalues: √ λ˜ j = a + 2 bc cos[(2j − 1)π/(2m + 1)],

j = 1, 2, · · · , m.

Therefore, by using the relationship between the matrix eigenvalue and its norm, we have B2,2 = C2,2 =| 2 − 2 cos[(2m − 1)π/(2m + 1)] | =| 2 − 2 cos[π − 2π/(2m + 1)] | =| 2 + 2 cos[2π/(2m + 1)] |< 4.

(1.2.15)

If we replace un of (1.2.14) with u∗n = T Aεn (n = 1, 2, · · · , L) and = α n (n = L + 1, L + 2, · · · , N ), we obtain the PODROEFD scheme as follows: ⎧ ⎪ u∗n = T Aεn , n = 1, 2, · · · , L, ⎪ ⎪ ⎨ t t (1.2.16) Bα n + Cα n + tF nm , α n+1 = α n + ⎪ 2 x y 2 ⎪ ⎪ ⎩ n = L, L + 1, · · · , N − 1,

u∗n

n )T are vectors yet to be determined. where α n = (α1n , α2n , · · · , αM By using the orthogonal vectors in T multiplied by Eq. (1.2.16), we obtain ⎧ n T n ⎪ ⎪ ⎨ α =  u , n = 1, 2, · · · , L, t t T T T n n n (1.2.17) α n+1 = α n + x 2  Bα + y 2  Cα + t F m , ⎪ ⎪ ⎩ n = L, L + 1, · · · , N − 1.

After α n (n = L, L+1, · · · , N ) are obtained from the system of Eq. (1.2.17), we can obtain the PODROEFD solution vectors for Eq. (1.2.1) as follows: u∗n = α n , n = 1, 2, · · · , L, L + 1, · · · , N.

(1.2.18)

Further, we can obtain the PODROEFD solution components for Eq. (1.2.1) as follows: ∗n u∗n j,k = ui , 0  j  J, 0  k  K,

1  i = k(J + 1) + j + 1  m = (K + 1)(J + 1).

(1.2.19)

Remark 1.2.1. It is easily seen that the classical FD scheme (1.2.2) at each time level contains m unknown quantities, whereas the system of Eqs. (1.2.17)– (1.2.18) at the same time level (when n > L) contains only M unknown quantities (M L m, usually M = 6, but m = O(104 ) ∼ O(106 )). Therefore, the PODROEFD model (1.2.17)–(1.2.18) includes very few degrees of freedom and does not involve repeated computations. Here, we extract the snapshots from the

Reduced-Order Extrapolation Finite Difference Schemes Chapter | 1

17

first L classical FD solutions; but when we solve real-world problems we may, instead, extract snapshots from the samples of experiments of physical system trajectories.

1.2.4 Error Estimates of the Reduced-Order Finite Difference Solutions for the 2D Parabolic Equation First, from (1.2.13), we obtain un − u∗n 2 = un − unM 2 = (A − T A)ε n 2 

 λM+1 ,

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

(1.2.20)

Rewrite the second equation of the system of Eqs. (1.2.16) as follows: t t Bu∗n + Cu∗n + tF nm , 2 x y 2 n = L, L + 1, · · · , N − 1.

u∗n+1 = u∗n +

(1.2.21)

Put δ = tB2,2 /x 2 + tC2,2 /y 2 . By subtracting (1.2.21) from (1.2.14) and taking norm, from (1.2.20), we have un+1 − u∗n+1 2  (1 + δ)un − u∗n 2  ···  (1 + δ)n+1−L uL − u∗L 2   (1 + δ)n+1−L λM+1 , n = L, L + 1, · · · , N − 1.

(1.2.22)

By summarizing the above discussion and noting that the absolute value of each vector component is less than the vector norm, from (1.2.3), we obtain the following theorem. Theorem 1.2.4. The errors between the solutions unjk from the classical FD scheme and the solutions u∗n j k from the PODROEFD scheme (1.2.17)–(1.2.18) satisfy the following estimates:  | unjk − u∗n j k | E(n) λM+1 , 1  j  J, 1  k  K, 1  n  N, (1.2.23) where E(n) = 1 (1  n  L), whereas E(n) = (1 + δ)n−L (L + 1  n  N ), δ = tB2,2 /x 2 + tC2,2 /y 2 . Further, the errors between the accurate solution u(x, y, t) from Eq. (1.2.1) and the solutions u∗n j k from the POD-based reduced-order FD scheme (1.2.17)–(1.2.18) satisfy the following estimates:  2 2 | u(xj , yk , tn ) − u∗n j k |= O(E(n) λM+1 , t, x , y ), 1  j  J, 1  k  K, 1  n  N.

(1.2.24)

18 Proper Orthogonal Decomposition Methods for Partial Differential Equations

√ Remark 1.2.2. The error terms containing λM+1 in Theorem 1.2.4 are caused by the POD-based reduced-order for the classical FD scheme, which could be used to select the number of the POD basis, i.e., it is necessary to √ take M such that λM+1 = O(t, x 2 , y 2 ). Whereas E(n) = (1 + δ)n−L (L + 1  n  N ) are caused by extrapolating iterations, which√could be used as the guide for renewing the POD basis, i.e., if (1 + δ)n−L λM+1 > necessary to update the POD basis. If we take λM+1 max(t, x 2 , y 2 ), it is √ that satisfies (1 + δ)N−L λM+1 = O(t, x 2 , y 2 ), then the PODROEFD scheme (1.2.17)–(1.2.18) is convergent and, thus, we do not have to update the POD basis.

1.2.5 The Implementation of the Algorithm of the POD-Based Reduced-Order Finite Difference Scheme for the 2D Parabolic Equation In order to facilitate the use of the PODROEFD scheme for the 2D parabolic equation, the following implementation steps of the algorithm for the PODROEFD scheme (1.2.17)–(1.2.18) are helpful. Step 1. Classical FD computation and extraction of snapshots Write the classical FD scheme (1.2.2) in vector form (1.2.14) and find the solution vectors un = (un1 , un2 , · · · , unm )T (n = 1, 2, · · · , L) of (1.2.14) at the first few L steps (in the following, say, take L = 20 in Section 1.2.6). Step 2. Snapshot matrix A and eigenvalues of AT A Formulate snapshot matrices A = (uni )m×L and compute the eigenvalues λ1  ˜ of λ2  · · ·  λl > 0 (l = rankA) and the eigenvectors φ j (j = 1, 2, · · · , M) T matrices A A. Step 3. Choice of POD basis For the error tolerance μ = O(t, x 2 , y 2 ), decide the numbers M (M  ˜ M) of POD bases such that λu(M+1)  μ and formulate the POD bases  =  (ϕ 1 , ϕ 2 , · · · , ϕ M ) (where ϕ j = Aφ j / λj , j = 1, 2, · · · , M). Step 4. Solve/compute the PODROEFD model Solve the PODROEFD scheme (1.2.17)–(1.2.18) to obtain the reduced-order so∗n ∗n T lution vectors u∗n = (u∗n 1 , u2 , · · · , um ) ; further, obtain the component forms ∗n ∗n uj,k = ui (0  j  J , 0  k  K, i = k(J + 1) + j + 1, 1  i  m = (K + 1)(J + 1)). Step 5. Check accuracy and renew POD basis to continue  Set δ = tB2,2 /x 2 + tC2,2 /y 2 . If (1 + δ)n−L λu(M+1)  μ. Then ∗n ∗n T u∗n = (u∗n 1 , u2 , · · · , um ) is just the solution vectors for the PODROEFD scheme (1.2.17)–(1.2.18) that satisfy the accuracy requirement. Else, i.e., if (1 + δ)n−L λu(M+1) > μ, put u1 = u∗(n−L) , u2 = u∗(n−L+1) , . . . , uL = u∗(n−1) and return to Step 2.

Reduced-Order Extrapolation Finite Difference Schemes Chapter | 1

19

FIGURE 1.2.1 The solution u obtained via the classical FD scheme.

1.2.6 A Numerical Example for the 2D Parabolic Equation A numerical example here is presented to demonstrate the advantage of the PODROEFD scheme (1.2.17)–(1.2.18). To this end, we take f (x, y, t) = 0, u(x, y, 0) = s(x, y) = sin πx sin πy, and u(x, y, t) = 0 in Eq. (1.2.1). The computational domain is taken as = {(x, y); 0  x  2, 0  y  2}. The spatial steps are taken as x = y = 0.02 and the time step t = 0.00005. Thus, we have 8t/x 2  1 and 8t/y 2  1 such that δ < 1. We first find the classical FD solution at t = 2 by means of the classical FD scheme (1.2.2), which is depicted graphically in Fig. 1.2.1. Next, we take 20 classical FD solutions of the classical FD scheme (1.2.2) √ in the first 20 steps as snapshots. By computing directly, we have achieved λ7  4 × 10−4 . It can be shown that as long as we take the first six POD bases, the theoretical accuracy requirement can be satisfied. Next, we find the reduced-order FD solution at t = 2 (n = 40 000) by means of the PODROEFD scheme (1.2.17)–(1.2.18), where it is unnecessary to update the POD basis and the solution is depicted graphically in Fig. 1.2.2. Because the PODROEFD scheme (1.2.17)–(1.2.18) only includes six unknown quantities and uses six optimizing data of the first 20 classical solutions as initial values, it saves computing time, reduces the degrees of freedom in the numerical computations, and alleviates the TE accumulation and, therefore, the PODROEFD solution is more efficient. Fig. 1.2.3 is the (log 10) error chart between the classical FD solutions and the reduced-order FD solution with different number of up to 20 POD bases at t = 2. It is shown that as long as we take M = 6, the reduced-order FD solutions obtained satisfy the accuracy requirement (i.e., its error does not exceed 4 × 10−4 ). Thus, the theoretical results are consistent with the ones from numerical calculations (they are all O(10−4 )). Hence, it is shown that the PODROEFD

20 Proper Orthogonal Decomposition Methods for Partial Differential Equations

FIGURE 1.2.2 The solution u obtained via the PODROEFD scheme.

FIGURE 1.2.3 The (log10) error plot between the classical FD solutions and the reduced-order FD solution with different number of up to 20 POD bases at t = 2.

scheme (1.2.17)–(1.2.18) is effective for the 2D parabolic equation. See the advantages and benefits of POD in Foreword and Introduction of the book.

1.3 A POD-BASED REDUCED-ORDER EXTRAPOLATION FINITE DIFFERENCE SCHEME FOR THE 2D NONSTATIONARY STOKES EQUATION In this section, a PODROEFD scheme is given for the 2D nonstationary Stokes equation. The PODROEFD scheme to produce the solutions on the time span [T0 , T ] (T0 T ) is obtained from the information on a short time span [0, T0 ] by extrapolation and iteration. The guidelines to choose the number of POD bases and renew the POD bases are provided, and an implementation for solving the PODROEFD scheme is given. Some numerical experiments are provided to illustrate the feasibility and efficiency of the PODROEFD scheme for simulating a channel flow with local expansion. The PODROEFD scheme for the 2D nonstationary Stokes equation is based on the work in [117].

Reduced-Order Extrapolation Finite Difference Schemes Chapter | 1

21

FIGURE 1.3.1 The domain of Problem 1.3.1.

1.3.1 Background for the 2D Nonstationary Stokes Equation The channel flow with local expansion in this section is motivated by applications, for example, for the case of the capillary blood vessel flow in the human body, where there is the microchannel flow with expansion geometries. It can be simplified into a nonstationary Stokes channel flow with local expansion [16,34, 144]. Its geometry is approximately made by two squares at the top and bottom of the channel. The flow domain is shown in Fig. 1.3.1. Thus, the mathematical model for the channel flow with local expansion is described by the following system of the nonstationary Stokes equation. Problem 1.3.1. Find (u, v) and p such that, for T > 0, ⎧   ⎪ 1 ∂ 2u ∂ 2u ∂p ∂u ⎪ ⎪ + − = f, (x, y, t) ∈ × (0, T ], + ⎪ ⎪ 2 2 ⎪ ∂t Re ∂x ∂x ∂y ⎪ ⎪   ⎪ ⎪ 1 ∂ 2v ∂ 2v ∂p ⎪ ∂v ⎪ ⎨ + − = g, (x, y, t) ∈ × (0, T ], + ∂t Re ∂x 2 ∂y 2 ∂y ⎪ ∂u ∂v ⎪ ⎪ + = 0, (x, y, t) ∈ × (0, T ], ⎪ ⎪ ⎪ ∂x ∂y ⎪ ⎪ ⎪ ⎪ u(x, y, t) = ϕ1 (x, y, t), v(x, y, t) = ϕ2 (x, y, t), (x, y, t) ∈ ∂ × (0, T ], ⎪ ⎪ ⎩ (x, y) ∈ , u(x, y, 0) = u0 (x, y), v(x, y, 0) = v0 (x, y), (1.3.1) where (u, v) is the velocity vector, p is the pressure, T is the total time duration, Re is the Reynolds number, f (x, y, t) and g(x, y, t) are the given body forces in the x-direction and y-direction, respectively, and ϕ1 (x, y, t), ϕ2 (x, y, t), u0 (x, y), and v0 (x, y) are four given functions. The nonstationary Stokes equation constitutes an important system of equations in fluid dynamics with many applications beyond the channel flow with local expansion [8,16,18,34,42,144,146,183]. For Problem 1.3.1, generally there are no analytical solutions. One has to rely on numerical solutions. The classical FD scheme with second-order accuracy in time and spatial variable discretizations is one of the simplest and most convenient high accuracy methods for

22 Proper Orthogonal Decomposition Methods for Partial Differential Equations

finding numerical solutions of the nonstationary Stokes equation. However, the approach involves a large number of degrees of freedom (unknown quantities). Especially, due to the TE accumulation in the computational process, it may also appear not to converge after several computation steps. Thus, an important task is to lessen the degrees of freedom so as to reduce the computational load and save CPU time in the process in a way that guarantees sufficiently accurate numerical solutions. Here, we employ the POD method to reduce its degrees of freedom.

1.3.2 A Classical Finite Difference Scheme for the 2D Nonstationary Stokes Equation and the Generation of Snapshots Let N be a positive integer, t = T /N the time step increment, x and y the spatial step increments in the x- and y-directions, respectively, and n denote the value of function u, v, and p at the points un 1 , v n 1 , and pj,k j + 2 ,k

j,k+ 2

(xj + 1 , yk , tn ), (xj , yk+ 1 , tn ), and (xj , yk , tn ) (0  j  J, 0  k  K, 0  n  2 2 N), respectively. Then Problem 1.3.1 is known to have the following classical FD scheme with second-order accuracy: 

uj + 1 ,k − uj − 1 ,k 2

2

x un+11

j + 2 ,k

2

2

2

2

y

n+1 = 0,

2t n n − pj,k ) + 2tfjn+ 1 ,k , (p x j +1,k 2 2t n n n − − pj,k ) + 2tgj,k+ (p 1, x j,k+1 2

= Fjn+ 1 ,k −

v n+1 1 = Gnj,k+ 1 j,k+ 2

+

vj,k+ 1 − vj,k− 1

(1.3.2) (1.3.3) (1.3.4)

where Fjn+ 1 ,k = un−11 j + 2 ,k 2  n 2t uj + 12 ,k−1 − 2uj + 12 ,k + uj + 12 ,k+1 uj − 12 ,k − 2uj + 12 ,k + uj + 32 ,k + + , Re y 2 x 2 Gnj,k+ 1 = v n−1 1 j,k+ 2 2   vj,k− 1 − 2vj,k+ 1 + vj,k+ 3 n 2t vj −1,k+ 12 − 2vj,k+ 12 + vj +1,k+ 12 2 2 2 + + . Re x 2 y 2 Inserting (1.3.3) and (1.3.4) into (1.3.2), one can obtain the approximate FD scheme of Poisson equations for p as follows:   pj −1,k − 2pj,k + pj +1,k pj,k−1 − 2pj,k + pj,k+1 n + = R, (1.3.5) x 2 y 2

Reduced-Order Extrapolation Finite Difference Schemes Chapter | 1 1 where R = tx [2t (fj + 1 ,k − fj − 1 ,k ) + Fj + 1 ,k − Fj − 1 ,k ]n + 2 2 2 2 [2t (gj,k+ 1 − gj,k− 1 ) + Gj,k+ 1 − Gj,k− 1 ]n . 2

2

2

23

1 ty

×

2

If tRe  4 and 4t  min{Rex 2 , Rey 2 } or t = o(Rex 2 , Rey 2 ), by taking the same approach as in the proof of Theorem 1.2.1, we easily prove that the FD scheme (1.3.3)–(1.3.5) is stable (see also [34,192]). If the solution triple (u, v, p) to Problem 3.2.1 is sufficiently regular, we have the following error estimates: n n (u(xj + 1 , yk , tn ), v(xj , yk+ 1 , tn ), p(xj , yk , tn )) − (unj+ 1 ,k , vj,k+ 1 , pj,k )2 2

2

2

2

= O(t , x , y ), n = 1, 2, · · · , N, 2

2

2

(1.3.6)

where  · 2 denotes the l 2 norm of the vector. If the Reynolds number Re, body force f (x, y, t) in the x-direction and body force g(x, y, t) in the y-direction, boundary value functions ϕ1 (x, y, t) and ϕ2 (x, y, t), initial value functions u0 (x, y) and v0 (x, y), time step increment t, and spatial step increments x and y are given, let u0 1 = u1

j + 12 ,k

= u0 (xj + 1 , yk ), v 0 2

j,k+ 12

= v1

j,k+ 12

j + 2 ,k

= v0 (xj , yk+ 1 ). By solving the FD 2

scheme (1.3.3)–(1.3.5), we can obtain the classical FD solutions un

j + 12 ,k

, vn

j,k+ 12

,

and (0  j  J, 0  k  K, 1  n  N ). n (i = kJ + j + 1, 1  i  m, Set uni = un 1 , vin = v n 1 , and pin = pj,k n pj,k

j + 2 ,k

j,k+ 2

m = (J + 1)(K + 1), 0  j  J − 1, 0  k  K − 1), respectively. We may choose the first L solutions from the set {uni , vin , pin }N n=1 (1  i  m) including (1  i  m, L N ) containing N ×m elements to construct a set {uli , vil , pil }L l=1 L × m elements, which are just snapshots. Remark 1.3.1. The snapshots are drawn from the first L FD solutions here. However, when one computes the real-world problems, one may also choose the ensemble of snapshots from physical system trajectories by drawing samples from the experiments.

1.3.3 Formulations of the POD Basis and the POD-Based Reduced-Order Extrapolating Finite Difference Scheme In the following, we first formulate a set of POD bases, and then establish the PODROEFD scheme with second-order accuracy for the nonstationary Stokes equation. The set of snapshots {uli , vil , pil }L l=1 (1  i  m) in Section 1.3.2 can be expressed in the following three m × L matrices As = (sil )m×L (s = u, v, p). Let λs1  λs1  · · ·  λs M˜ s > 0 (M˜ s = rank(As )) be the positive eigenvalues of the matrices As ATs (s = u, v, p) and let the matrices U s = (φ s1 , φ s2 , · · · , φ sm ) and the matrices V s = (ϕ s1 , ϕ s2 , · · · , ϕ sL ) consist of the

24 Proper Orthogonal Decomposition Methods for Partial Differential Equations

orthonormal eigenvectors of the matrices As ATs and ATs As , respectively. Then, it follows easily from linear algebra that As = U s D M˜ s V Ts (where D M˜ s = √ √ √ diag{ λs1 , λs2 , · · · , λs M˜ s , 0, · · · , 0}). Thus, U s = (φ s1 , φ s2 , · · · , φ sm ) (s = u, v, p) make up three sets of POD bases. The number of mesh points is far larger than that of the snapshots drawn, that is, m  L; the order m for matrices As ATs is far larger than the order L for matrices ATs As . But the numbers of their positive eigenvalues are identical, so we may first solve the eigenequation corresponding to matrices ATs As to find the eigenvectors ϕ sj (j = 1, 2, . . . , M˜ s ), and then, by the relationship  φ sj = As ϕ sj / λsj (j = 1, 2, . . . , M˜ s , s = u, v, p), we obtain the eigenvectors φ sj (j = 1, 2, . . . , M˜ s ) corresponding to the nonzero eigenvalues for matrix As ATs . Thus, by using the same technique as in Section 1.2.2, three optimizing orthonormal POD bases s = (φ s1 , φ s2 , · · · , φ sMs ) (Ms L, s = u, v, p) are formed by the first Ms (0 < Ms  M˜ s  L) columns from three groups of POD bases U s = (φ s1 , φ s2 , · · · , φ sm ). Further, by the properties of the norm of a matrix (see (1.2.12) in Section 1.2.2), the following error estimates hold:  As − AMs 2,2 = As − s Ts As 2,2 = λs(Ms +1) , s = u, v, p, (1.3.7) where A2,2 = supx Ax2 /x  · 2 is √ the l 2 vector norm, AMs = √ 2, √ T U s D Ms V s , and D Ms = diag{ λs1 , λs2 , · · · , λsMs , 0, · · · , 0} (0 < Ms  M˜ s  L). Let ε l (l = 1, 2, . . . , L) denote unit column vectors whose lth component is 1. Set s nm = (s1n , s2n , · · · , unm )T , s = u, v, p, n = 1, 2, · · · , N.

(1.3.8)

Then, from (1.3.7), the following error estimates hold: s lm − s Ts s lm 2 = (As − s Ts As )ε l 2  As − s Ts As 2,2 εl 2  = λs(Ms +1) , s = u, v, p, l = 1, 2, · · · , L,

(1.3.9)

T l l which  shows that s s s m is the optimal approximation for s m and the errors are λs(Ms +1) (s = u, v, p). By using the notation of (1.3.8), the classical FD scheme (1.3.3)–(1.3.5) is now written in the following vector scheme:

pnm = F˜1 (unm , v nm ), 0  n  N, (1.3.10)  n+1 n+1 T  n−1 n−1 T um , v m = um , v m + F˜ (unm , v nm , pnm ), 1  n  N − 1, where F˜ and F˜1 are determined from (1.3.3)–(1.3.4) and (1.3.5), respectively.

Reduced-Order Extrapolation Finite Difference Schemes Chapter | 1

25

Set !T  ∗n ∗n ∗n T um , v m , p m = u α nMu , v β nMv , p γ nMp ,

(1.3.11)

∗n ∗n ∗n ∗n ∗n ∗n T ∗n ∗n T ∗n where u∗n m = (u1 , u2 , · · · , um ) , v m = (v1 , v2 , · · · , vm ) , p m = (p1 , ∗n ∗n T p2 , · · · , pm ) are three column vectors corresponding to u, v, and p, respectively, and α nMu = (α1 , α2 , · · · , αMu )T , β nMv = (β1 , β2 , · · · , βMv )T , and γ nMp = (γ1 , γ2 , · · · , γMp )T . If unm , v nm , p nm in (1.3.10) are approximately re∗n ∗n placed with u∗n m , v m , p m in (1.3.11) (n = 0, 1, 2, · · · , N ), by noting that three matrices u , v , and p are formed with the orthonormal eigenvectors, the PODROEFD scheme for the channel flow with local expansion, with Mu + Mv + Mp (Mu , Mv , Mp L m) unknowns, is denoted by ⎧ n n n T n T n T n ⎪ ⎪ ⎪ α Mu = u um , β Mv = v v m , γ Mp = p pm , n = 1, 2, · · · , L, ⎪ ⎪ ⎪ ⎨ γ n = T F˜1 (u α n−1 , v β n−1 ), n = L, L + 1, · · · , N, p Mp Mu Mv (1.3.12) !T !T ⎪ n+1 n+1 n−1 n−1 ⎪ n n n ˜ ⎪ αM , β M = α Mu , β Mv + G(α Mu , β Mv , γ Mp ), ⎪ u v ⎪ ⎪ ⎩ n = L, L + 1, · · · , N − 1,

˜ n , β n , γ n ) = (Tu , Tv )T F˜ (u α n , v β n , p γ n ). where G(α Mv Mv Mu Mp Mu Mp After α nMu , β nMv , and γ nMp are obtained from (1.3.12), the reduced-order FD solutions for the PODROEFD scheme are obtained by n u∗n m = u α Mu ,

n v ∗n m = v β Mv ,

Thus, in component forms: u∗n

j + 12 ,k

p ∗n m = p γ Mp , n = 0, 1, 2, · · · , N. (1.3.13) ∗n = u∗n i ,v

j,k+ 12

= vi∗n , and p ∗n

j,k+ 12

= pi∗n (0 

j  J , 0  k  K, i = k(J + 1) + j + 1, 1  i  m = (K + 1)(J + 1)); the PODROEFD scheme with second-order accuracy is attained. Remark 1.3.2. Since the classical FD scheme (1.3.3)–(1.3.5) at each time level includes 3m degrees of freedom, while the system of Eqs. (1.3.12) and (1.3.13) at each time level (when n > L) contains only (Mu + Mv + Mp ) degrees of freedom (Mu , Mv , Mp L m, for example, in Section 1.3.6, L = 20, Mu = Mv = Mp = 6, but m = 136 × 104 ), the system of Eqs. (1.3.12) and (1.3.13) is the PODROEFD scheme with much fewer degrees of freedom and the second-order accuracy and has no repetitive computations. Thus, it has the same advantage and efficiency as the PODROEFD scheme, just as in Section 1.2.

1.3.4 Error Estimates and a Criterion for Renewing the POD Basis In the following, we estimate the errors between the reduced-order FD solutions for (1.3.12) and (1.3.13) and the classical FD solutions for (1.3.3)–(1.3.5).

26 Proper Orthogonal Decomposition Methods for Partial Differential Equations

By using (1.3.13), we write the second and third equations of (1.3.12) in the following vector form: ˜ ∗n−1 , v ∗n−1 p∗n ), L + 1  n  N, m = F1 (um m  ∗n+1 ∗n+1 T  ∗n−1 ∗n−1 T um , v m = um , v m ∗n ∗n +F˜ (u∗n m , v m , p m ), L  n  N − 1,

(1.3.14)

where their stability conditions are also made to satisfy tRe  4 and 4t  min{Rex 2 , Rey 2 } (see [34] or [192]). Set en = (unm , v nm , p nm )T − ∗n ∗n T (u∗n m , v m , p m ) . By the first equation of (1.3.12) and (1.3.13), we obtain en = (un , v n , pn )T − (u Tu un , v Tv v n , p Tp p n )T ,

(1.3.15)

where n = 1, 2, · · · , L. From (1.3.9) and (1.3.15), we obtain en 2 



λu(Mu +1) +



λv(Mv +1) +

" λp(Mp +1) , n = 1, 2, · · · , L. (1.3.16)

By (1.3.10) and (1.3.14), we have en+1 2  en−1 2 + Men 2 , n = L, L + 1, · · · , N − 1,

(1.3.17)

where M = max{tRe/4, 4t/(Re x 2 ), 4t/(Re y 2 )}  1 under the stability conditions tRe  4 and 4t  min{Rex 2 , Rey 2 }. If t = o(x, y, Re x 2 , Re y 2 ), M is a small positive constant. Summing (1.3.17) from L to n − 1 yields en 2  eL−1 2 + eL 2 + M

n−1 

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

i=L

(1.3.18)

n−1

If we write ξn = M i=L ei 2 + eL−1 2 + eL 2 and δ = M + 1, then we have en 2  ξn and ξn − ξn−1 = Men−1 2 (n  2). Therefore, we get ξn  (M + 1)ξn−1 = δξn−1  δ 2 ξn−2  · · ·  δ n−L ξL = δ n−L (eL−1 2 + eL 2 ).

(1.3.19)

Set C(δ n ) = 2δ n−L . We obtain from (1.3.16)

"   en 2  2δ n−L [ λu(Mu +1) + λv(Mv +1) + λp(Mp +1) ] "   = C(δ n )[ λu(Mu +1) + λv(Mv +1) + λp(Mp +1) ], (1.3.20)

where n = L + 1, L + 2, ..., N . Synthesizing the above discussions yields the following result.

Reduced-Order Extrapolation Finite Difference Schemes Chapter | 1

27

Theorem 1.3.1. If (unm , v nm , p nm )T (n = 1, 2, · · · , N ) are the solution vectors formed by the solutions for the classical FD scheme (1.3.3)–(1.3.5), ∗n ∗n T (u∗n m , v m , p m ) (n = 1, 2, · · · , N ) are the reduced-order FD solutions for the PODROEFD scheme (1.3.12) and (1.3.13), we have the following error estimates: ∗n ∗n (unm , v nm , p nm ) − (u∗n m , v m , p m )2 # $ "   C(δ n ) λu(Mu +1) + λv(Mv +1) + λp(Mp +1) , n = 1, 2, · · · , N,

where C(δ n ) = 1 (1  n  L) and C(δ n ) = 2(1 + M)n−L (L + 1  n  N ) together with M = max{tRe/4, 4t/(Re x 2 ), 4t/(Re y 2 )}. Since the absolute value of each vector component is less than its norm, combining (1.3.6) with Theorem 1.3.1 yields the following result. Theorem 1.3.2. The exact solution for the nonstationary Stokes equation and the reduced-order FD solutions obtained from the PODROEFD scheme (1.3.12) and (1.3.13) satisfy the following error estimates: for 1  n  N , | u(xj + 1 , yk , tn ) − u∗n

| + | v(xj , yk+ 1 , tn ) − v ∗n 1 | j + 12 ,k j,k+ 2 2 ∗n | + | p(xj , yk , tn ) − pj,k

     = O C(δ n ) λu(Mu +1) + λv(Mv +1) + λp(Mp +1) , t 2 , x 2 , y 2 . 2

Remark 1.3.3. The error estimates of Theorem 1.3.2 give a guide for choosing  the number  of POD bases,  namely, taking Mu , Mv , and Mp such that λu(Mu +1) + λv(Mv +1) + λp(Mp +1) = O(t 2 , x 2 , y 2 ). Here, C(δ n ) = 2(1 + M)n−L (L + 1  n  N ) are caused by extrapolating iteration  and they may act as a guide for renewing POD bases, namely, when C(δ n )[ λu(Mu +1) +   λv(Mv +1) + λp(Mp +1) ] > max(t 2 , x 2 , y 2 ), it is time for the renewal of POD bases.

1.3.5 Implementation for the POD-Based Reduced-Order Extrapolating Finite Difference Scheme The implementation for the PODROEFD scheme (1.3.12) and (1.3.13) consists of the following five steps. Step 1. Classical FD computation and formulation of snapshots Solving the classical FD scheme (1.3.3)–(1.3.5) at the first few L steps (empirically, say, take L = 20) yields the classical FD solutions un 1 , v n 1 , and j + 2 ,k

j,k+ 2

n (0  j  J, 0  k  K, 1  n  L) and further produces a set of snapshots pj,k n n n {uli , vil , pil }L l=1 (1  i  m) with L × m elements, where ui = u 1 , vi =

vn

j,k+ 12

j + 2 ,k

n (i = kJ + j + 1, 1  i  m, m = (J + 1)(K + 1), 0  , and pin = pj,k

j  J − 1, 0  k  K − 1), respectively.

28 Proper Orthogonal Decomposition Methods for Partial Differential Equations

Step 2. Snapshot matrices As and eigenvalues of ATs As n )T (s = u, v, p, n = 1, 2, · · · , N ). Formulate the snapLet s nm = (s1n , s2n , · · · , sm shot matrices As = (sil )m×L (s = u, v, p) and find the eigenvalues λs1  λs2  · · ·  λs M˜ s > 0 (M˜ s = rank As ) and corresponding eigenvectors ϕ sj (j = 1, 2, · · · , M˜ s , s = u, v, p) of ATs As (s = u, v, p).

Step 3. Choice of POD bases For the error tolerance μ = O(t 2 , x 2 , y 2 ) desired,  decide thenumbers Ms (Ms  M˜ s , s = u, v, p) of POD bases such that λu(Mu +1) + λv(Mv +1) +  λp(Mp +1)  μ, and formulate the POD bases s = (φ s1 , φ s2 , · · · , φ sMs )  (where φ sj = As ϕ sj / λsj , j = 1, 2, · · · , Ms , s = u, v, p). Step 4. Solve/compute the PODROEFD model Solving the PODROEFD scheme (1.3.12) and (1.3.13) yields the reduced∗n ∗n ∗n ∗n ∗n ∗n ∗n order solution vectors u∗n m = (u1 , u2 , · · · , um ), v m = (v1 , v2 , · · · , vm ), ∗n ∗n ∗n ∗n and p m = (p1 , p2 , · · · , pm ); further, the process produces the components ∗n ∗n ∗n ∗n (0  j  J , 0  k  K, u∗n 1 = u∗n 1 = vi , and p 1 = pi i , v j + 2 ,k

j,k+ 2

j,k+ 2

i = k(J + 1) + j + 1, 1  i  m = (K + 1)(J + 1)). Step 5. Check accuracy and renew POD bases to continue Set M = max{0.25tRe, 4t/(Re x 2 ), 4t/(Re y 2 )}. If 2(1 + M)n−L

#

λu(Mu +1) +



λv(Mv +1) +

$ " λp(Mp +1)  μ,

∗n ∗n ∗n ∗n ∗n ∗n ∗n ∗n ∗n then u∗n m = (u1 , u2 , · · · , um ), v m = (v1 , v2 , · · · , vm ), and p m = (p1 , ∗n ∗n p2 , · · · , pm ) (n = 1, 2, · · · , N ) are exactly the solutions satisfying the desirable accuracy. Else, namely, if

2(1 + M)n−L

#

λu(Mu +1) + ∗(n−l)

l ) = (s put (s1l , s2l , · · · , sm 1 p) and return to Step 2.

$ "  λv(Mv +1) + λp(Mp +1) > μ,

∗(n−l)

, s2

∗(n−l)

, · · · , sm

) (l = 1, 2, · · · , L, s = u, v,

∗n Remark 1.3.4. Step 5 could be adapted as follows: if u∗n−1 − u∗n m m 0  um − ∗n+1 ∗n−1 ∗n ∗n ∗n+1 ∗n−1 ∗n ∗n um 2 , v m − v m 0  v m − v m 2 , and pm − pm 2  pm − ∗n ∗n p∗n+1 2 (n = L, L + 1, · · · , N − 1), then (u∗n m m , v m , p m ) (n = 1, 2, · · · , N ) are the reduced-order solution vectors for PODROEFD scheme (1.3.12) and (1.3.13) satisfying the desirable accuracy. Else, namely, if u∗n−1 − u∗n m m 0 < ∗n ∗n+1 ∗n−1 ∗n ∗n ∗n+1 ∗n−1 um − um 2 , v m − v m 0 < v m − v m 2 , or pm − p ∗n m 2 < ∗(n−i) ∗n ∗n+1 i pm − pm 2 (n = L, L + 1, · · · , N − 1), let s m = s m (i = 1, 2, · · · , L, s = u, v, p) and return to Step 2.

Reduced-Order Extrapolation Finite Difference Schemes Chapter | 1

29

FIGURE 1.3.2 When Re = 1000, the top chart (A) and the bottom chart (B) are the contours of the classical FD solution and the reduced-order FD solution of the velocity (u, v) at time instant t = 9, respectively.

1.3.6 Some Numerical Experiments for the 2D Nonstationary Stokes Equation In the following, we present some numerical experiments of the channel flow with local expansions, i.e., with two squared protrusions at the middle top and bottom of the channel to validate the feasibility and efficiency of the PODROEFD scheme with second-order accuracy and to show that the numerical results are consistent with theoretical estimates. ¯ is given as in Fig. 1.3.1. Take We assume that the computational domain Re = 1000, f = g = 0. Except for the inflow on the left boundary with the periodic flow velocity of (u, v) = (0.1(y − 2)(8 − y) sin 2πt, 0) (2  y  8) and outflow on the right boundary with velocity of (u, v) satisfying v = 0 and ∂u/∂x = 0, all initial and boundary value conditions are taken as 0. We divide ¯ into mesh by taking spatial step increments as x = y = 10−2 and then the taking time step increment as t = 0.001. We have obtained a numerical solution (un 1 , v n 1 ) of velocity (u, v) and a numerical solution p n

j,k+ 12

j + 2 ,k

j,k+ 2

of pressure p by the classical FD schemes

(1.3.3)–(1.3.5) with n = 9000 (i.e., t = 9), which are depicted graphically as the contours in Figs. 1.3.2(A) and 1.3.3(A), respectively. We have only employed the first L = 20 numerical solutions (un 1 , v n 1 , pn

j,k+ 12

j + 2 ,k

j,k+ 2

) (n = 1, 2, · · · , 20) from the classical FD scheme, forming 20 snap-

shot vectors (unm , v nm , p nm ) (n = 1, 2, · · · , 20, m = 136 × 104 ). Afterwards, by Step 2 in Section 1.3.5, we have found three groups of 20 eigenvalues λuj , λvj , and λpj (j = 1, 2, · · · , 20) which are arranged in a nonincreasing order and three groups of 20 eigenvectors ϕ uj , ϕ vj , and ϕ pj (j = 1, 2, · · · , 20) corresponding to the 20 eigenvalues. By computing, we have achieved that

30 Proper Orthogonal Decomposition Methods for Partial Differential Equations

FIGURE 1.3.3 When Re = 1000, the top chart (A) and the bottom chart (B) are the classical FD solution and the reduced-order FD solution of the pressure field p at the time instant t = 9, respectively.

 √ √ the eigenvalues satisfy λu7 + λv7 + λp7  4 × 10−4 , which guides us to take the first three groups of six eigenvectors as POD bases by Step 3 in Section 1.3.5, i.e.,the POD bases are taken as u = (φ u1 , φ u2 , · · · , φ u6 ) (with φ uj = Au ϕ uj / λuj , j = 1, 2, · · · , 6), v = (φ v1 , φ v2 , · · · , φ v6 ) (with  φ vj = Av ϕ vj / λvj , j = 1, 2, · · · , 6), and p = (φ p1 , φ p2 , · · · , φ p6 ) (with  φ pj = Ap ϕ pj / λpj , j = 1, 2, · · · , 6). Finally, when n = 9000 (i.e., at t = 9) and Re = 103 , the numerical solutions u∗n 1 , v ∗n 1 , and p ∗n 1 are obtained j + 2 ,k

j,k+ 2

j,k+ 2

from the PODROEFD scheme by Step 4 in Section 1.3.5. It is necessary to renew POD bases once t = 5. The reduced-order FD solution at t = 9 is obtained and depicted graphically in Figs. 1.3.2(B) and 1.3.3(B), respectively. Each pair of figures (A) and (B) in Figs. 1.3.2(A) and 1.3.3(B) has exhibited close similarity, respectively, but the reduced-order FD solutions obtained from the PODROEFD scheme are computed with higher efficiency than the classical FD solutions since the PODROEFD scheme needs much fewer degrees of freedom and, thus, it can also significantly reduce the TE accumulation in the computational process. Fig. 1.3.4 shows the mean absolute errors (MAEs) between the reducedorder FD solutions obtained from the PODROEFD scheme with second-order accuracy and a different number of POD bases and the solutions obtained from classical FD schemes (1.3.3)–(1.3.5) when Re = 1000 at t = 9. Comparing the classical FD scheme with the PODROEFD scheme containing six optimal bases and implementing the numerical simulation computations with Re = 1000 at t = 9, we have found that, for the classical FD schemes (1.3.3)–(1.3.5) containing 3 × 136 × 104 unknown quantities at each time level, the required computing time was about 48 minutes on a laptop, while, for the PODROEFD scheme with

Reduced-Order Extrapolation Finite Difference Schemes Chapter | 1

31

FIGURE 1.3.4 The MAEs between the reduced-order FD solutions with different number of POD bases and the classical FD solution when Re = 1000 at t = 9.

six optimal bases only including 3 × 6 unknown quantities at the same time level, the corresponding time was only 16 seconds on the same laptop, so the ratio of speed is 180:1, while the errors between their solutions do not exceed 4 × 10−4 . Though our experiments are in a sense recomputing what we have already computed by the classical FD scheme at the first few L = 20 steps, when we compute actual problems, we may construct the snapshots and POD basis by drawing samples from experiments and then solve directly the PODROEFD scheme with second-order accuracy, while it is unnecessary to solve the classical FD schemes (1.3.3)–(1.3.5). Thus, the time consuming calculations and resource demands in the computational process have resulted in great saving. It has also shown that finding the approximate solutions for the nonstationary Stokes equation using the PODROEFD scheme with second-order accuracy is computationally effective. The numerical results are consistent with those obtained for the theoretical case as neither the theoretical nor the numerical errors exceeds 4 × 10−4 . In addition, if one uses the reduced-order FD scheme with first-order time accuracy as in [113] to find the reduced-order FD solution at t = 9, it is necessary to take the time step as k = 10−4 and carry out 9 × 104 steps in order to obtain the same accuracy as here. Thus, its computing load is 10 times that of the PODROEFD scheme with second-order accuracy, and its TE accumulation in the computational process increases significantly as well as it is repeating computations of the classical first-order time accuracy FD scheme on the same time interval [0, T ]. Therefore, the PODROEFD scheme with second-order accuracy here is different from the existing reduced-order schemes and it offers improvement over the existing reduced-order methods [2,5,11,19,38,40,45,48,51,52,57,58,62, 63,91,113,118,122,128,135,137,141,155,156,166,167,170,171,184–186,199].

32 Proper Orthogonal Decomposition Methods for Partial Differential Equations

1.4 POD-BASED REDUCED-ORDER EXTRAPOLATING FINITE DIFFERENCE SCHEME FOR 2D SHALLOW WATER EQUATION In this section, we employ the POD method to establish a PODROEFD scheme with very few degrees of freedom for the 2D shallow water equations (SWEs) with sediment concentration. We also provide the error estimates between the accurate solution and the classical FD solutions as well as those between the exact solution and the PODROEFD solutions. Moreover, we present two numerical simulations to illustrate that the PODROEFD scheme can greatly reduce the computational load. Thus, both the feasibility and efficiency of the PODROEFD scheme are validated. The main reference for the work here is [96].

1.4.1 Model Background and Survey for the 2D Shallow Water Equation A system of SWEs can be used to describe the propagation and evolution of short waves in shallow waters, also referred to as the Saint-Venant system (see [36]). It has extensive applications in ocean, environmental, and hydraulic engineering. In particular, in coastal engineering, [46] discusses applications to various problems in open-channel flows in rivers and reservoirs, tidal flows in estuary and coastal water regions, bore wave propagation, and the stationary hydraulic jump in river. Because SWEs are a system of nonlinear PDEs, they generally have no analytical solutions. One has to rely on numerical solutions. Here we mention a number of references on the study of numerical solutions for the 2D SWEs including only the continuity equation and the momentum equation, modeling the effects of the water depth and the velocity of fluid; for example, the finite volume (FV) method on unstructured triangular meshes in Anatasiou and Chan [9], the upwind methods in Bermudez and Vazquez [17], the parallel block preconditioning techniques in Cai and Navon [22], the optimal control technique of finite element (FE) limited-area in Chen and Navon [30], the least-squares FE method in Liang and Hsu [74], the FD Lax–Wendroff weighted essentially nonoscillatory (WENO) schemes in Lu and Qiu [77], the FE simulation technique in Navon [131], the FD WENO schemes in Qiu and Shu [136], the Roe approximate Riemann solver technique in Rogers et al. [142], the essentially nonoscillatory and WENO schemes with the exact conservation property in Vukovic and Sopta [172], the explicit multiconservation FD scheme in Wang [173], the composite FV method on unstructured meshes in Wang and Liu [174], the high-order FD WENO schemes in Xing and Shu [180], the high-order well-balanced FV WENO schemes and discontinuous Galerkin (DG) methods in Xing and Shu [181], the positivity preserving high-order wellbalanced DG methods in Xing et al. [182], the dispersion–correction FD scheme in Yoon et al. [189], the nonoscillatory FV method in Yuan and Song [190], the

Reduced-Order Extrapolation Finite Difference Schemes Chapter | 1

33

surface gradient method in Zhou et al. [196], and the total variation diminishing FD scheme in Wang et al. [175]. Nevertheless, the transport and sedimentation of silt and sand are some important processes in causing change to the natural environment such as the formation and evolution of a delta, expansion of alluvial plains, detouring of rivers, etc. They also cause some serious problems that should be carefully addressed to in hydraulic problems, such as irrigation systems, transportation channels, hydroelectric stations, ports, and other coastal engineering works. A model for the 2D SWEs including sediment concentration is available in [191] with some numerical methods based on an optimal control approach (see [198]) and a mixed FE technique (see [125, 126]). It is well known that the model based on a classical FD scheme in [191] is one of the simplest and most convenient methods for solving 2D SWEs with sediment concentration. However, it also contains many degrees of freedom, i.e., unknown quantities. Therefore, its computational complexity is also higher with previously mentioned adverse effects in generating errors. It is advantageous to build the reduced-order FD scheme with sufficiently high accuracy and very few degrees of freedom. Here, we continue with the PODROEFD scheme for the 2D SWEs with sediment concentration. Some POD-based reduced-order models for the 2D SWEs already exist (see, e.g., [29,151,152,199]), but these POD-based reduced-order models for the 2D SWEs have not included the sediment concentration effect. They also employ the numerical solutions obtained from the classical numerical methods on the entire time span [0, T ] to formulate the POD basis, build the PODbased reduced-order models, and recompute the solutions on the same time span [0, T ]; these also belong to repeated computations. Here we thoroughly improve the existing methods, where we only adopt the first few snapshots of given classical FD numerical solutions for the 2D SWEs on a very short time span [0, T0 ] (T0 T ) to formulate the POD basis and build the PODROEFD scheme, before finding the numerical solutions on the total time span [0, T ] by extrapolation and iteration as well as the POD basis update. Thus, an important advantage of the POD method is kept, i.e., using the given data on a very short time span [0, T0 ] to predict future physics on the whole time span [T0 , T ]. So this study is useful as a motivating example for real-world computational problems and big data. In the following, we first devote ourselves to the formulation of the snapshots and the POD basis from the classical FD solutions for the 2D SWEs with sediment concentration and the PODROEFD scheme, and then we provide the error estimates of solutions and the implementation of the algorithm for the PODROEFD scheme. We will provide two numerical simulation examples to verify the reliability and effectiveness for our PODROEFD scheme.

34 Proper Orthogonal Decomposition Methods for Partial Differential Equations

1.4.2 The Governing Equations and the Classical FD Scheme for the 2D Shallow Water Equation Including Sediment Concentration Let ⊂ R 2 be a bounded and connected domain. The governing equations for the 2D SWEs with sediment concentration are known as follows (see [46,191], with notational adaptation): ∂Z ∂(Zu) ∂(Zv) + + ∂t ∂x ∂y   2 ∂ Z ∂ 2Z + 2 , (x, y, t) ∈ × (0, T ), (1.4.1) =γ ∂x 2 ∂y  2  ∂ u ∂ 2u ∂u u∂u v∂u + + + −fv =A ∂t ∂x ∂y ∂x 2 ∂y 2 √ ∂(Z + zb ) CD u u2 + v 2 − , (x, y, t) ∈ × (0, T ), (1.4.2) −g ∂x Z  2  ∂ v ∂ 2v ∂v u∂v v∂v + + + +fu=A ∂t ∂x ∂y ∂x 2 ∂y 2 √ ∂(Z + zb ) CD v u2 + v 2 − , (x, y, t) ∈ × (0, T ), (1.4.3) −g ∂y Z ∂S u∂S v∂S + + ∂t ∂x ∂y   2 αω(S − S ∗ ) ∂ S ∂ 2S + , (x, y, t) ∈ × (0, T ), (1.4.4) + =ε Z ∂x 2 ∂y 2   ∂zb αω(S − S ∗ ) ∂u ∂v + gb + = , (x, y, t) ∈ × (0, T ), (1.4.5) ∂t ∂x ∂y ρ where γ (m2 /s) and A (m2 /s) are two viscosity coefficients, (u, v) (m/s) is the velocity vector, Z = z − zb (m) the water depth, z (m) the surface height, zb (m) the height of the river bed (see Fig. 1.4.1), f (1/s) the Coriolis constant, g (m/s2 ) the gravitational constant, CD (nondimensional) the coefficient of bottom drag, ε (m2 /s) the diffusion coefficient of sand, ω (m/s) the falling speed of suspended sediment particles, S (kg/m3 ) the concentration of sediment in water, ρ (kg/m3 ) the density of dry sand (taken as a constant), α (nondimensional) the constant of sediment variety, S ∗ = K[(u2 + v 2 )3/2 /(gωZ)]l the capability of sediment transport in a bottom bed (a given empirical function), gb = (u2 + v 2 )3/2 Z p d q [1 − vc /(u2 + v 2 )1/2 ] also a given empirical function, vc (m/s) the velocity of sediment mass transport (a given function, too), d (m) the diameter of sediment, and K (kg/m3 ), l (nondimension),  (s3 /m2 ), p (nondimensional), and q = −p are all empirical constants. The boundary conditions are assumed as follows: Z(x, y, t) = Z0 (x, y, t), u(x, y, t) = u0 (x, y, t), v(x, y, t) = v0 (x, y, t),

Reduced-Order Extrapolation Finite Difference Schemes Chapter | 1

35

FIGURE 1.4.1 Water profile.

S(x, y, t) = S0 (x, y, t), zb (x, y, t) = zb0 (x, y, t), (x, y, t) ∈ ∂ × (0, T ),

(1.4.6)

where Z0 (x, y, t), u0 (x, y, t), v0 (x, y, t), S0 (x, y, t), and zb0 (x, y, t) are all given functions. The initial conditions are assumed as follows: Z(x, y, 0) = Z 0 (x, y), u(x, y, 0) = u0 (x, y), v(x, y, 0) = v 0 (x, y), S(x, y, 0) = S 0 (x, y), zb (x, y, 0) = zb0 (x, y), x ∈ ∂ × (0, T ), (1.4.7) where Z 0 (x, y), u0 (x, y), v 0 (x, y), S 0 (x, y), and zb0 (x, y) also are all given functions. Let  t be the time step, let x and y be the spatial steps, and let N = T /t . By discretizing (1.4.1), (1.4.4), and (1.4.5) at a reference point (xj , yk , tn ), (1.4.2) at a reference point (xj + 1 , yk , tn ), and (1.4.3) at a reference 2 point (xj , yk+ 1 , tn ), we obtain the classical FD scheme for the 2D SWEs with 2 sediment concentration as follows:  n  n + Zn n n + Zn − 2Zj,k Zj +1,k − 2Zj,k Zj,k+1 j −1,k j,k−1 n+1 Zj,k = tγ + x 2 y 2 ⎛ n ⎞ u 1 Z n 1 − un 1 Z n 1 vn 1 Zn 1 − vn 1 Zn 1 j + 2 ,k j + 2 ,k j − 2 ,k j − 2 ,k j,k+ 2 j,k+ 2 j,k− 2 j,k− 2 ⎠ − t ⎝ + x y n , + Zj,k

tCD un 1 j + 2 ,k

(1.4.8)

"

(un 1 )2 + (v n 1 )2 j + 2 ,k j + 2 ,k n+1 u 1 =− n j + 2 ,k Z 1 j + 2 ,k ⎛ n n u 3 − 2u 1 + un 1 un 1 j + 2 ,k j + 2 ,k j − 2 ,k j + 2 ,k+1 + tA ⎝ + 2 x ⎛

− t ⎝

un

j + 12 ,k

(unj+1,k − unj,k ) x

vn +

j + 12 ,k

(un

j + 12 ,k+ 12

y

− 2un

+ un

− un

)

j + 12 ,k y 2

j + 12 ,k− 12

j + 12 ,k−1

⎞ ⎠

⎞ ⎠

36 Proper Orthogonal Decomposition Methods for Partial Differential Equations

− gt

n n n Zjn+1,k + zb,j +1,k − Zj,k − zb,j,k

x "

v n+1 1 = − + tA ⎝ ⎛ − t ⎝ − gt

vn

j +1,k+ 12

un

j,k+ 12

+

(v n

j + 12 ,k+ 12

− vn

j − 12 ,k+ 12

vn

)

y  n n + Sn Sj +1,k − 2Sj,k j −1,k x 2

tunj,k (S n 1 j + ,k 2

− Sn 1 ) j − 2 ,k

x n − S ∗n ) αωt (Sj,k j,k n Zj,k

n αωt(Sj,k

ρ

−

j,k− 12

⎠

⎞ ⎠

− t(f u)nj,k+ 1 ,

(1.4.10)

2

+

⎞

+ vn

n n ) (vj,k+1 − vj,k

y

n n + Sn − 2Sj,k Sj,k+1 j,k−1



y 2

n (S n tvj,k j,k+ 1

2

− Sn 1 ) j,k− 2

y

n + Sj,k ,

(1.4.11)

n (un tgb,j,k

∗n ) − Sj,k

j,k+ 12

−

n n n − zn + zb,j,k+1 − Zj,k Zj,k+1 b,j,k

n+1 n = zb,j,k − zb,j,k

+

2

(un 1 )2 + (v n 1 )2 j,k+ 2 j,k+ 2 n + vj,k+ 1 n Z 2 1 j,k+ 2 − 2v n 1 + v n v n 3 − 2v n 1 j,k+ 2 j −1,k+ 12 j,k+ 2 j,k+ 2 + 2 x y 2

x

n+1 = tε Sj,k

−

2

tCD v n 1 j,k+ 2

j,k+ 2

⎛

+ unj+ 1 ,k + t(f v)nj+ 1 ,k , (1.4.9)

j + 12 ,k

− un

j − 12 ,k

x

) −

n (v n tgb,j,k

j,k+ 12

− vn

j,k− 12

)

y (1.4.12)

,

where n = 1, 2, · · · , N , j = 1, 2, · · · , J , k = 1, 2, · · · , K, J = max{[| x1 − x2 |]: (x1 , y), (x2 , y) ∈ }, and K = max{[| y1 − y2 |] : (x, y1 ), (x, y2 ) ∈ }. In order to prove the stability of the FD scheme (1.4.8)–(1.4.12), it is necessary to introduce the following discrete Gronwall lemma (see [80]). Lemma 1.4.1 (The discrete Gronwall lemma). If {an } and {bn } are two nonnegative sequences, and {cn } is a positive monotone sequence, that satisfy an + bn  cn + λ¯

n−1 

ai (λ¯ > 0),

a 0 + b0  c 0 ,

i=0

then ¯ an + bn  cn exp(nλ),

n = 0, 1, 2, · · · .

For the classical FD scheme (1.4.8)–(1.4.12), we have the following result.

Reduced-Order Extrapolation Finite Difference Schemes Chapter | 1

37

Theorem 1.4.2. Under the conditions t · (| u | + | v |)  min{4γ , 4ε, 4A} and 4t max{γ , A, ε}  min{x 2 , y 2 }, the classical FD scheme (1.4.8)–(1.4.12) is locally stable. Further, we have the following error estimates: n | Z(xj , yk , tn ) − Zj,k | + | u(xj + 1 , yk , tn ) − unj+ 1 ,k | 2

2

n n + | v(xj , yk+ 1 , tn ) − vj,k+ 1 | + | S(xj , yk , tn ) − Sj,k | 2

2

n + | zb (xj , yk , tn ) − zb,j,k |= O(t, x 2 , y 2 ),

1  n  N, 1  j  J, 1  k  K.

(1.4.13)

Proof. If γ t/x 2  1/4, γ t/y 2  1/4, and 4t (| u | + | v |)  min{γ , ε, A}, implying 4t (u∞ + v∞ )  min{γ , ε, A}, by (1.4.8), we have   2γ t 2γ t n+1 n | Zj,k | 1 − − | | Zj,k x 2 y 2 γ t γ t n n + (| Zjn+1,k | + | Zjn−1,k |) + (| Zj,k+1 | + | Zj,k−1 |) x 2 y 2 t (| unj+ 1 ,k | · | Zjn+ 1 ,k | + | unj− 1 ,k | · | Zjn− 1 ,k |) + x 2 2 2 2 t n n n n | + | vj,k− |) (| vj,k+ + 1 |·|Z 1 |·|Z j,k+ 12 j,k− 12 y 2 2   2t 2t u∞ + v∞ Z n ∞  1+ x y   γ γ  1+ (1.4.14) + Z n ∞ , 2x 2y where  · ∞ is the L∞ ( ) norm. Thus, from (1.4.14), we obtain   γ γ Z n ∞  1 + + Z n−1 ∞ , n = 1, 2, · · · , N. 2x 2y

(1.4.15)

By summing (1.4.15) from 1 to n, we obtain  Z n ∞  Z 0 ∞ +

 n−1 γ γ Z j ∞ , n = 1, 2, · · · , N. (1.4.16) + 2x 2y j =0

By applying the discrete Gronwall lemma, Lemma 1.4.1, to (1.4.16), we obtain   nγ nγ Z n ∞  Z 0 ∞ exp + , n = 1, 2, · · · , N, (1.4.17) 2x 2y & % showing that the series Z n+1 is locally stable when the time interval [0, T ] is finite. Further, it is convergent from the stability theories of FD schemes (see

38 Proper Orthogonal Decomposition Methods for Partial Differential Equations

[34] or [76]). As the water depth is positive, there are two positive constants β1 and β2 such that β1  Z n ∞  β2 , n = 0, 1, 2, · · · , N.

(1.4.18)

If 4At  min{x 2 , y 2 } and 4εt  min{x 2 , y 2 }, by using the same technique as when proving (1.4.15), from (1.4.9)–(1.4.12) and (1.4.18), we obtain   A 2gt n A un+1 ∞  1 + + un ∞ + Z ∞ 2x 2y x 2gt n CD A + (un ∞ + v n ∞ ), (1.4.19) zb ∞ + tf ∞ v n ∞ + x 4β1   A 2gt n A + v n ∞ + Z ∞ v n+1 ∞  1 + 2x 2y y 2gt n CD A + (un ∞ + v n ∞ ), (1.4.20) zb ∞ + tf ∞ un ∞ + y 4β1   ε ε n+1 + S n ∞ S ∞  1 + 2x 2y αωt + (S n ∞ + S ∗n ∞ ), (1.4.21) β1   n u ∞ v n ∞ + zbn+1 ∞  2tgbn ∞ x y αωt + (1.4.22) (S n ∞ + S ∗n ∞ ), ρ where n = 0, 1, 2, · · · , N − 1. Note that S ∗n ∞  K[(un 2∞ + v n 2∞ )3/2 /(gωβ1 )]m  K[(A/t)2m /(gωβ1 )m ](un ∞ + v n ∞ ) and p

p

gbn ∞  β2 d q (un 2∞ + v n 2∞ )3/2  β2 d q (A/t)3 . Set  = max{K[(A/t)2m /(gωβ1 )m ]αωt/(β1 + ρ) + A/(2x + 2y) p

+ tf ∞ + 2CD A/(4β1 ) + 2β2 d q A3 /(xt 2 ), 2gt/(x + y) + K[(A/t)2m /(gωβ1 )m ]αωt/(β1 + ρ), p

A/(2x + tf ∞ + 2y) + 2β2 d q A3 /(yt 2 ) + 2CD A/(4β1 ), ε/(2x) + ε/(2y) + αωt/(β1 + ρ)}.

Reduced-Order Extrapolation Finite Difference Schemes Chapter | 1

39

By (1.4.19)–(1.4.22), we obtain un ∞ + v n ∞ + S n ∞ + zbn ∞  (1 +  )(un−1 ∞ + v n−1 ∞ + S n−1 ∞ + zbn−1 ∞ )     2gt 2gt Nγ Nγ + + Z 0 ∞ exp + , (1.4.23) x y 2x 2y where n = 1, 2, · · · , N . By summing (1.4.23) from 1 to n and using the discrete Gronwall lemma, Lemma 1.4.1, we obtain un ∞ + v n ∞ + S n ∞ + zbn ∞  (u0 ∞ + v 0 ∞ + S 0 ∞ + zb0 ∞ ) exp(n )     Nγ Nγ 2gnt 2gnt 0 + Z ∞ exp + exp(n ), + x y 2x 2y n = 1, 2, · · · , N. (1.4.24) As the time interval [0, T ] is finite, the RHS of (1.4.24) is bounded. Thus, from the stability theories of FD schemes (see Theorem 1.1.1 or [34,76]) and (1.4.24), we conclude that the classical FD scheme (1.4.8)–(1.4.12) is locally stable, and so its FD solutions are convergent. By Taylor’s formula, expanding (1.4.8) and (1.4.11) as well as (1.4.12) at the reference point (xj , yk , tn ), (1.4.9) at the reference point (xj + 1 , yk , tn ), and 2 (1.4.10) at the reference point (xj , yk+ 1 , tn ), or from the approximations to the 2 difference quotient to the derivatives, we obtain the error estimates (1.4.13). Remark 1.4.1. The classical FD scheme (1.4.8)–(1.4.12) is only first-order accurate in time. If one wants to obtain higher-order time approximation accuracy, it is necessary to change the time difference coefficients on the LHSs in (1.4.8)–(1.4.12) into higher order (for example, central difference or secondorder difference). Remark 1.4.2. The Coriolis constant f , the gravity acceleration g, the viscosity coefficients γ and A, the bottom drag coefficient CD , the sand diffusion coefficient ε, the falling speed of suspended sediment particles ω, the sediment mass transport velocity vc , the sediment diameter d, and empirical constants K, m, n, p, and q, boundary value functions Z0 (x, y, t), u0 (x, y, t), v0 (x, y, t), S0 (x, y, t), and zb0 (x, y, t), initial value functions Z 0 (x, y), u0 (x, y), v 0 (x, y), S 0 (x, y), and zb0 (x, y), the time step increment t, and the spatial step increments x and y are the requisites for solving the classical FD solun , Z n , and zn tions un 1 , v n 1 , Sj,k j,k b,j,k (0  j  J , 0  k  K, 1  n  j + 2 ,k

j,k+ 2

N ) for the 2D SWEs with sediment concentration through the FD schemes (1.4.8)–(1.4.12).

40 Proper Orthogonal Decomposition Methods for Partial Differential Equations

1.4.3 Establishment of the POD-Based Reduced-Order Extrapolating Finite Difference Scheme To formulate the POD basis, we extract the first L sequence of solutions l , Z l , zl (ul 1 , v l 1 , Sj,k j,k b,j,k ) (l = 1, 2, · · · , L) for the classical FD scheme j + 2 ,k

j,k+ 2

(1.4.8)–(1.4.12) as snapshots and set uli = ul

j + 12 ,k

, vil = v l

j,k+ 12

l , , Sil = Sj,k

l , and zl = zl Zil = Zj,k bi b,j,k (i = kJ + j + 1, 1  i  m, m = (J + 1)(K + 1), 0  j  J − 1, 0  k  K − 1), respectively. Thus, we formulate five m × L matrices Ar = (ril )m×L (r = u, v, S, Z, zb ) denoted by

⎛ ⎜ ⎜ ⎜ Ar = ⎜ ⎜ ⎝

r11

r12

···

r1L

r21 .. .

r22 .. .

··· .. .

r2L .. .

1 rm

2 rm

···

L rm

⎞ ⎟ ⎟ ⎟ ⎟ , r = u, v, S, Z, zb . ⎟ ⎠

It is obvious that the number of mesh points m is much larger than that of the L snapshots extracted. Thereby, the degree m for matrices Ar ATr is much larger than the degree L for matrices ATr Ar , but their positive eigenvalues are identical. Therefore, we may first find the eigenvalues λr1  λr2  · · ·  λr M˜ r > 0 (M˜ r = rankAr ) for matrices ATr Ar and corresponding eigenvectors ϕ rj . Then, by the relation  φ rj = Ar ϕ rj / λrj , j = 1, 2, . . . , M˜ r , r = u, v, S, Z, zb , we can obtain the eigenvectors φ rj (j = 1, 2, . . . , M˜ r ) corresponding to the nonzero eigenvalues for matrix Ar ATr (r = u, v, S, Z, zb ). We use the first Mr (0 < Mr  M˜ r  L) columns of the eigenmatrices U r = (φ r1 , φ r2 , · · · , φ r M˜ r ) to formulate five sets of orthonormal POD bases (see [113, 118]) r = (φ r1 , φ r2 , · · · , φ rMr ) (r = u, v, S, Z, zb ). Set n T r nm = (r1n , r2n , · · · , rm ) , n = 1, 2, · · · , N, r = u, v, S, Z, zb .

(1.4.25)

From (1.1.12) in Section 1.2.2, we easily obtain the following error estimates: r lm − r Tr r lm 2 

 λr(Mr +1) , l = 1, 2, · · · , L, r = u, v, S, Z, zb , (1.4.26)

2 1/2 is the standard l 2 norm of a vector a = (a , a , · · · , where a2 = ( m 1 2 i=1 ai ) T am ) . Thus, the classical FD scheme (1.4.8)–(1.4.12) can be rewritten in the

Reduced-Order Extrapolation Finite Difference Schemes Chapter | 1

41

following vector form: !T n+1 n+1 n+1 n+1 , v , S , Z , z un+1 m m m m b,m  n n n n n T = um , v m , S m , Z m , zb,m + F˜ (unm , v nm , S nm , Z nm , znb,m ), (1.4.27) where n = 1, 2, · · · , N − 1 and F˜ is defined by the classical FD scheme (1.4.8)–(1.4.12). In fact, in order to facilitate the programming computation, when we solve the classical FD scheme, it is also usually necessary to write the classical FD scheme in vector form. Let  ∗n ∗n ∗n ∗n ∗n T um , v m , S m , Z m , zb,m !T = u β nMu , v β nMv , R β nMS , Z β nMZ , zb β nMz , (1.4.28) b

∗n ∗n ∗n T where r ∗n m = (r1 , r2 , · · · , rm ) (r = u, v, S, Z, zb ) are five column vectors ∗n corresponding to r (r = u, v, S, Z, zb ), respectively. If we substitute u∗n m , vm , ∗n n ∗n n ∗n n n n S m , Z m , and zb,m in (1.4.28) for um , v m , S m , Z m , and zb,m in (1.4.27) (n = 0, 1, 2, · · · , N ) and note that five matrices r (r = u, v, S, Z, zb ) are formed with the orthonormal vectors, then we obtain the PODROEFD scheme for the 2D shallow water equation with sediment concentration as follows:

β nMu = Tu unm , β nMv = Tv v nm , β nMS = TS S nm , β nMZ = TZ Z nm , β nMz = Tzb znb,m , 1  n  L,

(1.4.29)

b

n+1 n+1 n+1 n+1 β n+1 Mu , β Mv , β MS , β MZ , β Mz b

!T

= β nMu , β nMv , β nMS , β nMZ , β nMz

b

!T

+

(Tu , Tv , TS , TZ , Tzb )T F˜ (u β nMu , v β nMv , S β nMS , Z β nMZ , zb β nMz ), b

L  n  N − 1.

(1.4.30)

The system of Eqs. (1.4.29)–(1.4.30) has only Mu + Mv + MS + MZ + Mzb (Mu , Mv , MS , MZ , Mzb L m) degrees of freedom at each time level and has no repetitive computations. After we have gotten β nMu , β nMv , β nMS , β nMZ , and β nMz from (1.4.29) and b (1.4.30), we will obtain the PODROEFD solutions as follows: n ∗n n u∗n m = u β Mu , v m = v β Mv , n z∗n b,m = zb β Mz , n = 1, 2, · · · , N. b

(1.4.31)

Further, we obtain the component forms of the PODROEFD solutions, denoted ∗n ∗n ∗n ∗n ∗n ∗n ∗n ∗n by u∗n 1 = u∗n 1 = vi , Sj,k = Si , Zj,k = Zi , and zb,j,k = zb,i (0  i ,v j + 2 ,k

j,k+ 2

j  J − 1, 0  k  K − 1, i = k(J + 1) + j + 1, 1  i  m = (K + 1)(J + 1)).

42 Proper Orthogonal Decomposition Methods for Partial Differential Equations

Remark 1.4.3. It is easily seen that the classical FD scheme (1.4.8)–(1.4.12) at each time level contains 5m unknowns, whereas the system of Eqs. (1.4.29)– (1.4.31) at each time level (when n > L) contains only Mu + Mv + MS + MZ + Mzb unknowns (Mu , Mv , MS , MZ , Mzb L m, for example, in Section 1.4.6, L = 20, Mu = Mv = MS = MZ = Mzb = 6, but m = 7000 or 25 × 106 ). Therefore, the PODROEFD scheme (1.4.29)–(1.4.31) contains much fewer degrees of freedom and does not involve repetitive computations. Here, we extract the snapshots from the first L classical FD solutions, but if we are faced with real-world problems we may instead extract snapshots from the samples of experiments of the relevant physical system’s trajectories.

1.4.4 Error Estimates for the POD-Based Reduced-Order Extrapolating Finite Difference Solutions In the following, we deduce the error estimates of the solutions for the PODROEFD scheme. By (1.4.31), we may write the PODROEFD scheme (1.4.29) and (1.4.30) in the following vector form: T n ∗n T n ∗n T n u∗n m = u u um , v m = v v v m , S m = S S S m , T n ∗n T n Z ∗n (1.4.32) m = Z Z Z m , zb,m = zb zb zb,m , n = 1, 2, · · · , L, !T , v ∗n+1 , S ∗n+1 , Z ∗n+1 , z∗n+1 u∗n+1 m m m m b,m  ∗n ∗n ∗n ∗n ∗n T ∗n ∗n ∗n ∗n = um , v m , S m , Z m , zb,m + F˜ (u∗n m , v m , S m , Z m , zb,m ),

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

(1.4.33)

whose stability conditions are also the same as those in (1.4.27), i.e., (1.4.8)– ∗n ∗n T (1.4.12). Set en = (unm , v nm , p nm )T − (u∗n m , v m , p m ) . By (1.4.26) and (1.4.32), we have en 2 = (unm , v nm , S nm , Z nm , znb,m )T − (u Tu unm , v Tv v nm , S TS S nm , Z TZ Z nm , zb Tzb znb,m )T 2     λu(Mu +1) + λv(Mv +1) + λS(MS +1) "  + λZ(MZ +1) + λzb (Mzb +1) , n = 1, 2, · · · , L. (1.4.34) By (1.4.27), (1.4.33), and (1.4.34), we obtain en 2  en−1 2 + δen−1 2 = (1 + δ)en−1 2  ···  (1 + δ)n−L eL 2

Reduced-Order Extrapolation Finite Difference Schemes Chapter | 1

43

#    (1 + δ)n−L λu(Mu +1) + λv(Mv +1) + λS(MS +1) $ "  + λZ(MZ +1) + λzb (Mzb +1) , n = L + 1, · · · , N, (1.4.35) where δ = t max{γ , A, ε}/min{x 2 , y 2 } and δ  1/4 under the stability conditions (1.4.8)–(1.4.12) and (1.4.29)–(1.4.31). Synthesizing the abovementioned discussions yields the error estimates of the reduced-order FD solutions for the PODROEFD scheme (1.4.29)–(1.4.31), as follows. Theorem 1.4.3. If (unm , v nm , S nm , Z nm , znb,m )T (n = 1, 2, · · · , N) are the solution vectors formed from the classical FD scheme (1.4.8)–(1.4.12) and ∗n ∗n ∗n T ∗n (u∗n m , v m , S m , Z m , zb,m ) (n = 1, 2, · · · , N) are the reduced-order FD solutions for the PODROEFD scheme (1.4.29)–(1.4.31), we have the following error estimates: ∗n ∗n ∗n ∗n (unm , v nm , S nm , Z nm , znb,m ) − (u∗n m , v m , S m , Z m , zb,m )2     C(δ n ) λu(Mu +1) + λv(Mv +1) + λS(MS +1) $ "  + λZ(MZ +1) + λzb (Mzb +1) , n = 1, 2, · · · , N,

where C(δ n ) = 1 (1  n  L), C(δ n ) = (1 + δ)n−L (L + 1  n  N ), and δ = t max{γ , A, ε}/min{x 2 , y 2 }. It is well known that the absolute value of each component of the vector is not larger than the norm of the vector. Thus, by combining (1.4.13) with Theorem 1.4.3, we obtain the error estimates between the accurate solution for the 2D SWEs and the PODROEFD solutions as follows. Theorem 1.4.4. The error estimates between the accurate solution for the 2D SWEs and the reduced-order FD solutions obtained from the PODROEFD scheme (1.4.29)–(1.4.31) are as follows: ∗n | + | v(xj , yk+ 1 , tn ) − vj,k+ | u(xj + 1 , yk , tn ) − u∗n 1 | j + 1 ,k 2

2

2

2

∗n ∗n + | S(xj , yk , tn ) − Sj,k | + | Z(xj , yk , tn ) − Zj,k | #  ∗n + | zb (xj , yk , tn ) − zb,j,k |= O C(δ n ) λu(Mu +1) + λv(Mv +1) $ ! "   + λS(MS +1) + λZ(MZ +1) + λzb (Mzb +1) , t, x 2 , y 2 ,

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

(1.4.36)

Remark 1.4.4 (How to renew the POD basis). The " error estimate terms     λu(Mu +1) + λv(Mv +1) + λS(MS +1) + λZ(MZ +1) + λzb (Mzb +1) in Theorems 1.4.3 and 1.4.4 are caused by the POD-based reduced-order for the classical FD scheme, which could be used as a criterion for the selection of the number

44 Proper Orthogonal Decomposition Methods for Partial Differential Equations

of the POD bases, i.e., it is necessary to choose Mu , Mv , MS ,"MZ , and Mzb     such that λu(Mu +1) + λv(Mv +1) + λS(MS +1) + λZ(MZ +1) + λzb (Mzb +1) = O(t, x 2 , y 2 ). Whereas C(δ n ) = (1 + δ)n−L (L + 1  n  N ) are caused by extrapolating iteration, which can be used as a criterion for the renewal of the     POD basis, i.e., if C(δ n )[ λu(Mu +1) + λv(Mv +1) + λS(MS +1) + λZ(MZ +1) + " λzb (Mzb +1) ] > max(t, x 2 , y 2 ), it is necessary to update the POD basis.

1.4.5 Algorithm Implementation for the POD-Based Reduced-Order Extrapolating Finite Difference Scheme The algorithm implementation for the PODROEFD scheme (1.4.29)–(1.4.31) involves the following five steps. Step 1. Classical FD computation and extraction of snapshots Solving the classical FD scheme (1.4.8)–(1.4.12) at the first few L steps (in the following, say, take L = 20) yields the classical FD solutions un 1 , v n 1 , j + 2 ,k

j,k+ 2

n , Z n , and zn Sj,k j,k b,j,k (0  j  J, 0  k  K, 1  n  L) and further forms a l }L (1  i  m) with L × m elements, where set of snapshots {uli , vil , Sil , Zil , zb,i l=1 n , Z n = Z n , and zn = zn uni = un 1 , vin = v n 1 , Sin = Sj,k i j,k b,i b,j,k (i = kJ + j + 2 ,k

j,k+ 2

j + 1, 1  i  m, m = (J + 1)(K + 1), 0  j  J − 1, 0  k  K − 1). Step 2. Snapshot matrices Ar and eigenvalues of ATr Ar Formulate the snapshot matrices Ar = (ril )m×L (r = u, v, S, Z, zb ) and compute the eigenvalues λr1  λr2  · · ·  λr M˜ r > 0 (M˜ r = rank Ar ) and the eigenvectors ϕ rj (j = 1, 2, · · · , M˜ r , r = u, v, S, Z, zb ) of matrices ATr Ar . Step 3. Choice of POD basis For the tolerance error μ = O(t, x 2 , y 2 ), decide the numbers Mr (Mr    M˜ r , r = u, v, S, Z, zb ) of POD bases such that λu(Mu +1) + λv(Mv +1) + "   λS(MS +1) + λZ(MZ +1) + λzb (Mzb +1)  μ and formulate the POD bases  r = (φ r1 , φ r2 , · · · , φ rMr ) (where φ sj = Ar ϕ rj / λrj , j = 1, 2, · · · , Mr , r = u, v, S, Z, zb ). Step 4. Solve/compute the PODROEFD model Solving the PODROEFD scheme (1.4.29)–(1.4.31) yields the reduced-order ∗n ∗n ∗n ∗n ∗n ∗n ∗n ∗n solution vectors u∗n m = (u1 , u2 , · · · , um ), v m = (v1 , v2 , · · · , vm ), S m = ∗n ∗n ∗n ∗n ∗n ∗n ∗n ∗n ∗n ∗n (S1 , S2 , · · · , Sm ), Z m = (Z1 , Z2 , · · · , Zm ), and zb,m = (zb,1 , zb,2 , · · · , ∗n ), which further leads to the component forms u∗n ∗n ∗n zb,m = u∗n 1 1 = vi , i ,v j + 2 ,k

j,k+ 2

∗n = S ∗n , Z ∗n = Z ∗n , and z∗n ∗n Sj,k i j,k i b,j,k = zb,i (0  j  J , 0  k  K, i = k(J + 1) + j + 1, 1  i  m = (K + 1)(J + 1)).

Reduced-Order Extrapolation Finite Difference Schemes Chapter | 1

45

Step 5. Check accuracy and renew POD basis to continue Set δ = t max{γ , A, ε}/min{x 2 , y 2 }. If   (1 + δ)n−L λu(Mu +1) + λv(Mv +1) $ "   + λS(MS +1) + λZ(MZ +1) + λzb (Mzb +1)  μ, ∗n ∗n ∗n ∗n ∗n ∗n ∗n ∗n ∗n ∗n then u∗n m = (u1 , u2 , · · · , um ), v m = (v1 , v2 , · · · , vm ), S m = (S1 , S2 ,· · ·, ∗n ∗n ∗n ∗n ∗n ∗n ∗n ∗n ∗n Sm ), Z m = (Z1 , Z2 , · · · , Zm ), and zb,m = (zb,1 , zb,2 , · · · , zb,m ) (n = 1, 2, · · · , N) are just the solution vectors for the PODROEFD scheme (1.4.29)– (1.4.31) that satisfy the desirable accuracy. Else, i.e., if

  (1 + δ)n−L λu(Mu +1) + λv(Mv +1) $ "   + λS(MS +1) + λZ(MZ +1) + λzb (Mzb +1) > μ, l ) = (r ∗l , r ∗l , · · · , r ∗l )(r = u, v, S, Z, z ; l = n − L, n − L − set (r1l , r2l , · · · , rm b m 1 2 1, · · · , n − 1); return to Step 2.

− u∗n Remark 1.4.5. Step 5 could also be updated as follows: if u∗n−1 m m 0  ∗n ∗n+1 ∗n−1 ∗n ∗n ∗n+1 ∗n−1 ∗n um − um 2 , v m − v m 0  v m − v m 2 , Z m − Z m 2  Z ∗n m − ∗n−1 ∗n−1 ∗n ∗n ∗n+1 ∗n ∗n Z ∗n+1  , S − S   S − S  , z − z   z 2 2 m m m 2 m m b,m 2 b,m − b,m ∗n+1 ∗n ∗n ∗n ∗n ∗n zb,m 2 (n = L, L + 1, · · · , N − 1), then (um , v m , Z m , S m , zb,m ) (n = 1, 2, · · · , N) are the reduced-order solution vectors for PODROEFD scheme (1.4.29)–(1.4.31) satisfying the desirable accuracy. Else, namely, if u∗n−1 − m ∗n ∗n+1 ∗n−1 ∗n ∗n ∗n+1 ∗n−1 ∗n  < u −u  , v −v  < v −v  , Z −Z  u∗n 2 2 m 0 m m m m 0 m m m m 2< ∗n−1 ∗n+1 ∗n−1 ∗n ∗n ∗n+1 ∗n Z ∗n − Z  , S − S  < S − S  , or z − z 2 2 m m m m 2 m m b,m 2 < b,m ∗n+1 ∗n i ∗i zb,m − zb,m 2 (n = L, L + 1, · · · , N − 1), let r m = r m (i = n − L, n − L − 1, · · · , n − 1, r = u, v, Z, S, zb ); return to Step 2.

1.4.6 Some Numerical Experiments for the 2D Shallow Water Equation With Sediment Concentration In the following, we present two sets of numerical experiments that demonstrate the feasibility and efficiency of the PODROEFD scheme for the 2D SWEs with sediment concentration.

1.4.6.1 Simulation of Sediment Transport and Flow in an Estuary ¯ = {(x, y) : 23 − 23x/25  y  27 + The computational domain is taken as 33x/25, 0  x  25} ∪ {(x, y) : 25  x  40, 0  y  50} (the unit of x and y is km). The depth at the entrance is 10 m (i.e., Z0 |x=0 = 0.01 km). The sediment thickness at the entrance is 2 m (i.e., zb0 |x=0 = 0.002 km). Velocity u0 of fluid in the x-direction at the entrance is 2 m/s (i.e., u0 |x=0 = 7.2 km/h), but v0 = 0.

46 Proper Orthogonal Decomposition Methods for Partial Differential Equations

FIGURE 1.4.2 The states of velocity, sediment concentration, and water depth in the delta of an estuary at the end of the first year.

Reduced-Order Extrapolation Finite Difference Schemes Chapter | 1

47

FIGURE 1.4.3 The states of velocity, sediment concentration, and water depth in the delta of an estuary at the end of the third year.

48 Proper Orthogonal Decomposition Methods for Partial Differential Equations

FIGURE 1.4.4 The states of velocity, sediment concentration, and water depth in the delta of an estuary at the end of the fifth year.

Reduced-Order Extrapolation Finite Difference Schemes Chapter | 1

49

FIGURE 1.4.5 The changes of the dominant 20 eigenvalues for a short period and the first dominant 20 eigenvalues for the full period of velocity (A, B), sediment concentration (C), and water depth in the delta of an estuary (D).

The sediment concentration in water flow is 1.2 kg/m3 (i.e., S0 = S 0 = 1.2 × 10−3 kg/km3 ). The change of bottom topography is every 100 km falling 1 m along the flow direction (i.e., zb0 = zb0 = 10−5 x + 2, 0  x  40). The bilateral boundaries of water flow are two solid borders, i.e., u0 = v0 = 0 on the set {(x, y) : y = 23 − 23x/25, 0  x  25} ∪ {(x, y) : y = 27 + 33x/25, 0  x  25} ∪ {(x, 0) : 25  x  40} ∪ {(x, 50) : 25  x  40}. The time step is t = 3600 s = 1 h. The spatial steps are x = y = 200 m = 0.2 km. We refer to [191], taking f = 1.1 × 10−4 , γ = 0.001, A = 7.5 × 10−3 , CD = 0.01, d = 0.001, ω = 0.01, vc = 0, α = 0.3, K = 0.35,  = 5, m = 0.92, p = −0.25, q = 0.25, and ρ = 1.5 × 103 . By means of the classical FD scheme (1.4.8)–(1.4.12), we obtain the clasn , and Z n of the velocity u in the sical FD solutions un 1 and v n 1 , Sj,k j,k j + 2 ,k

j,k+ 2

x-direction and v in the y-direction, the sediment concentration S, and the water depth Z (since the change of zb is very small, it has not been described) when n = 8760, 26 280, and 43 800 (namely, at the ends of the first year, third year, and fifth year, respectively), which are depicted graphically in the graphs on the left column in Figs. 1.4.2, 1.4.3, and respectively. √ √ 1.4.4,√ √ The computation is done with λu7 + λv7 + λS7 + λZ7  4.5 × 10−3 with L = 20. Thus, it is only necessary to choose the first six POD bases. The changes (black link lines) of the dominant 20 singular eigenvalues in Fig. 1.4.5 also verifies the fact. In addition, by comparing the dominant 20 singular eigenvalues (see Fig. 1.4.5) for a short period (on the first time interval 0  t  20 h)

50 Proper Orthogonal Decomposition Methods for Partial Differential Equations

FIGURE 1.4.6 The changes of relative deviations of numerical solutions of velocity (A, B), sediment concentration (C), and water depth in the delta of estuary (D) during 5 years.

and those for the full period (on the total time interval 0  t  43800 h), we find that the dominant singular eigenvalues for the short period are smaller than those for the full period. It implies that, if one takes the same modes (for example, six POD bases), the accuracy of solutions for the PODROEFD scheme (1.4.29)–(1.4.31) is higher than that for the usual POD FD scheme, where the usual POD FD scheme implies the need to use all classical numerical solutions on the time interval [0, T ] to form the snapshots and formulate the POD basis, and then to repeat computation of the POD-based reduced-order numerical solutions on the entire time interval [0, T ], as in [107,113]. Thus, we have obtained the POD-based reduced-order FD solutions of the velocity u in the x-direction and v in the y-direction, the sediment concentration S, and the water depth Z (since the change of zb is very small, it is omitted) when n = 8760, 26 280, and 43 800 (namely, at the ends of the first year, third year, and fifth year) by five steps of the algorithm implementation for the PODROEFD scheme (1.4.29)–(1.4.31) in Section 1.4.5. In the process, one needs to update the POD basis four times. The solution states are depicted graphically in the right column in Figs. 1.4.2, 1.4.3, and 1.4.4, respectively. Corresponding graphics in the left and right columns in Figs. 1.4.2, 1.4.3, and 1.4.4 exhibit close similarity.

N k k By using the relative deviations formula [r k − N k=1 (r /N)]/ k=1 (r /N) (r = u, v, S, Z, zb ), we have computed the relative deviations of the POD-based reduced-order FD solutions. The relative deviations of the POD-based reducedorder FD solutions on the starting time interval [0, T0 ] are slightly larger than those of the classical FD solutions. But since the PODROEFD scheme at each

Reduced-Order Extrapolation Finite Difference Schemes Chapter | 1

51

TABLE 1.4.1 RMSEs Between the Usual POD FD and the PODROEFD Solutions u

v

S

Z

N = 8760

1.48E−4

1.24E−4

1.37E−4

1.68E−4

N = 26 280

1.12E−3

1.35E−3

1.64E−3

1.72E−3

N = 26 280

4.38E−3

5.23E−3

6.36E−3

6.82E−3

time level only has 5 × 6 degrees of freedom, whereas the classical FD scheme has 5 × 7000 degrees of freedom, i.e., the number of degrees of freedom for the PODROEFD scheme is much smaller than that for the classical FD scheme, the PODROEFD scheme is efficient. After a specific time, the relative numerical deviations of solutions of the PODROEFD scheme are fewer than those of the classical FD scheme (see Fig. 1.4.6). In fact, Fig. 1.4.6 exhibits the TE accumulation on 0  t  5 (years), where the relative deviations of the classical FD solutions are far larger than those of the reduced-order solutions obtained from the PODROEFD scheme. According to the changes of relative deviations of classical FD solutions, the classical FD solutions appear to have bigger deviations than the POD-based reduced-order solutions do after some computing steps. Also, the error accumulation rate of the PODROEFD scheme is very slow, so it can continue to simulate the evolutions of water flow for a longer period. We have also found that the PODROEFD scheme is highly effective in manifesting the effect of sediment concentration. In order to quantify the efficiency of the PODROEFD scheme, we use the root mean square error (RMSE) and the correlation coefficient (CORCOE) between the usual POD FD solutions and the PODROEFD solutions at the ends of the first year, third year, and fifth year. RMSE and CORCOE are, respectively, obtained by the following formulas: ' ( N (1  n |2 , r = u, v, S, Z, j = 1, 3, 5, RMSE(rj ) = ) | r˜jn − rdj N n=1

N

CORCOE(rj ) = *

n=1 N

n=1

n − r¯ n ) (˜rjn − r¯˜jn )(rdj dj

(˜rjn

− r¯˜jn )2

N

n (rdj n=1

, r = u, v, S, Z, j = 1, 3, 5, n )2 − r¯dj

where r˜j n (r = u, v, S, Z, j = 1, 3, 5) are the usual POD solutions at the end of n are the PODROEFD solutions at the end of the j th year, and the j th year, rdj N = 8760, 26 280, and 43 800. Tables 1.4.1 and 1.4.2 are, respectively, RMSEs and CORCOEs between the usual POD FD solutions and the PODROEFD solutions at the ends of the first year, third year, and fifth year, i.e., N = 8760, 26 280, and 43 800 with six POD

52 Proper Orthogonal Decomposition Methods for Partial Differential Equations

TABLE 1.4.2 CORCOEs Between the Usual POD FD and the PODROEFD Solutions u

v

S

Z

N = 8760

1.57E−4

2.38E−4

2.46E−4

2.59E−4

N = 26 280

1.46E−6

2.35E−6

2.65E−6

2.71E−6

N = 26 280

1.43E−8

2.24E−8

2.58E−8

2.86E−8

FIGURE 1.4.7 The left and right graphics are the classical FD solutions and the PODROEFD solutions of the dam-break flow and the sediment concentration at t = 1 s, respectively.

bases. Table 1.4.1 shows that the numerically computed RMSEs are consistent with theoretical errors even if they increase with time step numbers. Table 1.4.2 also shows that the CORCOEs of the numerical solutions for two cases of the usual POD FD solutions and the PODROEFD solutions get smaller and smaller with time increasing, which is reasonable since the PODROEFD scheme only takes the first 20 solutions as snapshots. However, the errors are within the tolerance range. Therefore, the PODROEFD scheme is an improvement over the usual POD FD scheme. By comparing the classical FD scheme with the PODROEFD scheme containing six bases in implementing the numerical simulations when t = 5 years, we have found that the classical FD scheme at each time level has 5 × 7000 de-

Reduced-Order Extrapolation Finite Difference Schemes Chapter | 1

53

FIGURE 1.4.8 The left and right graphics are the classical FD solutions and the PODROEFD solutions of the dam-break flow and the sediment concentration at t = 3 s, respectively.

grees of freedom and requires 120 minutes computing time on a ThinkPad E530 PC, whereas the PODROEFD scheme with six POD bases at each time level only involves 5 × 6 degrees of freedom and the corresponding time is about 30 seconds on the same PC, i.e., the computing time of the classical FD scheme is about 240 times that of the PODROEFD scheme with six POD bases. We have also shown that the PODROEFD scheme can greatly reduce the TE accumulation in the process, diminish the calculation load, save time of calculations, and improve the accuracy of numerical solutions.

1.4.6.2 Simulation of Dam-Break Flow The dam-break flow, standing for an uncontrolled release of water when a vertical barrier is suddenly removed, provides one of the simplest available models for several important phenomena, such as break-out floods, sheet flow events, and the formative stages of lahars or debris flows. An idealized model of the dam-break flow may be described as follows: a barrier at x = 50 and 0  y  100 divides fluids of different depths (10 m and 5 m), until at time t = 0 a gate with width 15 m, on x = 50 and 50  y  75, of the barrier is instantaneously removed, and fluid (depth 10 m) floods into the shallower (depth 5 m) region. Thus, the computational domain for the dam-

54 Proper Orthogonal Decomposition Methods for Partial Differential Equations

FIGURE 1.4.9 The left and right graphics are the classical FD solutions and the PODROEFD solutions of the dam-break flow and the sediment concentration at t = 5 s, respectively.

¯ = [0, 100]×[0, 100]; here break flow is a square with area 100×100 m2 , i.e., the water depths are 10 m on the subdomain [0, 50] × [0, 100] and 5 m on the n = 0 in (1.4.12)–(1.4.29) subdomain [50, 100] × [0, 100]. Since zb = 0, zb,j,k ∗n and the corresponding zb,m = 0 in (1.4.29)–(1.4.31). In order to solve the classical FD scheme (1.4.12)–(1.4.29) and the PODROEFD scheme for the 2D SWEs including sediment concentration, we proceed by taking t = 0.01 s, x = y = 0.02 m, and we use the parameter values f = 1.1 × 10−4 , γ = 0.001, A = 7.5 × 10−3 , CD = 0.01, ω = 0.01, α = 0.3, K = 0.35, m = 0.92, and ρ = 1.5 × 103 (see [190,191]). By the classical FD scheme (1.4.8)–(1.4.12), we obtain the solutions of the dam-break flow and the sediment concentration (with zb = 0, which will not come into play) when n = 100, 300, and 500 (namely, at 1 s, 3 s, and 5 s), depicted graphically in the left column in Figs. 1.4.7, 1.4.8, and 1.4.9, respectively. We have also obtained the classical FD solutions of u and v at t = 5 s, depicted in the two graphics in the left column in Fig. 1.4.10. √ √ 20; this is achieved by computing λu7 + λv7 + √ Similarly, √ we use L = −3 λS7 + λZ7  3.5 × 10 . Thus, it is also only necessary to choose the first six POD bases. Then, by the five steps of the algorithm implementation for the

Reduced-Order Extrapolation Finite Difference Schemes Chapter | 1

55

FIGURE 1.4.10 The left and right graphics are the classical FD solutions and the PODROEFD solutions of the dam-break flow velocity u and v at t = 5 s, respectively.

PODROEFD scheme (1.4.29)–(1.4.31) in Section 1.4.5 without the need for POD-basis renewal, we have obtained the PODROEFD solutions for the dambreak flow and sediment concentration (as zb = 0, it is omitted), for n = 100, 300, and 500 (i.e., at 1 s, 3 s, and 5 s), depicted graphically in the right column in Figs. 1.4.7, 1.4.8, and 1.4.9, respectively. We have also obtained the PODROEFD solutions of u and v at t = 5 s, depicted graphically in the right column in Fig. 1.4.10. Each pair of corresponding graphics in the left and right columns in Figs. 1.4.7, 1.4.8, 1.4.9, and 1.4.10 manifests close similarity. Since the classical FD scheme at each time level includes 4 × 25 × 106 degrees of freedom in the numerical simulation of dam-break flow, whereas the PODROEFD scheme with six POD bases at each time level involves only 4 × 6 degrees of freedom, the PODROEFD scheme has the same advantages as before. We have also shown that the PODROEFD scheme is successful in simulating the dambreak flow. Remark 1.4.6. In this section, we have used the POD method to establish the PODROEFD scheme with very few degrees of freedom for the 2D SWEs with sediment concentration. We have also provided the error estimates between the exact solution and the classical FD solutions and the ones between the accu-

56 Proper Orthogonal Decomposition Methods for Partial Differential Equations

rate solution and the PODROEFD solutions. In particular, we have provided two numerical examples to illustrate that the PODROEFD scheme is effective in finding the numerical solutions for the 2D SWEs containing sediment concentration. Our PODROEFD model for the 2D SWEs including sediment concentration is completely different from existing POD-based reduced-order models for SWEs (see, e.g., [29,151,152]) and it constitutes an improvement over and a development for these existing models and other POD-based reducedorder methods or reduced-basis ones.

1.5 CONCLUSIONS AND DISCUSSION ABOUT POD-BASED REDUCED-ORDER EXTRAPOLATION FINITE DIFFERENCE SCHEMES In this chapter, we have described the principles and methods of the PODROEFD schemes from the simple to the complex. We first introduce the PODROEFD schemes for the 2D parabolic equation so that the readers can understand well the basic principle and overall process how the POD method is used to reduce the order of the classical FD schemes, which should be helpful for the beginning readers. After the readers have gotten the basic ideas of the PODROEFD schemes, we provide the method of the PODROEFD schemes for the more complex fluid mechanics equations, i.e., the nonstationary Stokes equation, so that the readers can further understand the techniques of the PODROEFD schemes. Finally, we provide the method of the PODROEFD schemes for even more complex problems, i.e., the 2D SWEs with sediment concentration so that readers are aware that the PODROEFD schemes hold what might be called “magic” power in the reduction of the degrees of freedom for the classical FD schemes, the decrease in the TE accumulation and the computational load, and the savings of time for calculations and resource demands during the computational process. These have been our objectives for this Chapter 1. In fact, if one can construct some advantageous classical FD schemes, the PODROEFD schemes can not only maintain the advantages of one’s classical FD schemes, but also greatly lessen the degrees of freedom of the classical FD schemes so that they could alleviate the TE accumulation and the computational load, and save time for calculations and resource demands during the computational process. Therefore, the PODROEFD schemes here have demonstrated themselves to be one of the most outstanding reduced-order methods. For more examples, see references [5–7,38,40,79,81,91,94,96,102,111,113, 117,118,121,122,127,154,155,158]. Remarkably, the principles and methods of the PODROEFD schemes were developed less than 5 years ago. Their theory, methods, and applications should represent much more further work worth studying. We hope the readers can develop their own methods. If these methods are applied to real-world problems, they are sure to result in great efficiency.