Uniform Bounds with Difference Quotients for Proper Orthogonal Decomposition Reduced Order Models of the Burgers Equation

Journal of Scientific Computing (2023) 95:43 https://doi.org/10.1007/s10915-023-02160-2 Uniform Bounds with Difference Quotients for Proper Orthogonal Decomposition Reduced Order Models of the Burgers Equation Birgul Koc1·Tomás Chacón Rebollo2·Samuele Rubino2 Received: 17 June 2022 / Revised: 22 December 2022 / Accepted: 16 February 2023 / Published online: 20 March 2023 © The Author(s) 2023 Abstract In this paper, we prove uniform error bounds for proper orthogonal decomposition (POD) reduced order modeling (ROM) of Burgers equation, considering difference quotients (DQs), introduced in Kunisch and Volkwein (Numer Math 90(1):117–148, 2001). In particular, we study the behavior of the DQ ROM error bounds by considering L2() and H1 0() POD spaces and l∞(L2)and natural-norm errors. We present some meaningful numerical tests checking the behavior of error bounds. Based on our numerical results, DQ ROM errors are several orders of magnitude smaller than noDQ ones (in which the POD is constructed in a standard way, i.e., without the DQ approach) in terms of the energy kept by the ROM basis. Further, noDQ ROM errors have an optimal behavior, while DQ ROM errors, where the DQ is added to the POD process, demonstrate an optimality/super-optimality behavior. It is conjectured that this possibly occurs because the DQ inner products allow the time dependency in the ROM spaces to make an impact. Keywords Difference quotients ·Proper orthogonal decomposition ·Reduced order models ·Error analysis ·Optimality 1 Introduction Reducedordermodels(ROMs)areoneofthemostpopularlow-dimensionalsurrogatemodels to obtain the numerical simulation of linear and nonlinear systems [10,13,15–18,23,32– 35,42,44,45]. To build a low dimensional ROM, one can use the following most popular BBirgul Koc [email protected] Tomás Chacón Rebollo [email protected] Samuele Rubino [email protected] 1IMUS, Universidad de Sevilla, Sevilla, Spain 2Departamento EDAN and IMUS, Universidad de Sevilla, Sevilla, Spain 123 43 Page 2 of 27 Journal of Scientific Computing (2023) 95 :43 frameworks such as proper orthogonal decomposition (POD) [1,3–7,12,14,17,21,25, 31,33,46,47], reduced basis methods (RBM) [8,36,40], empirical interpolation method (EIM), and discrete empirical interpolation method (DEIM) [2,11,37,39]. In this work, we specifically use the POD framework to build the ROM. The ROM error analysis for the parabolic problems has been worked in [19,26,27,29, 41,43]. Difference quotients (DQs) (i.e., scaled snapshots of the form (un−un−1)/t,n= 1,...,N) were proposed by Kunisch and Volkwein in Remark 1 of [26] as a means to achieve time discretization optimality. The effect of the DQs in linear applications and the ROM error analysis for the parabolic problems are investigated in [9,28]and[19,22,26, 29], respectively. The authors in [19] provide the DQ ROM numerical results for the Burgers equation without providing any numerical analysis for the results. Thus, we emphasize that, to our knowledge, the ROM error analysis considering the DQs for nonlinear problems has never been proven. The main aim of the paper is to provide an error analysis for the ROM approximation of Burgers Eqs. (1) considering DQ, L2() and H1 0() POD space frameworks and the l∞(L2()) andnatural-norms. Weprovideuniform ROM errorbounds andpresent numerical tests in which we observe a much smaller error of the DQ ROM errors than the noDQ ones, with respect to energy kept in the ROM basis. This indicates that the DQ ROM keeps better-selected information for the same amount of energy due to the difference quotients. The rest of the paper is organized as follows: In Sect.2, we briefly describe the noDQ/DQ POD methodology and provide error bounds for the POD projection error. DQ Crank- Nicolson (CN) ROM for the Burgers Eq. (1) is defined in Sect.3. The error analysis of the DQ CN ROM is discussed in Sect.4. Specifically, in Sect.4.1, we classify the optimality type for the different ROM discretization errors. Furthermore, in Sects.4.2 and 4.3,weprovide l∞(L2)and natural-norm, i.e., l∞(L2)∩l2(H1 0)DQ ROM error bounds considering L2() and H1 0() POD bases and their optimality behavior, respectively. As a mathematical model, we use viscous Burgers equation: ⎧ ⎪ ⎨ ⎪ ⎩ ut−νuxx +uux=f,x∈, t∈(0,T), u(0,t)=u(1,t)=0,t∈(0,T), u(x,0)=u0(x), x∈. (1) In Sect.5, we provide ROM error bounds considering noDQ, DQ framework, L2(),H1 0() POD spaces, and l∞(L2), natural-norm errors. Specifically, in Sects.5.1 and 5.2, we numerically discuss the optimality behavior of the noDQ and DQ ROM errors, respectively, considering L2(),H1 0() POD spaces, and l∞(L2), natural-norm errors. In Sect.5.3, we numerically compare and discuss the noDQ and DQ ROM errors considering both POD spaces and norm errors. Finally, Sect.6presents the conclusions and future research directions. 2 Proper Orthogonal Decomposition (POD) This section builds a general POD framework with/without DQs. The construction of the noDQ/DQ POD basis is straightforward and can be summarized in the following steps: (i) We collect a snapshot data set S:={u0,u1, ..., uN}that is contained in a real Hilbert space Hand by solving (1) for equispaced parameter tj=jt,∀j=0, .., Nwhere t=T/N. (ii) We obtain noDQ/DQ orthonormal POD basis functions (usually called POD modes), i.e., {ϕ1, ..., ϕr}⊂Swith fixed r>0 value by solving the following generalized minimization 123 Journal of Scientific Computing (2023) 95 :43 Page 3 of 27 43 problem: min (ϕi,ϕj)H=δij u−Pru2 ∗(2) where the ∗norm is defined as v∗:= 1 MξN  i=0vi2 H+ξ N  i=1∂vi2 H1/2 ,(3) with the weight Mξand tuning parameter ξ, and the DQs, i.e., ∂vk,definedby ∂vk:= vk−vk−1 t,(4) and Pr:H→His the orthogonal projection onto Xr:= span{ϕi}r i=1given by Pru= r  i=1 (u,ϕ i)Hϕi,u∈H.(5) Depending on the choice of the tuning parameter ξ(we let ξtobe0or1)in(3), (7) yields the noDQ/DQ POD. To be more specific, the choice of ξ=0 discards the DQ summation in (3); thus, the minimization problem (2) is solved for just the given snapshot data set Swith the weight Wξ=N+1 to obtain the standard POD (noDQ POD) basis. On the other hand, choosing ξ=1 (taking into account the DQ summation) leads to (2)tobesolvedforthe different snapshot data set S={u0, .., uN,∂u1, ..., ∂uN}with Mξ=2N+1, which results in the DQ POD basis. In order to solve (2), one considers the eigenvalue problem Kv=λv, (6) where Kji := (uj,ui)H,i,j=1, .., Mξis the snapshot correlation matrix, λ1≥λ2≥ , ..., ≥λd>0 are the positive eigenvalues, and vk,k=1, ..., d, are the associated eigenvectors. Then, the solution of (2)isgivenby ϕi=1 √λi Mξ  j=1 (vi)juj.(7) To obtain the L2or H1 0POD space framework, one needs to choose the Hilbert space Hin (3)and(6) as either L2or H1 0, respectively. For example, if one needs to create the DQ-H01 POD, one should choose H=H1 0in (3)and(6)andξ=1in(3) to solve the minimization problem (2) and eigenvalue problem (6). Now, we provide DQ POD approximation errors in Lemma 2.1 proven in [19] by considering different norms and projections onto Xr=span{ϕi}r i=1⊂H. Furthermore, we present the uniform DQ POD projection error bounds in Theorem 2.2 provenin[24], Theorem 3.7. These results are necessary to prove DQ ROM error bounds and show their optimality behavior in Sect.3. Lemma 2.1 Let Xr=span{ϕi}r i=1⊂H,letP r:H→Hbe the orthogonal projection onto Xras defined in (5), and let d be the number of positive POD eigenvalues, where λDQ i represents the DQ POD eigenvalues for the collection UDQ ={un}N n=0∪{∂un}N n=1described 123 43 Page 4 of 27 Journal of Scientific Computing (2023) 95 :43 above. If W is a real Hilbert space with UDQ ⊂W and Rr:W→W is a bounded linear projection onto Xr,then 1 2N+1N  n=0 un−Prun  2 W+ N  n=1 ∂un−Pr∂un  2 W= d  i=r+1 λDQ iϕi2 W,(8a) 1 2N+1N  n=0 un−Rrun  2 W+ N  n=1 ∂un−Rr∂un  2 W= d  i=r+1 λDQ iϕi−Rrϕi2 W. (8b) Theorem 2.2 Let Xr=span{ϕi}r i=1⊂H,letP r:H→Hbe the orthogonal projection onto Xras defined in (5), let d be the number of positive POD eigenvalues, where λDQ irepresents the DQ POD eigenvalues for UDQ, and let tj:= jt,∀j=0, ..., Nwheret=T/N. If W is a real Hilbert space with UDQ ⊂W and Rr:W→W is a bounded linear projection onto Xr,then max 0≤k≤N  uk−Pruk   2 H≤C d  i=r+1 λDQ i,(9a) max 0≤k≤N  uk−Pruk   2 W≤C d  i=r+1 λDQ iϕi2 W,(9b) max 0≤k≤N  uk−Rruk   2 W≤C d  i=r+1 λDQ iϕi−Rrϕi2 W,(9c) where C=6max{1,T2}. Remark 2.3 The bounds in Lemma 2.1 and Theorem 2.2 are still valid if one replaces the snapshots data {u0,u1, ..., uN}, which are the continuous solution data in this paper as in [24], with finite element (FE) solutions {u0 h,u1 h, ..., uN h},wherehis the spatial discretization parameter. Furthermore, the data set used to generate the POD basis in (2) should be the same as the data set used in the POD projection error bounds in Lemma 2.1 and Theorem 2.2. Remark 2.4 The construction of all ROMs, which are used in the following sections, are obtained by using Crank-Nicolson and Galerkin time and space discretizations, respectively. However, they differ from each other based on two main criteria: noDQ/DQ and L2/H1 0POD frameworks. Thus, when we label the name of the models for brevity, we drop CN, POD, and ROM acronyms, and, for clarity, we consider noDQ/DQ and L2/H1 0acronyms. 3 Reduced Order Modeling (ROM) In this section, we present a numerical method for the Burgers Eq. (1), which is the proper orthogonal decomposition reduced order model. First,wedefinethefunctionspace X=H1 0() endowedwiththeinnerproduct(u,v)H1 0= (ux,v x)L2.Wetakeu(·,t)∈X,t∈[0,T]to be the weak solution of the weak formulation of the Burgers equation with homogeneous Dirichlet boundary conditions: (∂tu,v)L2+ν(ux,v x)L2+(uux,v)L2=(f,v)L2∀v∈X.(10) 123 Journal of Scientific Computing (2023) 95 :43 Page 5 of 27 43 Applying Crank-Nicolson and Galerkin discretizations in time and space, respectively to the weak formulation of the Burgers Eq. (10) results in the CN ROM: ∀vr∈Xr, ∂un+1 r,v rL2+ν((ur)n+1/2 x,(v r)x)L2+(un+1/2 r(ur)n+1/2 x,v r)L2=(fn+1/2,v r)L2, (11) where ∂un+1 r:= (un+1 r−un r)/t. Remark 3.1 We use the notation zn+1/2for any discrete-time function zto denote the average zn+1/2:= 1 2zn+1+zn. However, for a continuous time function, we use fn+1/2to denote f(tn+t/2). 4 Error Analysis In this section, we prove uniform error bounds for the DQ ROM for the Burgers Eq. (1). Specifically, we provide l∞(L2)and natural-norm (l∞(L2)∩l2(H1 0)) error bounds considering L2() and H1 0() POD bases and their optimality behavior in Sect.4.2 and Sect.4.3, respectively. We start the analysis by applying the CN time discretization to (10), which yields the following: ∀v∈X, ∂un+1,vL2+ν(un+1/2 x,v x)L2+(un+1/2un+1/2 x,v)L2=(fn+1/2,v)L2+τn(v), (12) with the corresponding consistency error τn(v) := ∂un+1−∂tu(tn+t/2), vL2+νuxx(tn+t/2)−un+1/2 xx ,vL2 +un+1/2un+1/2 x−u(tn+t/2)ux(tn+t/2), vL2. (13) The consistency error (13) does not depend on the fterm because of Remark 3.1.Furthermore, we assume the following regularity conditions on the continuous solution uand the terms in (13): u∈L∞(H1 0()), (14a) uttt,(utt)xx,(u(utt)x+uttux)∈L2(0,T;L2()). (14b) Now, we define the regularity constants, which are the bounds for the terms in (14b)as In,1(u):= utttL2(tn,tn+1;L2)+(utt)xxL2(tn,tn+1;L2) +u(utt)x+uttuxL2(tn,tn+1;L2), In(u):= uttt2 L2(tn,tn+1;L2)+(utt)xx2 L2(tn,tn+1;L2) +u(utt)x+uttux2 L2(tn,tn+1;L2), I(u):= uttt2 L2(0,T;L2)+(utt)xx2 L2(0,T;L2) +u(utt)x+uttux2 L2(0,T;L2). (15) Now, we subtract (11) from (12) by choosing v=vrin (12)(sinceXr⊂X), and label the discretized error en+1:= un+1−un+1 r. Then, one gets the following error equation: 123 43 Page 6 of 27 Journal of Scientific Computing (2023) 95 :43 ∀vr∈Xr, ∂en+1,v rL2+ν(en+1/2 x,(v r)x)L2+(un+1/2un+1/2 x,v r)L2 −(un+1/2 r(ur)n+1/2 x,v r)L2=τn(vr). (16) Then, we split the discretized error en+1into two parts as en+1=un+1−un+1 r=(un+1−wn+1 r)−(un+1 r−wn+1 r) =ηn+1−φn+1 r,(17) where wn+1 r:= Rrun+1is chosen as the Ritz projection of un+1on Xr(for different analyses, wrcould be chosen differently), which is defined as (u−Rru)x,(v r)xL2=0∀vr∈Xr.(18) The POD projection error (η) and the discretization error (φ)in(17) are defined as ηn+1:= un+1−wn+1 r,(19a) φn+1 r:= un+1 r−wn+1 r.(19b) During the analysis, we need a standard stability estimate of un rfor the CN scheme, which is max 0≤n≤Nun rL2≤C.(20) Now, by using the error splitting (17), one can rewrite the error equation (16)as (∂φn+1 r,v r)L2+ν((φr)n+1/2 x,(v r)x)L2=(∂ηn+1,v r)L2+ν(ηn+1/2 x,(v r)x)L2 +(un+1/2un+1/2 x,v r)L2−(un+1/2 r(ur)n+1/2 x,v r)L2−τn(vr). (21) Then, (18) leads to (ηx,(v r)x)L2=0,∀vr∈Xr; thus, the second term (ηn+1/2 x,(v r)x)L2 on the right-hand side of (21) vanishes. We continue the error analysis by choosing vr:= φn+1/2 r,then(21) is rewritten as 1 2tφn+1 r2 L2−φn r2 L2+ν(φr)n+1/2 x2 L2=(∂ηn+1,φn+1/2 r)L2 +(un+1/2un+1/2 x,φn+1/2 r)L2−(un+1/2 r(ur)n+1/2 x,φn+1/2 r)L2−τn(φn+1/2 r). (22) Now, we individually bound the terms in (22). During the analysis, Cis a generic constant that only depends on the data. By using the Cauchy-Schwarz inequality and the Young’s inequality, the first term on the right-hand side of (22), can be bounded as (∂ηn+1,φn+1/2 r)L2≤∂ηn+1L2φn+1/2 rL2 ≤1 4∂ηn+12 L2+φn+1/2 r2 L2. (23) Next, we arrange the nonlinear terms in (22) by adding and subtracting the term (un+1/2 run+1/2 x,φn+1/2 r)L2. Then, we rewrite the nonlinear terms as (un+1/2un+1/2 x,φn+1/2 r)L2−(un+1/2 r(ur)n+1/2 x,φn+1/2 r)L2 =((u−ur)n+1/2un+1/2 x,φn+1/2 r)L2+(un+1/2 r(u−ur)n+1/2 x,φn+1/2 r)L2 =(ηn+1/2un+1/2 x,φn+1/2 r)L2−(φn+1/2 run+1/2 x,φn+1/2 r)L2 +(un+1/2 rηn+1/2 x,φn+1/2 r)L2−(un+1/2 r(φr)n+1/2 x,φn+1/2 r)L2. (24) 123 Journal of Scientific Computing (2023) 95 :43 Page 7 of 27 43 Now, we individually bound nonlinear terms in (24). By using the Hölder’s inequality, the regularity condition of the continuous solution in (14a), and Young’s inequality, we can bound the first term in (24)as (ηn+1/2un+1/2 x,φn+1/2 r)L2≤ηn+1/2L2un+1/2 xL∞φn+1/2 rL2 ≤Cηn+1/2L2φn+1/2 rL2 ≤Cηn+1/22 L2+φn+1/2 r2 L2. (25) For the second term in (24), we use the Hölder’s inequality and the regularity condition of the continuous solution in (14a)as (φn+1/2 run+1/2 x,φn+1/2 r)L2≤φn+1/2 rL2un+1/2 xL∞φn+1/2 rL2 ≤Cφn+1/2 r2 L2. (26) For the third term in (24), we use the Hölder’s inequality, the standard stability estimate of un rin l∞(L2)for CN scheme in (20), the Sobolev embedding, and the Young’s inequality as (un+1/2 rηn+1/2 x,φn+1/2 r)L2≤un+1/2 rL2ηn+1/2 xL2φn+1/2 rL∞ ≤Cηn+1/2 x2 L2+C1(φr)n+1/2 x2 L2. (27) Finally, for the last nonlinear term in (24), we use the Hölder’s inequality, the stability estimate of un rin (20), the Agmon’s inequality (eq. after (45) in [26]): ϕL∞≤Cϕ1/2 L2ϕx1/2 L2,∀ϕ∈H1 0, and the Young’s inequality (p=4,q=4/3), then we get (un+1/2 r(φr)n+1/2 x,φn+1/2 r)L2≤un+1/2 rL2(φr)n+1/2 xL2φn+1/2 rL∞ ≤Cφn+1/2 r1/2 L2(φr)n+1/2 x3/2 L2 ≤Cφn+1/2 r2 L2+C2(φr)n+1/2 x2 L2. (28) To bound the consistency error (13), we use the Taylor’s theorem, the Young’s inequality, and the property (a+b+c)2≤3(a2+b2+c2)for In,1(u)in (15), then we have τn(φn+1/2 r)≤t3/2In,1(u)φn+1/2 rL2 ≤3 4t3In(u)+φn+1/2 r2 L2. (29) Now, we choose coefficient C1=C2=ν/4, and insert all bounds (23)-(29)into(22), then 1 2tφn+1 r2 L2−φn r2 L2+ν 2(φr)n+1/2 x2 L2≤Cφn+1/2 r2 L2+∂ηn+12 L2 +ηn+1/22 L2 +ηn+1/2 x2 L2+t3In(u). (30) To derive a valid error bound from (30), there should be a relation between the time step tand the viscosity coefficient νthat is explained in the following lemma. 123 43 Page 8 of 27 Journal of Scientific Computing (2023) 95 :43 Lemma 4.1 Let t<4Cν3 27 , then it holds (1−Ct)φn+1 r2 L2+νt(φr)n+1/2 x2 L2≤1+Ctφn r2 L2+Ct∂ηn+12 L2 +ηn+1/22 L2+ηn+1/2 x2 L2 +t3In(u), (31) where the constant Cabove is independent of all discretization parameters but depends on the problem data. Proof First, bounding the term φn+1/2 r2 L2yields 1 2φn+1 r2 L2+φn r2 L2in (30). Then, if one computes the constant coefficients in (23)-(29), the generic constant Cin (30) will be C:= max 1 4,un+1/2 x2 L∞ 4,un+1/2 xL∞ 2, C2 ν,27 C4 8ν3,3 4,3 2.(32) Then the constant Cdoes not depend on the ROM solution. Without loss of generality, in this paper, we assume that 27C4 8ν3dominates the terms in (32). Finally, multiplying the resulting equation with 2t, one obtains (31) with the constraint t<4Cν3 27 ,whereC=C−4. Remark 4.2 One can bound the term (φr)n+1/2 x2 L2in (27)-(28)byusingFEinverse estimates, as in Theorem 9.2 in [38]. Remark 4.3 For any φn+1 r∈Xh, which is the FE space that contains Xr, then the following FE inverse estimate holds: (φn+1 r)xL2≤Ch −1φn+1 rL2.(33) However, by using the FE inverse estimate to bound the term (φr)n+1/2 x2 L2in (27)-(28), the generic constant Cin (31) depends on the space discretization, and we eventually loose hconvergence order, assuming uniformly regular grid. Furthermore,thebound in(31)obtained by usingtheGalerkin method leadstoa restriction on the time step t(see Lemma 4.1). Using stabilized methods would allow relaxing this restriction. In this paper, our concern is analyzing the error estimates optimality with respect to the different POD setting strategies, so we will consider moderate values of ν. Remark 4.4 In this paper, we construct the POD basis by using the snapshots, which are the FE solutions; thus, ηnbounds in Lemma 2.1 and Theorem 2.2 should be expressed in terms of the FE solution data. We next bound the projection error ηnin terms of the finite element solution. We next bound the projection error ηnin terms of the finite element solution. Lemma 4.5 Let Rru and Rruhbe the Ritz projection, which is defined in (18), of the continuous solution u and FE solution uh, respectively, then the following estimates hold 1 N+1 N  n=0un−Rrun2 W≤C(h2l+t4)+1 N+1 N  n=0un h−Rrun h2 W(34a) 123 Journal of Scientific Computing (2023) 95 :43 Page 9 of 27 43 1 N N  n=1∂(un−Rrun)2 L2≤C(h2l+1+t3)+1 N N  n=1∂(un h−Rrun h)2 L2,(34b) where h and t are space and time discretization parameters, respectively, and l is the FE interpolation order. Proof We start proving with W=H1 0in (34a). By using the definition of the Ritz projection (18), we have ∇(un−Rrun)2 L2=(∇(un−Rrun), ∇(un−Rrun)), =(∇(un−Rrun), ∇(un−Rrun h)), ≤∇(un−Rrun)L2∇(un−Rrun h)L2. (35) Then, by adding and subtracting the FE solution un hin (35) and summing the resulting inequality from n=0toN,weget ∇(un−Rrun)2 L2≤∇(un−un h)2 L2+∇(un h−Rrun h)2 L2, 1 N+1 N  n=0∇(un−Rrun)2 L2≤C(h2l+t4)+1 N+1 N  n=0∇(un h−Rrun h)2 L2. (36) For W=L2case, we still revisit the definition of the Ritz projection (18) and choose vr=Rru,thenwehave ∇RruL2≤∇uL2.(37) Then, by adding and subtracting the FE solution un hin un−RrunL2, applying the Poincaré inequality, the relation in (37), and summing from n=0toNgive un−Rrun2 L2≤un−un h2 L2+Rr(un−un h)2 L2+un h−Rrun h2 L2, ≤un−un h2 L2+Cp∇Rr(un−un h)2 L2+un h−Rrun h2 L2, 1 N+1 N  n=0un−Rrun2 L2≤C(h2l+t4)+1 N+1 N  n=0un h−Rrun hL2. (38) Finally, to bound the term ∂(un−Rrun), we follow similar steps as in (38) and get the following: 1 N N  n=1∂(un−Rrun)2 L2≤C(h2l+t3)+1 N N  n=1∂(un h−Rrun h)L2.(39)  Now, in the following two sections, by using Lemma 4.1, we continue to derive and discuss DQ ROM errors in different norms considering L2() and H1 0() POD bases that will lead to different consistency error estimates. 123 43 Page 16 of 27 Journal of Scientific Computing (2023) 95 :43 Fig. 3 The behavior of the l∞(L2)and natural-norm noDQ-H01 errors For l∞(L2)noDQ-L2 error bound, applying (42a) in Lemma 4.8 to the second term, which is in (ϕi)x2 L2,in(58a) yields a coefficient CL2 inv:= √Sr2which depends on the ROM dimension r. We numerically investigate the behavior of √Sr2, and observe that it increases as r increases (see the top-right plot in Fig.2). Since the l∞(L2)noDQ-L2 error bound does meet the suboptimality criteria in Remark 4.7,l∞(L2)noDQ-L2 error bound (58a) is suboptimal. For natural-norm noDQ-L2 error bound, if l∞(L2)and l2(H1 0)parts of the natural-norm in (58a) are controlled by ϕi2 L2and (ϕi)x2 L2, respectively, then based on Definition 4.6, one can conclude that the natural-norm noDQ-L2 (58a) is optimal-I. For l∞(L2)noDQ-H01 error bound, applying (42b) in Lemma 4.8 to ϕiH1 0in (58b) yields the coefficient CH1 0 inv=M−1 r2. We numerically investigate the behavior of CH1 0 inv and observe that it increases as r increases (see the top-right plot in Fig.3). Since, the l∞(L2)noDQ-H01 error bound does meet the suboptimality criteria in Remark 4.7,l∞(L2) noDQ-H01 error bound (58b) is suboptimal. For natural-norm noDQ-H01 error bound, if l∞(L2)and l2(H1 0)parts of the naturalnorm in (58b) are controlled by ϕi2 L2and 1 =(ϕi)x2 L2, respectively, then based on Definition 4.6, one can conclude that the natural-norm noDQ-H01 (58a) is optimal-I. The behavior of the error bounds, which are theoretically derived in Sects.4.2 and 4.3, are summarized in Table 1. Considering L2and H1 0POD basis, error norms, and noDQ/DQ frameworks, we observe that for the DQ errors bounds with all cases and the natural norm noDQ-L2 and noDQ-H01 are optimal; whereas the l∞(L2)noDQ-L2 and noDQ-H01 error bounds are suboptimal. 123 Journal of Scientific Computing (2023) 95 :43 Page 17 of 27 43 Fig. 4 Comparison of the noDQ-L2 and noDQ-H01 solution plots with two different rvalues Table 1 Theoretical results: the behavior of the noDQ and DQ ROMs with L2() and H1 0() POD basis and l∞(L2)and natural norm errors l∞(L2)Natural norm noDQ-L2 Suboptimal Optimal-I Sect.4.3 Sect.4.3 (58a) in Remark 4.14 (58a) in Remark 4.14 noDQ-H01 Suboptimal Optimal-I Sect.4.3 Sect.4.3 (58b) in Remark 4.14 (58b) in Remark 4.14 DQ-L2 Optimal-II Optimal-II Sect.4.2 Sect.4.3 (43) in Theorem 4.9 (54) in Theorem 4.12 DQ-H01 Optimal-I Optimal-I Sect.4.2 Sect.4.3 (52) in Theorem 4.11 (57) in Theorem 4.13 (41a) in Definition 4.6 (41a) in Definition 4.6 5 Numerical Results In this section, we provide numerical results for the Burgers equation (1) with the following initial condition u0(x)=1,x∈(0,1/2], 0,x∈(1/2,1].(59) 123 43 Page 18 of 27 Journal of Scientific Computing (2023) 95 :43 Fig. 5 The behavior of the l∞(L2)and natural-norm DQ-L2 errors This condition generates a smooth solution at any t>0, with an infinite time gradient at t=0. This allows us to test the role of the difference quotients. To obtain full order model (FOM) data (the FOM simulation is provided in Fig.1), we solve (1) by using the finite element method considering ν=10−2,f=0, mesh size h=1/512, piecewise linear finite elements for the spatial discretization, and Crank-Nicolson time discretization. A small time step t=10−3is taken to obtain the errors due to the POD discretization. For all test cases, we compute two different absolute norm ROM errors: El∞(L2)=max 0≤k≤Nek2 L2,(60a) El∞(L2)∩l2(H1 0)=max 0≤k≤Nek2 L2+νt N−1  n=0en+1/2 x2 L2,(60b) being e=uh−ur. Remark 5.1 Since the ROM initial condition is chosen as u0 r:= Rru0, the discretization error (43)att=t0, i.e., φ0 r=0 in all the error bounds derivation in numerical results. 5.1 noDQ ROM Results In this section, we numerically discuss the behavior of the l∞(L2)and natural-norm noDQ- L2 and noDQ-H01 error bounds. Based on the noDQ error estimates in (58a)-(58b)in 123 Journal of Scientific Computing (2023) 95 :43 Page 19 of 27 43 Remark 4.14, we define the following RHS terms: noDQ-RHS1 = d  i=r+1 λnoDQ iϕi2 L2+(ϕi)x2 L2+t2+t4I(u), (61a) noDQ-RHS2 = d  i=r+1 λnoDQ i(1+ϕi2 L2)+t2+t4I(u). (61b) In the top left plot of Fig.2, we plot the noDQ-RHS1 defined in (61a)andl∞(L2), naturalnormnoDQ-L2 errorsin (58a).Weobservethat thel∞(L2)and natural-normnoDQ-L2 errors stay below the noDQ-RHS1. The top right plot shows how the scaling of Sr2,whichis defined above Lemma 4.8, changes as rchanges. For the bottom plots in Fig.2,weplot the linear regression (LR) orders for l∞(L2)and natural-norm noDQ-L2 errors, from left to right, respectively. Since the LR orders in the bottom plots are close to 1, we numerically observe that l∞(L2)and natural-norm noDQ-L2 errors are optimal. However, theoretical discussions in Remark 4.14 show that the l∞(L2)and natural-norm noDQ-L2 error bounds are suboptimal and optimal, respectively. In the top left plot of Fig.3, we plot the noDQ-RHS2 defined in (61b)andl∞(L2), naturalnorm noDQ-H01 errors in (58b). The top right plot shows how the scaling of M−1 r2,which is defined above Lemma 4.8, changes as rchanges. For the bottom plots in Fig.3,weplot the linear regression (LR) orders for l∞(L2)and natural-norm noDQ-H01 errors, from left to right, respectively. Since the LR orders in the bottom plots are close to 1, we conclude that l∞(L2)and natural-norm noDQ-H01 errors are optimal. However, theoretical discussions in Remark 4.14 show that the l∞(L2)and natural-norm noDQ-H01 error bounds are suboptimal and optimal, respectively. In Fig.4, we plot the noDQ-L2 and noDQ-H01 solutions with two different rvalues, i.e., r=5,20. For both rvalues, we observe that the noDQ-H01 yields slightly more accurate results than the noDQ-L2, especially for low rvalues such as r=5. 5.2 DQ ROM Results In this section, we numerically discuss the behavior of the l∞(L2)and natural-norm DQ ROM errors considering L2() and H1 0() POD bases. Based on the DQ-L2 error estimates in (43)and(54), we define the following RHS term: DQ-RHS1 = d  i=r+1 λDQ iϕi−Rrϕi2 L2+(ϕi−Rrϕi)x2 L2+t2+t4I(u), (62a) to discuss the behavior of the DQ-L2 error estimates in (43)and(54). In the top left plot of Fig.5, we plot the DQ-RHS1 defined in (62a)andl∞(L2), natural-norm noDQ-L2 errors in (43)and(54). We observe that the l∞(L2)and natural-norm DQ-L2 errors stay below the DQ-RHS1. The top right plot shows how the scaling of Sr2, which is defined above Lemma 4.8, changes as rchanges. For the bottom plots in Fig.5, we plot the linear regression (LR) orders for l∞(L2)and natural-norm DQ-L2 errors, from left to right, respectively. The LR order for the l∞(L2)DQ-L2 error bound is more than 1.5 (1 is considered optimal); thus, we conclude that the l∞(L2)DQ-L2 is superoptimal; whereas we theoretically prove that it is optimal (see the discussion after Theorem 4.9). For the natural-norm DQ-L2 error bound, the LR order is around 1; thus, we conclude that the natural-norm DQ-L2 is optimal. 123 43 Page 20 of 27 Journal of Scientific Computing (2023) 95 :43 Fig. 6 The behavior of the l∞(L2)and natural-norm DQ-H01 errors The behavior of the natural-norm error bound theoretically and numerically match (see the discussion after Theorem 4.12 for the theoretical conclusion). Based on the DQ-H01 error estimates in (52)and(57), we define the following RHS term: DQ-RHS2 = d  i=r+1 λDQ i1+ϕi2 L2+t2+t4I(u), (63a) to discuss the behavior of the DQ-H01 error estimates in (52)and(57). In the top left plot of Fig.6, we plot the DQ-RHS2 defined in (63a)andl∞(L2), natural-norm noDQ-H01 errors in (52)and(57). We observe that the l∞(L2)and natural-norm DQ-H01 errors stay below the DQ-RHS2. Furthermore, a sudden decrease in the l∞(L2), and natural-norm DQ-H01 errors arises when a large enough number of POD modes is attempted. We think this may be related to the sharp gradients of the solution that are well represented on the DQ ROM basis only when its dimension is large enough. The top right plot shows how the scaling of M−1 r2, which is defined above Lemma 4.8, changes as rchanges. For the bottom plots in Fig.6, we plot the linear regression (LR) orders for l∞(L2)and natural-norm DQ-H01 errors, from left to right, respectively. The LR order for the l∞(L2)DQ-H01 error bound is more than 1.5 (1 is considered optimal); thus, we conclude that the l∞(L2)DQ-H01 is superoptimal; whereas we theoretically prove that is is optimal (see the discussion after Theorem 4.11). For the natural-norm DQ-H01 error bound, the LR order is around 1; thus, we conclude that the natural-norm DQ-H01 is optimal. The behavior of the natural-norm error bound theoretically and numerically match (see the discussion after Theorem 4.13 for the theoretical conclusion). InFig.7, weplot the DQ-L2and DQ-H01solutionswith differentrvalues,i.e.,r=13,28. For r=13, the DQ-H01 solution is less accurate than the DQ-L2 one since the ROM 123 Journal of Scientific Computing (2023) 95 :43 Page 21 of 27 43 Fig. 7 Comparison of the DQ-L2 and DQ-H01 solution plots with two different rvalues Table 2 Numerical results: The behavior of the noDQ and DQ ROMs with L2() and H1 0() POD basis and l∞(L2)and natural norm errors l∞(L2)Natural norm noDQ-L2 Optimal Optimal Sect.5.1 Sect.5.1 Bottom-left plots in Fig.2Bottom-right plots in Fig.2 noDQ-H01 Optimal Optimal Sect.5.1 Sect.5.1 Bottom-left plots in Fig.3Bottom-right plots in Fig.3 DQ-L2 Superoptimal Optimal Sect.5.2 Sect.5.2 Bottom-left plots in Fig.5Bottom-right plots in Fig.5 DQ-H01 Superoptimal Optimal Sect.5.2 Sect.5.2 Bottom-left plots in Fig.6Bottom-right plots in Fig.6 dimension rdoes not exceed the threshold value, which is r=27 as can be observed in Fig.6. Furthermore, the plots show that when enough POD modes are guaranteed, the DQ-H01 solution rapidly yields an accurate solution. The behavior of the error bounds, which are numerically derived in Sects.5.1 and 5.2, are summarized in Table 2. Considering L2and H1 0basis, error norms, and noDQ/DQ frameworks, we observe that the noDQ errors bounds with all cases and the natural norm 123 43 Page 22 of 27 Journal of Scientific Computing (2023) 95 :43 Fig. 8 l∞(L2)(left) and natural-norm (right) ROM errors decay with respect to the ratio of energy kept by the ROM Rλ DQ-L2 and DQ-H01 are optimal; whereas the l∞(L2)DQ-L2 and DQ-H01 error bounds are superoptimal. 5.3 noDQ and DQ ROM Comparison In this section, we numerically compare how all noDQ and DQ ROM errors change based on (i) the ROM dimension, i.e., r, (ii) the residual eigenvalues, i.e., rλ, and (iii) the ratio of residual eigenvalues, i.e., Rλ, which are defined as rλ(i):= d  i=r+1 λi,(64a) Rλ(i):=d  i=r+1 λi/d  i=1 λi,(64b) where λrepresents the noDQ/DQ eigenvalues, which solve (6), and drepresents the number of positive eigenvalues. In Fig.8, we plot all l∞(L2)ROM errors in the left plot and all natural-norm errors in the right plot and compare the behavior of these errors with respect to the ratio of energy kept by the ROM Rλin (64b). Both plots in Fig.8show that the DQ-L2 recovers a much smaller error than the noDQ-L2 for the same ratio of energy considered. Furthermore, the DQ-H01 has a better slope than the noDQ-H01, even if the DQ-H01 decays later, after stagnation. The H01-norm stagnation is possibly due to the hard test problem we have taken with very sharp gradients. In Fig.9, we plot all l∞(L2)ROM errors in the left plot and all natural-norm errors in the right plot and compare the behavior of these errors with respect to the energy kept by the ROM rλin (64a). For both plots in Fig.9, the DQ ROM errors are over three orders of magnitude smaller than noDQ ROM errors over the same rλinterval. Based on these results, we may interpret that for a given energy, the DQ ROMs hold much more information than the noDQ ROMs. In Fig.10, we plot all l∞(L2)ROM errors in the left plot and all natural-norm errors in the right plot and compare the behavior of these errors with respect to the ROM dimension r. For almost all rvalues, in both plots in Fig.10, noDQ errors are almost lower than the 123 Journal of Scientific Computing (2023) 95 :43 Page 23 of 27 43 Fig. 9 l∞(L2)(left) and natural-norm (right) ROM errors decay with respect to the energy kept by the ROM rλ Fig. 10 l∞(L2)(left) and natural-norm (right) ROM errors decay with respect to the ROM dimension r DQ errors. The natural-norm noDQ and DQ ROM error behaviors (right plot of Fig.10)are more explicit and clearer than the l∞(L2)noDQ and DQ ROM error behaviors (left plot of Fig.10). Furthermore, in both plots in Fig.10, we observe that the noDQ-L2, noDQ-H01, and DQ-L2 errors decrease progressively until r=47 and r=55, in the left and right plot, respectively, then they stagnate. However, the DQ-H01 error in both plots has a constant error behavior until r=27, and when it guarantees enough POD modes, it shows a drastic decrease from r=27 to r=28 and after r=47 in the left plot and r=55 in the right plot, it stagnates again. From Sects.5.1 and 5.2, we numerically observe that l∞(L2)noDQ-L2 and noDQ-H01 error bounds are optimal and DQ-L2 and DQ-H01 are superoptimal. If one shows that the slopes of the DQ-L2 and DQ-H01 errors are much sharper than the noDQ-L2 and noDQ-H01 errors, then superoptimal behavior makes sense. The slope comparison of the errors in Figs.8 and 9is clearer than Fig.10. In both Fig.8and 9, we observe that the decreasing rates of the DQ-L2 errors are better than the noDQ-L2 ones, being over three orders of magnitude smaller in a wide range of Rλand rλ, respectively, until stagnation occurs, likely due to round-off errors. For the H01- norm comparison, the DQ-H01 has a better slope than the noDQ-H01, even if the DQ-H01 decays later, after stagnation. Therefore, the DQ keeps better-selected information for the 123 43 Page 24 of 27 Journal of Scientific Computing (2023) 95 :43 same amount of energy due to the difference quotients. This supports the interest in using the DQs. 6 Conclusions In this paper, we provided uniform ROM error bounds of nonlinear PDEs, considering the Burgers equation as the first preliminary step considering the DQs. Overall, we theoretically proved and numerically investigated the behavior of the DQ ROM error bounds by considering L2() and H1 0() POD spaces and l∞(L2)and natural-norm errors. Furthermore, we provided the noDQ ROM errors without complete theoretical support to make a clear and complete conclusion. The main results of this paper can be summarized as follows: (i) At the theoretical level, we derived four different DQ ROM error bounds by considering two error norms, i.e., l∞(L2) and natural-norm, and two different POD space frameworks, i.e., L2() and H1 0() POD spaces. (ii) In Sect.4.2, we theoretically proved that both l∞(L2)DQ-L2 and DQ-H01 errors are optimal with respect to the ROM discretization. (iii) In Sect.4.3, we obtained the same theoretical results obtained in Sect.4.2 for the natural-norm DQ-L2 and DQ-H01 errors. (iv) In Sect.5, we observed that the DQ ROM errors are several orders of magnitude lower than the noDQ errors in terms of the rate of energy kept by the ROM basis. In Sect.5.1, we numerically showed that all noDQ errors have optimal behaviors. The numericalconvergenceratesofthel∞(L2)noDQ-L2andnoDQ-H01 errorsarebetter thanthe theoretical ones since the l∞(L2)noDQ-L2 (58a) and noDQ-H01 error (58b) are suboptimal theoretical supports. In Sect.5.2, we numerically showed that the l∞(L2)DQ-L2 and DQ-H01 error bounds are superoptimal; whereas, the natural-norm DQ-L2 and DQ-H01 error bounds are optimal. In Sects.4.2 and 4.3, we theoretically proved that all DQ ROM error bounds are optimal. Based on these results, we conclude that the numerical convergence rates for the l∞(L2) DQ-L2 and DQ-H01 are better than the theoretical ones. Furthermore, in Sect.5.3,weprovidedFigs.8-10 to discuss the superoptimality behavior of the l∞(L2)DQ-L2 and DQ-H01 errors and the optimality behavior of the l∞(L2)noDQ- L2 and noDQ-H01 errors by comparing their slopes. We believe that considering the time dependency by the DQ inner products increases the linear regression orders. Finally, in Figs.4and 7, we compared the noDQ-L2 with noDQ-H01 and DQ-L2 with DQ-H01 solution plots to understand how the POD space framework affects the accuracy of the solution for both noDQ and DQ cases. In Fig.4, we concluded that noDQ-L2 and noDQ- H01 yield similar results; however, for a lowrvalue, the noDQ-H01 solution is slightly more accurate than the noDQ-L2. However, the comparison between the DQ-L2 and DQ-H01 is quite different than the noDQ one. In Fig.7, we observed that the DQ-H01 yields an accurate solution after enough POD modes are guaranteed but reaching the accurate results is quicker than the DQ-L2. Extending and improving the effectiveness of the DQs on the Navier-Stokes equations will be our future research direction. Acknowledgements This work has been supported by Programma Operativo FEDER Andalucia 2014-2020 Grant US-1254587 and European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie Actions, grant agreement 872442 (ARIA). Author Contributions We confirm that the manuscript has been read and approved by all named authors and that there are no other persons who satisfied the criteria for authorship but are not listed. We further confirm 123 Journal of Scientific Computing (2023) 95 :43 Page 25 of 27 43 that all have approved the order of authors listed in the manuscript of us. We confirm that we have given due consideration to the protection of intellectual property associated with this work and that there are no impediments to publication, including the timing of publication, with respect to intellectual property. In so doing we confirm that we have followed the regulations of our institutions concerning intellectual property. We understand that the Corresponding Author is the sole contact for the Editorial process (including Editorial Manager and direct communications with the office). He/she is responsible for communicating with the other authors about progress, submissions of revisions, and final approval of proofs. We confirm that we have provided a current, correct email address accessible by the Corresponding Author. Signed by all authors as follows: BK, 12/21/2022, TCR, 12/21/2022, SR, 12/21/2022 Funding Funding for open access publishing: Universidad de Sevilla/CBUA This work has been supported by Programma Operativo FEDER Andalucia 2014-2020 Grant US-1254587 and European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie Actions, Grant agreement 872442 (ARIA). Data Availability Data sharing is not applicable to this article as no datasets were generated or analyzed during the current study. Declarations Conflict of interest We wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome. Furthermore, the authors have no relevant financial or non-financial interests to disclose. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. References 1. Afanasiev, K., Hinze, M.: Adaptive control of a wake flow using proper orthogonal decomposition. In: Lecture Notes in Pure and Applied Mathematics 216, 317–332 (2001) 2. Antil, H., Heinkenschloss, M., Sorensen, D. C.: Application of the discrete empirical interpolation method to reduced order modeling of nonlinear and parametric systems. In: Reduced Order Methods for Modeling and Computational Reduction, pp 101–136. Springer (2014) 3. Azaïez, M., Rebollo, T.C., Rubino, S.: A cure for instabilities due to advection-dominance in pod solution to advection-diffusion-reaction equations. J. Comput. Phys. 425, 109916 (2021) 4. Ballarin, F., Manzoni, A., Quarteroni, A., Rozza, G.: Supremizer stabilization of POD-Galerkin approximation of parametrized steady incompressible Navier–Stokes equations. Int. J. Numer. Meth. Eng. 102, 1136–1161 (2015) 5. Banks,H.T.,delRosario,R.C.,Smith,R.C.:Reduced ordermodelfeedbackcontroldesign:Computational studies for thin cylindrical shells. Technical report, North Carolina State University. Center for Research in Scientific Computation (1998) 6. Banks, H.T., Joyner, M.L., Wincheski, B., Winfree, W.P.: Nondestructive evaluation using a reduced-order computational methodology. Inverse Prob. 16(4), 929 (2000) 7. Bergmann, M., Bruneau, C.-H., Iollo, A.: Enablers for robust pod models. J. Comput. Phys. 228(2), 516–538 (2009) 8. Chacón Rebollo, T., Delgado Ávila, E., Gómez Mármol, M. M.: Reduced basis method for the Smagorinsky model. In: Recent Develop Numer. Methods Model Reduction (2016) 9. Chen, P., Quarteroni, A., Rozza, G.: Multilevel and weighted reduced basis method for stochastic optimal control problems constrained by stokes equations. Numer. Math. 133(1), 67–102 (2016) 123