Pointwise error bounds in POD methods without difference quotients

Journal of Scientific Computing (2025) 103:24 https://doi.org/10.1007/s10915-025-02838-9 Pointwise Error Bounds in POD Methods Without Difference Quotients Bosco García-Archilla1 ·Julia Novo2 Received: 2 September 2024 / Revised: 10 January 2025 / Accepted: 17 February 2025 / Published online: 8 March 2025 © The Author(s) 2025 Abstract In this paper we consider proper orthogonal decomposition (POD) methods that do not include difference quotients (DQs) of snapshots in the data set. The inclusion of DQs have been shown in the literature to be a key element in obtaining error bounds that do not degrade with the number of snapshots. More recently, the inclusion of DQs has allowed to obtain pointwise (as opposed to averaged) error bounds that decay with the same convergence rate (in terms of the POD singular values) as averaged ones. In the present paper, for POD methods not including DQs in their data set, we obtain error bounds that do not degrade with the number of snapshots if the function from where the snapshots are taken has certain degree of smoothness. Moreover, the rate of convergence is as close as that of methods including DQs as the smoothness of the function providing the snapshots allows. We do this by obtaining discrete counterparts of Agmon and interpolation inequalities in Sobolev spaces. Numerical experiments validating these estimates are also presented. Keywords Proper orthotonal decomposition ·Error analysis ·Pointwise error estimates in time Mathematics Subject Classification 65M15 ·65M60 1 Introduction There seems to be an unclosed debate on wether it is necessary or advisable to include the difference quotients (DQs) of the snapshots (or function values) in the data set in proper orthogonal decomposition (POD) methods. These allow for a remarkable dimension reduction in large-scale dynamical systems by projecting the equations onto smaller spaces spanned by the first elements of the POD basis. This basis is extracted from the data set, originally the set of snapshots, although, as mentioned, some authors include their DQs. Inclusion of BBosco García-Archilla [email protected] Julia Novo julia.nov[email protected] 1Departamento de Matemática Aplicada II, Universidad de Sevilla, Sevilla, Spain 2Departamento de Matemáticas, Universidad Autónoma de Madrid, Madrid, Spain 123 24 Page 2 of 26 Journal of Scientific Computing (2025) 103 :24 these has been essential to obtain optimal error bounds [8,12,15,19] for POD methods. More recently, pointwise (in time) error bounds have been proved in case DQs are added to the set of snapshots [12], or if the DQs are complemented with just one snapshot [4]. In fact, counterexamples are presented in [12] showing that if DQs are not included, pointwise projection errors degrade with the number of snapshots. However, this degradation is hard to find in many practical instances, and, at the same time, while some authors report improvement in their numerical simulations if DQs are included in the data set [8], others find just the opposite [10,11]. The present paper does not pretend to settle the argument, but hopes to shed some light on why, in agreement with the numerical experience of some authors, it may be possible to obtain good results in practice without including DQs. In particular, we show that if a function has first derivatives (with respect to time) square-integrable in time, POD projection errors do not degrade with the number of snapshots. Furthermore, if second or higher-order derivatives are square-integrable in time, POD methods are convergent and, the higher the order of the derivatives that are square integrable, the closer to optimal the rate of convergence becomes. We do this by obtaining discrete versions of Agmon and interpolation inequalities in Sobolev spaces, which allow to bound the L∞norm of a function in terms of the L2norm and higher-order Sobolev’s norms. We remark that our results do not contradict previous results in the literature. On the one hand, the error bounds we prove are as close to optimal as the smoothness of the solution allows, but not optimal. On the other hand, concerning the counterexamples in [12], we notice that they require the snapshots to be multiples of the POD basis and, hence, an orthogonal set, a property that is lost if a function is smooth and the snapshots are function values taken at increasingly closer times as the number of snapshots increases, because, on the one hand, the distance between two orthogonal vectors is not smaller than their norms and, on the other hand, for a continuous function, the distance between two consecutive function values of finer and finer partitions tends to zero. In conclusion, the snapshots of the counterexamples in [12] cannot be function values of the same function at denser and denser partitions of the same time interval. The reference [15] is the first one showing convergence results for POD methods for parabolic problems. Starting with this paper one can find several references in the literature developing the error analysis of POD methods. We mention some references that are not intended to be a complete list. In [16] the same authors of [15] provide error estimates for POD methods for nonlinear parabolic systems arising in fluid dynamics. The authors of [3] show that comparing the POD approximation with the L2projection over the reduced order space instead of the H1 0()-projection (as in [15]) DQs are not required in the set of snapshots. However, in that case, the norm of the gradient of the L2-projection has to be bounded. This leads to non-optimal estimates for the truncation errors, see [12]. POD approximation errors in different norms and using different projections can also be found in [19]. The rest of the paper is as follows. In Sect. 2we present the main results of this paper and their application to obtain pointwise error bounds in POD methods based only on snapshots. Section3present some numerical experiments. In Sect.4, the auxiliary results needed to prove the main results in Sect. 2are stated and proved, and the final section contains the conclusions. 123 Journal of Scientific Computing (2025) 103 :24 Page 3 of 26 24 2 Main Results Let be a bounded region in Rdwith a smooth boundary or a Lipschitz polygonal boundary. We use standard notation, so that Hs() denotes the standard Sobolev space of index s. We also consider the space H1 0() of functions vanishing on the boundary and with square integrable first derivatives. Let T>0andletu:×[0,T]→Rbe a function such that its partial derivative with respect to its second argument, denoted by thenceforth, satisfies that u∈Hm(0,T,X),whereXdenotes either L2() or H1 0(),and uHm(0,T,X)=m  k=0 1 T2(m−j)T 0 ∂k tu(·,t)  2 Xdt1/2 (1) and ·Xdenotes the norm in X,thatisfX=(f,f)1/2if X=L2() or fX= (∇f,∇f)1/2if X=H1 0(). Typically, in many practical instances, uis a finite element approximation to the solution of some time-dependent partial differential equation. The factors 1/T2(m−j)in the definition of ·Hm(0,T,X)are introduced so that the constants in the estimates that follow are scale-invariant (see e.g., [2, Chapter 4]). In what follows, if a constant is denoted with a lower case letter, this means that it is scale-invariant. Given a positive integer M>0 and the corresponding time-step t=T/M, we consider the time levels tn=nt,n=0,...,M, and the function values or snapshots of u un=u(·,tn), n=0,...,M. We denote by σ1≥...≥σJthe singular values and by ϕ1,...,ϕJthe left singular vectors of the operator from RM+1(endorsed with the standard Euclidean norm) which maps every vector x=[x1,...,xM+1]Tto Kx =1 √M+1 M  n=0 xn+1un,(2) operator which, being of finite rank, is compact. The singular values satisfy that λk= σ2 k,fork=1,...,Jare the positive eigenvalues of the correlation matrix M= 1 M+1((ui,uj)X)0≤i,j≤M. The left singular vectors, known as the POD basis, are mutually orthogonal and with unit norm in X.For1≤r≤dlet us denote by Vr=span(ϕ1,...,ϕr), and by Pr X:X→Vrthe orthogonal projection of Xonto Vr. It is well-known that 1 M+1 M  n=0 un−Pr Xun  2 X= k>r σ2 k.(3) In the present paper, the rate of decay given by the square root of the right-hand side above will be referred to as optimal (see [8]). In many instances of interest in practice, at least when uis a finite element approximation to the solution of a dissipative partial differential equation (PDE), the singular values decay exponentially fast, so that for very small values of r(compared with the dimension of the finite element space) the right-hand side of (3) is smaller than the error of uwith respect the solution of the PDE it approximates. POD methods are Galerkin methods where the approximating spaces are the spaces Vr. In POD methods, one obtains values un r∈Vrintended to be approximations un r≈un, 123 24 Page 4 of 26 Journal of Scientific Computing (2025) 103 :24 Table 1 Maximum pointwise projection error max0≤n≤M un−Pr Xun Xfor different values of the number of snapshots M+1andubeing the quadratic FE approximation to the periodic solution of the Brusselator problem described in [7], when X=L2(),r=25 and Tis one period M128 256 512 1024 max. err. 1.12 ×10−41.12 ×10−41.14 ×10−51.15 ×10−5 n=0,...,M,and(3) is key in obtaining estimates for the errors un r−un. Yet, the fact that the left hand side (3) is an average allows only for averaged estimates of the errors un r−un, and one would be interested in pointwise estimates. In order to obtain them, some authors assume (without a result to support it) that the individual errors  un−Pr Xun Xdecay with r as the square root of the right-hand side of (3) (see e.g., [9, (2.9)]). Yet, in principle, the only rigorous pointwise estimate that one can deduce from (3)is  un−Pr Xun X≤√M+1 k>r σ2 k1 2 ,n=0,...,M.(4) This estimate is sharp as shown in examples provided in [12], where for at least one value of n, the equality is reached. Estimate (4) degrades with the number M+1 of snapshots. One would expect then in practice that if Mis increased while maintaining the value of r the individual errors  un−Pr Xun Xshould increase accordingly. Yet this is not necessary the case as results in Table 1show, with correspond to the snapshots of the quadratic FE approximation to the periodic solution of the Brusselator problem described in [7](seealso Sect.3below) in the case where X=L2(),r=25 and Tis one period. The results in Table 1are in agreement with the following result. Theorem 2.1 Let u be bounded in H1(0,T,X). Then, for the constant cAdefined in (40) below, the following bound holds: max 0≤n≤M un−Pr Xun X≤cA2T k>r σ2 k1 4 (I−Pr X)∂tu 1 2 L2(0,T,X)+2 k>r σ2 k1 2 . If u0+···+uM=0, then the second term on the right-hand side above can be omitted. Proof The result is a direct consequence of the case m=1 in estimate (73) in Theorem 4.12 below when applied to f=u−Pr Xuand the identity (3) above, as we now explain. The norm ·0in (73)isdefinedin(38–39). With this norm and applying (3), we have  u−Pr Xu 0=T M M  n=0 u−Pr Xu  2 01 2 =M+1 MT k>r σ2 k1 2 ,(5) and the proof is finished by noticing that (M+1)/M≤2, for M≥1.  Observe that in the pointwise estimate in Theorem 2.1 no factor depending on Mappears on the right-hand side, except, perhaps, the singular values, but if u∈H1(0,T,X)then, σk→σc,kas M→∞,whereσc,1≥σc,2≥..., are the singular values of the operator Kc:L2(0,T)→Xgiven by [16, Section 3.2] Kcf=1 Tt 0 f(t)u(t)dt. 123 Journal of Scientific Computing (2025) 103 :24 Page 5 of 26 24 Observe also that, at worst,  (I−Pr X)∂tu L2(0,T,X)≤ ∂tuL2(0,T,X). Thus, Theorem 2.1 guarantees that pointwise projection errors do not degrade as the number of snapshots increases. Furthermore, they decay with r, although like k>rσ2 k1/4, which is the square root of the optimal rate. Yet, the following result shows that pointwise projection errors decay with a rate as close to optimal as the smoothness of uallows. In the sequel we define rM=M+1 M(6) Theorem 2.2 Let u be bounded in Hm(0,T,X). Then, for the constants cAdefined in (48) and cmin Lemma 4.10 below, the following bound holds: max 0≤n≤M un−Pr Xun X≤√2cAc 1 2 mrMT k>r σ2 k1 2−1 4m (I−Pr X)∂tu  1 2m Hm−1(0,T,X) +rM k>r σ2 k1 2 . If u0+···+uM=0, then the second term on the right-hand side above can be omitted. Proof SimilarlytoTheorem2.1, the result is a direct consequence of estimate (73) in Theorem 4.12 below applied to f=u−Pr Xu, the identity (3) above and (5).  The above results account for projection errors, but in the analysis of POD methods the error between the POD approximation and the snapshots, un r−undepends also on the difference quotients Dun=un−un−1 t,n=1,...,M.(7) for which we have the following result. Theorem 2.3 In the conditions of Theorem 2.1 the following bound holds: t M  n=1 (I−Pr X)Dun  2 X1 2 ≤2cmrMT k>r σ2 km−1 2m (I−Pr X)∂tu 1 m Hm−1(0,T,X). Proof The proof is a direct consequence of estimate (72) in Theorem 4.12 below applied to f=u−Pr Xu, the identity (3)and(5).  Remark 2.4 An important case in practice in that of periodic orbits, since they are key elements in bifurcation diagrams of many dynamical systems associated to PDEs. In this case, the snapshot u0can be omitted in (2)and(3), and M+1 can be replaced by Min previous formulae (see also Sect.4.1) and, consequently, the ratio rMin (6) can be replaced by 1. Also, in the estimates in Theorems 2.1 and 2.2 one can replace the constant cmby 1, which gives much more favourable estimates (see values of constants cmin Table 7as well as  (I−Pr X)∂tu Hm−1(0,T,X)can be replaced by  (I−Pr X)∂m tu L2(0,T,X),thatis 123 24 Page 6 of 26 Journal of Scientific Computing (2025) 103 :24 max 0≤n≤M un−Pr Xun X≤√2cAT k>r σ2 k1 2−1 4m (I−Pr X)∂m tu  1 2m L2(0,T,X) + k>r σ2 k1 2 ,(8) t M  n=1 (I−Pr X)Dun  2 X1 2 ≤2T k>r σ2 k1 2−1 2m (I−Pr X)∂m tu  1 m L2(0,T,X).(9) As in Theorem 2.1, the second term on the right-hand side of (8) can be omitted if the snapshots have zero mean. The proof of these estimates follows by applying Theorem 4.11 to (I−Pr X)u. Some reduction in the size of the constants cmcan be obtained also when computing quasi periodic orbits in an invariant tori if one uses the techniques described in [17] to compute them, since although u(T)= u(0), one still has u(T)−u(0)1, and one can modify the technique so that this is also the case with first derivatives or higher order ones. Remark 2.5 In the analysis of POD methods, as it will be the case in next section, one usually has to estimate (I−Pr X)in a norm ·other than that of X.If·is the norm of a Hilbert space Wcontaining the snapshots u0,...,uM, then one can use [12, Lemma 2.2] to get 1 M+1 m  n=0 (I−Pr X)un  2= k>r σ2 k ϕk  2.(10) Applying this estimate instead of the identity (3) in the proofs of Theorems 2.2 and 2.3 one gets the following estimates max 0≤n≤M un−Pr Xun ≤√2cAc 1 2 mrMT k>r σ2 k ϕk  21 2−1 4m (I−Pr X)∂tu  1 2m Hm−1(0,T,W) +rM k>r σ2 k ϕk  21 2 .(11) t M  n=1I−Pr X)Dun  21 2 ≤2cmrMT k>r σ2 k ϕk  2m−1 2m (I−Pr X)∂tu  1 m Hm−1(0,T,W).(12) As in Remark 2.4,whenuis periodic in twith period T, then, cmcan be replaced by 1,  (I−Pr X)∂tu Hm−1(0,T,W)by  (I−Pr X)∂m tuL2(0,T,W)and rMcan also be replaced by 1. 2.1 Application to a POD Method We now apply the previous results to obtain pointwise error bounds for a POD method applied to the heat equation where no DQ were included in the data set to obtain the POD basis {ϕ1,...,ϕJ}. To discretize the time variable we consider two methods, the implicit Euler method and the second-order backward differentiation formula (BDF2). Analysis similar to 123 Journal of Scientific Computing (2025) 103 :24 Page 7 of 26 24 part of computations done in the present section has been done before (see e.g., [6, Lemma 3]) but it is carried out here because better error constants are obtained. Let us assume that uis a semi-discrete FE approximation to the solution of the heat equation with a forcing term fso that u(t)∈Vfor some FE space Vof piecewise polynomials over a triangulation and satisfies (∂tu,ϕ)+ν(∇u,∇ϕ) =(f,ϕ), ∀ϕ∈V,(13) where ν>0 is the thermal diffusivity and, in this section, (·,·)denotes the standard inner product of L2(). We consider the POD method (Dun r,ϕ)+ν(∇un r,∇ϕ) =(f,ϕ), ∀ϕ∈Vr,n=1,...,M,(14) with u0 r=Rru(·,0), (15) where Rrdenotes the Ritz projection, (∇Rru,∇ϕ) =(∇u,∇ϕ), ∀ϕ∈Vr. Instead of using the first-order convergent Euler method considered above we may use the BDF2, with gives the the following method (Dun r,ϕ)+ν(∇un r,∇ϕ) =(f,ϕ), ∀ϕ∈Vr,n=2,...,M,(16) where Dun r=3 2Dun r−1 2Dun−1 r,n=2,...,M,(17) and, for simplicity, un r=Pr Xun,n=0,1.(18) We now prove the following result. Theorem 2.6 Let X be H1 0and let p =1in the case of method (14–15)and p =2otherwise. Assume that u ∈Hp+1(0,T,L2)∩Hm(0,T,H1 0)for some m ≥2. Then, the following bound holds for 0≤n≤M:  un r−un ≤4CP√Tc mrMT k>r σ2 k1 2−1 2m (I−Pr X)∂tu  1 m Hm−1(0,T,H1 0) +√2CPcAc1/2 mrMT k>r σ2 k1 2−1 4m (I−Pr X)∂tu  1 2m Hm−1(0,T,H1 0) +rM k>r σ2 k1 2 +p+3(t)p√T  ∂p+1 tu  L2(0,T,L2).(19) Proof Since the POD basis {ϕ1,...,ϕJ}has been computed with respect to the inner product in H1 0,wehavePr X=Rr. We first analyze method (14–15). We notice that (DPr Xun,ϕ)+ν(∇Pr Xun,∇ϕ) =(f,ϕ)+(τn,ϕ), ∀ϕ∈Vr,(20) where τn=DPr Xun−∂tu(·,tn), n=1,...,M. 123 24 Page 8 of 26 Journal of Scientific Computing (2025) 103 :24 Subtracting (20) from (14), for the error en r=un r−Pr Xun,n=0,...,M, we have the following relation, (Den r,ϕ)+ν(∇en r,∇ϕ) =(τn,ϕ), ∀ϕ∈Vr,n=1,...,M.(21) Taking ϕ=ten rand using that (en r−en−1 r,en r)=( en r  2− en−1 r  2+ en r−en−1 r  2)/2, one obtains  en r  2− en−1 r  2+2νt ∇en r  2≤2t τn  en r ,(22) If  en r > en−1 r ,then−1<− en−1 r / en r . Hence, dividing both sides of (22)by en r  we get  en r − en−1 r ≤2t τn .(23) If, on the contrary,  en r ≤ en−1 r , for the right-hand side in (22) we write  τn  en r ≤  τn ( en r + en−1 r )/2 so that dividing both sides of (22)by en r + en−1 r ,wealso obtain (23). Summing from n=1ton=mand using (15) we obtain  em r ≤2t m  n=1 τn ≤2√Tt m  n=1 τn  21 2 ,(24) where, in the last step we have applied Hölder’s inequality. With respect to τn, by adding ±Dunwe have  τn ≤ (I−Pr X)Dun + Dun−∂tu(·,tn) . If u∈H2(0,T,L2)Taylor expansion with integral remainder allows us to write  Dun−∂tu(·,tn) ≤1 t   tn tn−1 (t−tn−1)∂2 tu(·,t)dt   ≤√t ∂2 tu L2(tn−1,tn,L2). Thus, from (24)weget  em r ≤2√Tt m  n=1 (I−Pr X)Dun  21 2 +2t√T ∂2 tu L2(0,T,L2). Using Poincaré inequality  v ≤CP ∇v ,∀v∈H1 0(), (25) and applying Theorem 2.3 we have 123 Journal of Scientific Computing (2025) 103 :24 Page 9 of 26 24 max 0≤n≤M en r ≤4cmCP√TrMT k>r σ2 km−1 2m (I−Pr X)∂tu  1 m Hm−1(0,T,H1 0) +2t√T ∂2 tu L2(0,T,L2). Finally, to estimate the errors un−un rwe write un r−un=en r+Pr Xun−un, and apply (25) and Theorem 2.2 to conclude (19). The analysis of method (16–18) is very similar to that of (14), and we comment on the differences next. As it is well-known (and a simple calculation shows) (Den r,en r)=E2 n−E2 n−1, where En=1 2 en r  2+ 2en r−en−1 r  21/2 ,n=1,...,M. It is easy to check that  en r ≤En, so that one can repeat (with obvious changes) the analysis above with the implicit Euler method with  en r replaced by Ento reach  em r ≤Em≤t m  n=2 Tn ≤√Tt m  n=2 Tn  21 2 ,(26) instead of (24),wherewehaveusedthatE1=0 due to (18), and where Tn=Pr XDun−∂tu(·,tn)=(Pr X−I)Dun+Dun−∂tu(·,tn). For the last term above we have (see e.g., [6, Lemma 4])  Dun−∂tu(·,tn) ≤√5(t)3/2tn tn−2 ∂3 tu(·,t)  2dt1 2 . Consequently, from (26) and noticing that  Den r ≤(3/2) Den r +(1/2) Den−1 r ,wehave  em r ≤2√Tt m  n=1 (I−Pr X)Dun  21/2 +√5T(t)2 ∂3 tu L2(0,T,L2),(27) so that applying (25)andTheorem2.3 it follows that max 0≤n≤M en r ≤4cmCP√TrMT k>r σ2 km−1 2m (I−Pr X)∂tu  1 m Hm−1(0,T,H1 0) +(t)2√5T ∂3 tu L2(0,T,L2), and, for the error un r−u=en r+(I−Pr X)un, applying Theorem 2.2 we get (19).  Remark 2.7 In view of Remark 2.5, in the above estimates we may replace the Poincaré constant CPby 1, σ2 kby σ2 k ϕk  2and the norm ·Hm−1(0,T,H1 0)by ·Hm−1(0,T,L2). Also in view of Remark 2.4,whenuis periodic in twith period T, the constant cmcan be replaced by 1 and  (I−Pr X)∂tu Hm−1(0,T,H1 0)by  (I−Pr X)∂m tu L2(0,T,H1 0)(or by  (I−Pr X)∂m tu L2(0,T,L2) if we use the estimates in Remark 2.5). Remark 2.8 If, for the method (16), u1 ris computed from u0 rwith one step of the implicit Euler method, instead of setting u1 r=Pr Xu1, then one can use a general stability result for the BDF2 like [6, Lemma 3] and then, use standard estimates to express the term Du1−∂tu(·,t1) in τ1as  Du1−∂tu(·,t1) ≤t ∂2 tu L∞(0,t1,L2).Thus, one gets estimates similar to (19) 123 24 Page 16 of 26 Journal of Scientific Computing (2025) 103 :24 A similar argument shows that (43) also holds for m<l≤M. Now, since we are assuming that fτis of zero mean, we can write fm=fm−1 M+1 M  l=0 fl=1 M+1 l=m (fm−fl), and, thus, taking norms, applying (43)  fm ≤√τ M+1 l=m|m−l|1/2 Dfτ 0. Now, a simple calculation shows that 1 M+1 l=m|m−l|1/2≤1 M+1 M  l=1 l1/2 and it is not difficult to prove that the last expresion does not exceed (2/3)√M.Thus,it follows that  fm ≤2 3√T Dfτ 0.(44) We now take msuch that  fm =min 0≤l≤M fl .(45) We have that  fm =√T √T fm =1 √TτM fm  21/2≤1 √T fτ 0.(46) From this inequality and (44) it follows that for fmsatisfying (45)wehave  fm  2≤2 3 fτ 0 Dfτ 0. This, together with (42) finishes the proof.  Remark 4.2 Notice that (41) implies (42)alsointhecasewherefτis defined as  fτ 0=τ 2 f0  2+τ M−1  n=1 fn  2+τ 2 fM  21/2 ,(47) and that (46) also holds in this case. Thus, Lemma 4.1 is valid if fτis defined as in (47). Remark 4.3 If the Hilbert space Xis either Rdfor some d≥0, or L2(),forsome⊂Rd or (a subespace of) the Sobolev space Hs() for some integer s>0, then the constant cA in (40) can be set to cA=√2.(48) We now show that this is the case for X=L2() since the argument is easily adapted to the other cases. We first notice that since the sequence fτhas zero mean, then, for every x∈ (except maybe in a set of zero measure) not all values (fn(x))M n=0can have the same sign. 123 Journal of Scientific Computing (2025) 103 :24 Page 17 of 26 24 Consequently, for every x∈there is an integer mx∈[1,M]and a value λx∈[0,1]for which 0 =(1−λx)fmx(x)+λxfmx−1(x)=0. We notice then that f2 mx(x)=fmx(x)fmx(x)=fmx(x)( fmx(x)−((1−λx)fmx(x)+λxfmx−1(x)) =fmx(x)λx(fmx(x)−fmx−1(x)). Thus, for x∈such that mx≤n, arguing as in the proof of Lemma 4.1 we may write f2 n(x)=f2 mx(x)+τ n  j=mx+1 Re (fj(x)+fj−1(x))Dfj(x) ≤τλxfmx(x)Dfmx(x)+τ n  j=mx+1fj(x)+fj−1(x)Dfj(x). This argument is easily adapted to the case where mx>n. Observe that the last inequality implies f2 n(x)≤τ M  j=1fj(x)+fj−1(x)Dfj(x), and a similar argument shows that this inequality to be valid also when mx>n. Thus, integrating in it easily follows that fn2≤2fτ0Dfτ0. We now treat first the periodic case, since it is less technical than the general one and the main ideas can be presented more easily. 4.1 The Periodic Case In this subsection we assume that fM=f0and we extend periodically the sequence fτ,that is fm=fM+mfor m=±1,±2,.... Also, we will consider fτ0as defined in (47). We define  Dkfτ 0as  Dkfτ 0=τ M  n=1 Dkfn  21/2 ,k=2,...,M.(49) Notice that for k=1 the expression above coincides with Dfτdefined in (39). Also, since f0=fM, the expression in (49)fork=0 coincides with the alternative definition of fτ given in (47), for which, as observed in Remark 4.2 Lemma 4.1 is also valid. The proof of the following lemma is the discrete counterpart of integration by parts in a periodic function  f  2=T 0 (f(t), f(t)) dt=−T 0 (f(t), f(t)) dt≤fL2(0,T,X) f L2(0,T,X). Lemma 4.4 Let fτ=(fn)n∈Zbe a periodic sequence in X with period M. Then, for k = 2,...,M−1,  Dkfτ 0≤ Dk−1fτ  1/2 0 Dk+1fτ  1/2 0 123 24 Page 18 of 26 Journal of Scientific Computing (2025) 103 :24 Proof Since τ Dkfn  2=τ(Dkfn,Dkfn)=(Dkfn,Dk−1fn−Dk−1fn−1),wehave  Dkfτ  2 0= M  n=1 (Dkfn,Dk−1fn−Dk−1fn−1) =(DkfM,Dk−1fM)−(Dkf1,Dk−1f0)− M−1  n=1 (Dkfn+1−Dkfn,Dk−1fn). Now we notice that for the second term on the right-hand side above we have (Dkf1,Dk−1f0)= (DkfM+1,Dk−1fM), and thus  Dkfτ  2=− M  n=1 (Dkfn+1−Dkfn,Dk−1fn)=−τ M  n=1 (Dk+1fn+1,Dk−1fn) and the proof is finished by applying Hölder’s inequality.  The proof of the following reuslt follows the ideas in the proof of [1, Lemma 4.12] adapted to the discrete case. Lemma 4.5 Let fτ=(fn)n∈Za periodic sequence in X with period M. Then, for m = 2,...,M,  D1fτ 0≤ fτ  m−1 m 0 Dmfτ  1 m 0.(50) Proof We will prove below that (50) will be a consequence of the following inequality  Dkfτ 0≤ fτ  1 k+1 0 Dk+1fτ  k k+1 0,k=1,...,M−1,(51) which we will now proof by induction. Notice that (51) holds for k=1 as a direct application of Lemma 4.4. Now assume that (51) holds for k=1,...,m−1, and let us check that it also holds for k=m. Applying Lemma 4.4 we have  Dmfτ 0≤ Dm−1fτ  1/2 0 Dm+1fτ  1/2 0. Applying the induction hypothesis we have  Dmfτ 0≤ fτ  1 2m 0 Dmfτ  m−1 2m 0 Dm+1fτ  1/2 0, so that dividing by  Dmfτ  m−1 2m 0we obtain  Dmfτ  m+1 2m 0≤ fτ  1 2m 0 Dm+1fτ  1/2 0, from where (51) follows. Now, to proof (50), we start with the case m=2, which holds as a direct application of Lemma 4.4 for k=1. Then we apply (51) repeatedly to the factor with the highest-order difference, that is  D1fτ 0≤ f  1 2 D2fτ  1 2 0≤ f  1 2+1 6 D3fτ  1 3 0≤ f  1 2+1 6+1 12  D4fτ  1 4 0 ≤...≤ f  1 2+1 6+··· 1 m(m−1) Dmfτ  1 m 0, 123 Journal of Scientific Computing (2025) 103 :24 Page 19 of 26 24 and the proof is finished by noticing that 1 2+1 6+··· 1 m(m−1)=1 2+1 2−1 3+1 3−1 4+···+1 m−1−1 m ≤1−1 m=m−1 m.  From Lemmas 4.1 and 4.5 the following result follows. Lemma 4.6 Let fτ=(fn)n∈Za periodic sequence in X with period M and zero mean. Then, for k =1,...,M, the following estimate holds:  fn ≤cA fτ  1−1 2m 0 Dmfτ  1 2m 0,(52) where cAis the constant in (40). 4.2 The General Case In this section, and unless stated otherwise, fτ=(fn)M n=0is a sequence in Xwith zero mean. We define  Dkfτ 0=τ M  n=k Dkfn  21/2 ,k=1,...,M, and, in order to obtain scale invariant constants, for each k=0,...,M,  Dkfτ m=k+m  j=k 1 T2(m+k−j) j Djfτ  2 01/2 ,m=1,...,M−k,(53) where, here and in the sequel, T0=T,and Tk=(M+1−k)τ, k=1,...,M.(54) We also denote mk=1 M+1−k M  n=k Dkfn.(55) Observe that by applying Hölder’s inequality one gets  mk ≤1 √Tk Dkfτ 0(56) On the other hand, it is easy to check mk=minα∈R Dkfτ−α 0, where by Dkfτ−α we denote the sequence (Dkfn−α)M n=k. Consequently, for the sequence Dkfτ−mk= (Dkfn−mk)M n=k,wehave  Dkfτ−mk 0≤ Dkfτ 0.(57) Our first result extends Lemma 4.1 to the sequences (Dkfn)M n=k, whose means are not 0. 123 24 Page 20 of 26 Journal of Scientific Computing (2025) 103 :24 Lemma 4.7 The following bound holds for k =1,...M,  Dkfn ≤cA,1 Dkfτ  1/2 0 Dkfτ  1/2 1,n=k,...,M−1, where cA,1=1+c4/3 A3/4 ,(58) and cAis the constant in (40). Proof Applying Lemma 4.1 we have  Dkfn−mk ≤cA Dkfτ−mk  1/2 0 Dk+1fτ  1/2 0 ≤cA Dkfτ  1/2 0 Dk+1fτ  1/2 0,(59) where in the last inequality we have applied (57). Thus, by writing Dkfn=Dkfn−mk+mk andinviewof(56), we have  Dkfn ≤cA Dkfτ  1/2 0 Dk+1fτ  1/2 0+T−1/2 k Dkfτ 0 = Dkfτ  1/2 0cA Dk+1fτ  1/2 0+T−1/2 k Dkfτ  1/2 0. Now applying Hölder’s inequality, ab +cd ≤(ap+cp) 1 p(bq+dq) 1 q, to the second factor above, with p=4/3andq=4, the proof is finished.  Next,weextendLemma4.4 to the general case Lemma 4.8 For k =1,...,M−1,  Dkfτ 0≤cB,1 Dk−1fτ  1/2 0 Dkfτ  1/2 1, where cB,1=2+4(cAcA,1)21/2.(60) Proof We notice that Dkfτ=D(Dk−1fτ−mk−1)and let us denote yn=Dk−1fn−mk−1,n=k−1,...,M. Thus, arguing as in the proof of Lemma 4.4 we have  Dkfτ  2= M  n=k (Dkfn,yn−yn−1) =(DkfM,yM)−(Dkfk,yk−1)− M−1  n=k (Dkfn+1−Dkfn,yn) =(DkfM,yM)−(Dkfk,yk−1)−τ M−1  n=k (Dk+1fn+1,yn). We apply (59) and Lemma 4.7 to the first two terms on right-hand side above, and we apply (57) to the last one, to get,  Dkfτ  2≤2cAcA,1 Dk−1fτ  1 2 0 Dkfτ 0 Dkfτ  1 2 1+ Dk−1fτ 0 Dk+1fτ 0 ≤1 2 Dkfτ  2 0+2c2 Ac2 A,1 Dk−1fτ 0 Dkfτ 1+ Dk−1fτ 0 Dk+1fτ 0. 123 Journal of Scientific Computing (2025) 103 :24 Page 21 of 26 24 Thus,  Dkfτ  2≤ Dk−1fτ 04c2 Ac2 A,1 Dkfτ 1+2 Dk+1fτ 0 ≤ Dk−1fτ 04c2 Ac2 A,1+2 Dkfτ 1, and the proof is finished for k≥2.  Lemma 4.9 For m =2,...,M,  D1fτ 0≤cm fτ  m−1 m 0 D1fτ  1 m m−1,(61) where cm= m−2  j=0ˆc 1 j+1 j.(62) and the constants ˆcjare given in (66)below. Proof The result will be a consequence of the following inequality that will be proven below by induction:  Dkfτ j≤ˆcj Dk−1fτ  1 j+2 0 Dkfτ  j+1 j+2 j+1.(63) If (63) holds for j=0,...,m−2, then, applying it successively we have  D1fτ 0≤ˆc0 fτ  1 2 0 D1fτ  1 2 1≤ˆc0ˆc 1 2 1 fτ  1 2+1 6 0 D1fτ  1 3 2≤... ≤cm fτ  β 0 D1fτ  1 m m−1 where β= m−1  j=1 1 j(j+1)= m−1  j=11 j−1 j+1=1−1 m=m−1 m,(64) so that (61) follows. We now prove (63). Lemma 4.8 shows it holds for j=0 with ˆc0=cB,1. Let us assume that it holds for j=0,...,m−2, and let us show that it holds for m−1. We do this for k=1 since the argument is the same for k>1. We notice  D1fτ m−1=T−2(m−1) D1fτ  2 0+ D2fτ  2 m−21 2 and apply the induction hypotheses to the second term on the right-hand side above to get  D1fτ m−1=T−2(m−1) D1fτ  2 0+ˆc2 m−2 D1fτ  2 m 0 D2fτ  2(m−1) m m−11 2 = D1fτ  1 m 0T−2(m−1) D1fτ  2(m−1) m 0+ˆc2 m−2 D2fτ  2(m−1) m m−11 2 =AB. Now, we use that, as argued above, the induction hypothesis implies that (61) holds, so that we can write  D1fτ m−1=c 1 m m fτ  m−1 m2 0 D1fτ  1 m2 m−1B, 123 24 Page 22 of 26 Journal of Scientific Computing (2025) 103 :24 so that, dividing by  D1fτ  1 m2 m−1we have  D1fτ  m2−1 m2 m−1=c 1 m m fτ  m−1 m2 0B, which implies  D1fτ m−1=c m m2−1 m fτ  1 m+1 0B m2 m2−1.(65) Now we apply Holder’s inequality ab +cd ≤(ap+cp)1/p(bq+dq)1/qto B, with p=m and q=m/(m−1) B≤1+ˆc2m m−21 2mT−2m D1fτ  2 0+ D2fτ  2 m−1m−1 2m =1+ˆc2m m−21 2m D1fτ  m−1 m m, which, together with (65), implies  D1fτ m−1=1+ˆc2m m−2c2 mm 2(m+1)(m−1) fτ  1 m+1 0 D1fτ  m m+1 m. This is (63)for j=m−1 with ˆcj=1+ˆc2(j+1) j−1j+1 2(j+2)jj−1  i=0ˆc 1 i+1 ij+1 (j+2)j ,ˆc0=cB,1,(66) where cB,1is the constant given in (60).  Table 7shows the first 5 values of cmfor cAgiven by (48). Computation of further values suggests that cm≈γm−1with γ=5.9416. From Lemmas 4.1 and 4.9 the following result follows. Lemma 4.10 Let fτ=(fn)M n=0a sequence in X with zero mean. Then, for m =1,...,M−1, the following estimate holds:  fn ≤cAc 1 2 m fτ  1−1 2m 0 D1fτ  1 2m m−1,(67) where cAis the constant in (40)(or (48)for X being a space of functions or RN) and cmis cm=1if m =1and the constant in (62)for m ≥2. 4.3 Application to Sequences of Function Values We consider now the case where Xis either L2() or some Sobolev space Hs() for some domain ⊂Rd,·its norm, and fn=f(·,tn)for some function f:×[0,T]→R such that f∈Hm(0,T,X). We first notice that since for x∈we have Dk−1f(x,tn)= ∂k−1f(x,ξ)for some ξ∈(tn−k+1,tn), then it follows that  Dkfn ≤√k √τ ∂kf L2(tn−k,tn,X),k≤n≤M.(68) Consequently,  Dkfτ 0≤k ∂kf L2(0,T,X).(69) We have the following result. 123 Journal of Scientific Computing (2025) 103 :24 Page 23 of 26 24 Theorem 4.11 Let f :X×R→Ra function which is T -periodic in its last variable and such that f ∈Hm(0,T,X). Then, the sequence fτ=(fn)M n=1satisfies the following bounds for 1≤m≤M−1:  Dfτ 0≤2m−1 mm1 mfτ m−1 m 0 ∂mf  1 m L2(0,T,X),  fn ≤2m−1 2mm1 2mcAfτ1−1 2m 0 ∂mf  1 2m L2(0,T,X)+1 √T fτ 0,1≤n≤M, where cmis the constant in Lemma 4.10.If f 1+···+ fM=0, the the last term on the right-hand side above can be omitted. Proof Let us denote m0=(f1+···+fM)/M. Applying Hölder’s inequality is easy to check that |m0|≤T−1/2 fτ 0and, consequently, fτ−m00≤2fτ0. Thus, by expressing fn=(fn−m0)+m0, applying Lemmas 4.5 and 4.6 to fτ−m0and using (69)theresult follows for m≥2. For m=1 the result follows form Lemma 4.6 and (69).  For the general case, in view of the definition of  Dkfτ min (53)and(69)wehave  Dfτ  1 m m−1≤m−1  k=0 (k+1)2T Tk+12(m−k)1 T2(m−k) ∂k+1 tf  2 L2(0,T,X)1 2m ≤max 0≤k≤m−1(k+1)1 mT Tk+1m−k m ∂tf  1 m Hm−1(0,T,X),(70) where ·Hm−1(0,T,X)is defined in (1). We now estimate the maximum above. In view of the expression of Tkin (54), for 1 ≤k≤m−1wehave T Tk+1m−k m =M M−km−k m =1+k M−km−k m =1+1 (M/k)−1m−k m =1+1 (M/k)−1(M/k)−1m−k m((M/k)−1) ≤e m−k m((M/k)−1)=e k(m−k) m(M−k). Now the function x→ x(m−x)/(m(M−x)) is monotone increasing in the interval [1,m−1], unless m≤M/4 where it has a relative maximum at (M/2)−(M2/4)−Mm, but is easy to check that this value is strictly larger than m. Consequently, its maximum is achieved at x=m−1. Noticing that x1/x≤e1 efor x>0, we conclude that, for 2 ≤m≤M max 0≤k≤m−1(k+1)1 mT Tk+1m−k m ≤em−1 m(M−m+1)+1 e≤e1+1 e. One can check that e1+1 e≤3.93. Thus, from (70) it follows that  D1fτ  1 m m−1≤4 ∂tf  1 m Hm−1(0,T,X).(71) 123 24 Page 24 of 26 Journal of Scientific Computing (2025) 103 :24 Theorem 4.12 Let f ∈Hm(0,T,X). Then, the sequence fτ=(fn)M n=0satisfies the following bound for 1≤m≤M:  D1fτ ≤2m−1 mcmfτ m−1 m 0∂tf 1 m Hm−1(0,T,X),(72)  fn ≤2m−1 2mcAc1/2 m fτ  1−1 2m 0∂tf 1 2m Hm−1(0,T,X)+1 √T fτ 0, n=0,...,M.(73) If f0+···+ fM=0, the last term on the right-hand side above can be omitted. Proof Let us denote m0=(f0+···+fM)/(M+1). Applying Hölders inequality one sees that |m0|≤T−1/2 fτ 0so that fτ−m00≤2fτ0. Thus, by writing fn=(fn−m0)+m0, applying Lemmas 4.9 (if m≥2) and 4.10 to fτ−m0and using (71) the proof is finished.  5 Summary and Conclusions Several estimates, pointwise and averaged, have been proved for POD methods whose data set does not include DQs of the snapshots. These estimates include pointwise projection errors of snapshots (Theorems 2.1 and 2.2 in the general case and estimate (8) in the case of periodic functions in time) for orthogonal projections in L2and H1 0. They also include L2norm of H1 0-projection errors (estimate (11)), for which estimate (10) has been used. Our estimates also include (the square root of) the average of squared projection error of DQs (Theorem 2.3 for the general case and estimate (9) for periodic functions) again in the L2and H1 0cases, as well as the average of squares of L2norms of Ritz projection errors of DQs (estimate (12)). In all cases, when projections are onto the space spanned by the first relements of the POD basis, the rate of decay in terms of σ2 r+1+···+σ2 Jis as close to the optimal value (σ2 r+1+···+σ2 J)1/2as the smoothness of the function from where the snapshots are taken allows. These estimates have been used to obtain error bounds with the same rates of decay for POD methods applied to the heat equation, where both backward Euler method and the BDF2 are considered for discretization of the time variable (Theorem 2.6). The heat equation has been chosen for simplicity, but it is clear that, with techniques for nonlinear problems such as those, for example, in [5,7,9,13,16], it is possible to extend these estimates to POD methods not using DQs in the data set for nonlinear equations. Numerical experiments are also presented in this paper, where snapshots are taken from two challenging problems due to the size of derivatives with respect to time (see Fig. 1)and Table 4), derivatives which feature in the estimates commented above. In these experiments, pointwise projection errors are overestimated by our error bounds by factors not exceeding 28, and by a factor also not exceeding 28 in the case of the (square root of) the average of squares of L2norms of Ritz projection errors of DQs, if the estimates are those in Remark 2.7, although these factors are likely to be problem-dependent. This research was motivated by a better understanding of POD methods that do not include DQs of the snapshots in the data set. In particular by the apparent contradiction that while counterexamples in [12] clearly show that pointwise projection errors do degrade with the number of snapshots, computations like those in Table 1suggest that this may not necessarily be the case in practice. By taking into account the smoothness of functions from where the snapshots are taken, new estimates have been obtained and the apparent contradiction men123 Journal of Scientific Computing (2025) 103 :24 Page 25 of 26 24 tioned above has been explained, as well as better knowledge of the decay rate of pointwise errors when DQs are omitted from the data set. Funding Funding for open access publishing: Universidad de Sevilla/CBUA Bosco García-Archilla: Research is supported by grants PID2021-123200NB-I00 and PID2022-136550NB-I00 funded by MCIN/AEI/10.13039/ 501100011033 and by ERDF A way of making Europe, by the European Union. Julia Novo: Research is supported by grant PID2022-136550NB-I00 funded by MCINAEI/10.13039/501100011033 and by ERDF A way of making Europe, by the European Union. Data Availability The authors declare that the data supporting the findings of this study are available are availabe at https://doi.org/10.12795/11441/169465 and https://hdl.handle.net/11441/169465. Source code and part of data are also available at https://github.com/bgarchilla/pointwise/. Declarations Conflict of interest The authors declare that they have no Conflict of interest. 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. Adams, R. A.: Sobolev spaces, Academic Press [A subsidiary of Harcourt Brace Jovanovich, Publishers], New York-London, (1975) 2. Constantin, P., Foias, C.: Navier-Stokes Equations. Chicago Lectures on Mathematics. The University of Chicago Press, Chicago (1988) 3. Chapelle, D., Gariah, A., Sainte-Marie, J.: Galerkin approximation with proper orthogonal decomposition: new error estimates and illustrative examples. ESAIM: M2AN, 46 731-757 (2012) 4. Eskew, S.L., Singler, J.R.: A new approach to proper orthogonal decomposition with difference quotients. Adv. Comput. Math. 49(2), 13 (2023) 5. García-Archilla, B., John, V., Novo, J.: POD-ROMs for incompressible flows including snapshots of the temporal derivative of the full order solution. SIAM. J. Numer. Anal. 61(3), 1340–1368 (2023) 6. García-Archilla, B., John, V., Novo, J.: Second order error bounds for POD-ROM methods based on first order divided differences. Appl. Math. Letters 143, 108836 (2023) 7. García-Archilla, B., Novo, J.: POD-ROM methods: from a finite set of snapshots to continuous-intime approximations. SIAM J. Numer. Anal. (to apper) https://doi.org/24M1645681.arXiv:2403.06967 [math.NA], 8. Iliescu, T., Wang, Z.: Are the snapshot difference quotients needed in the proper orthogonal decomposition? SIAM J. Sci. Comput. 36, A1221–A1250 (2014). https://doi.org/10.1137/130925141 9. Iliescu, T., Wang, Z.: Variational multiscale proper orthogonal decomposition: Navier-Stokes equations. Numer. Meth. PDEs 30(2), 641–663 (2014) 10. John, V., Moreau, B., Novo, J.: Error analysis of a SUPG-stabilized POD-ROM method for convectiondiffusion-reaction equations. Comput. Math. Appl. 122, 48–60 (2022). https://doi.org/10.1016/j.camwa. 2022.07.017 11. Kean, K., Schneier, M.: Error analysis of supremizer pressure recovery for POD based reduced-order models of the time-dependent Navier-Stokes equations. SIAM J. Numer. Anal. 58, 2235–2264 (2020). https://doi.org/10.1137/19M128702X 12. Koc, B., Rubino, S., Schneier, M., Singler, J., Iliescu, T.: On optimal pointwise in time error bounds and difference quotients for the proper orthogonal decomposition. SIAM J. Numer. Anal. 59(4), 2163–2196 (2021) 123