Recursive POD expansion for the advection-diffusion-reaction equation

Commun. Comput. Phys. doi: 10.4208/cicp.OA-2017-0257 Vol. 24, No. 5, pp. 1556-1578 November 2018 Recursive POD Expansion for the Advection-Diffusion-Reaction Equation M. Aza¨ıez1,∗, T. Chac´on Rebollo2, E. Perracchione3and J. M. Vega4 1Bordeaux INP, I2M (UMR CNRS 5295), Universit´e de Bordeaux, 33607 Pessac, France. 2Departamento EDAN & IMUS, Universidad de Sevilla, C/Tarfia, s/n, 41012 Sevilla, Spain. 3Dipartimento di Matematica Tullio Levi-Civita, Universit`a di Padova, Padova, Italy. 4E.T.S.I. Aeron´autica y del Espacio, Universidad Polit´ecnica de Madrid, Madrid, Spain. Received 19 December 2017; Accepted (in revised version) 4 April 2018 Abstract. This paper deals with the approximation of advection-diffusion-reaction equation solution by reduced order methods. We use the Recursive POD approximation for multivariate functions introduced in [5] and applied to the low tensor representation of the solution of the reaction-diffusion partial differential equation. In this contribution we extend the Recursive POD approximation for multivariate functions with an arbitrary number of parameters, for which we prove general error estimates. The method is used to approximate the solutions of the advection-diffusion-reaction equation. We prove spectral error estimates, in which the spectral convergence rate depends only on the diffusion interval, while the error estimates are affected by a factor that grows exponentially with the advection velocity, and are independent of the reaction rate if this lives in a bounded set. These error estimates are based upon the analyticity of the solution of these equations as a function of the parameters (advection velocity, diffusion, reaction rate). We present several numerical tests, strongly consistent with the theoretical error estimates. AMS subject classifications: 65D15, 65M70, 74F10 Key words: Reduced order methods, recursive POD, multivariate functions. 1 Introduction Physical and real phenomena are often described by large scale problems for which a low dimensional approximation turns out to be essential. Therefore, in the last decades, ∗Corresponding author. Email addresses: [email protected] (M. Aza¨ıez), [email protected] (T. Chac´on), [email protected] (E. Perracchione), josemanuel.veg[email protected]s (J. M. Vega) http://www.global-sci.com/ 1556 c 2018 Global-Science Press M. Aza¨ıez et al. / Commun. Comput. Phys., 24 (2018), pp. 1556-1578 1557 ROMs, known as Karhunen-Lo`eve expansion (KLE), [21], have gained popularity. Such techniques are also known as Proper Orthogonal Decomposition (POD), [8], in mechanical computation, Principal Component Analysis (PCA) in statistics, [7, 18], and Singular Value Decomposition (SVD) in linear algebra, [16]. The POD is a very general information compression technique which finds its natural applications in a wide variety of fields, such as in digital image compression, bioinformatics, signal processing and resolution of PDEs, [3,4,6,8,19,25] and also the recent work of Cohen et al related to the numerical analysis of the convergence of polynomial approximation of parametric elliptic PDE’s [11,12]. Hence, the problem of the tensor representation of multivariate phenomena, by means of the KLE, is a challenging computational issue. Recently, techniques as the Proper Generalized Decomposition (PGD), [17], which gives a suitable approximation of multivariate functions by low dimensional varieties, or as the Higher Order SVD (HOSVD), [13, 14], have been developed. The latter provides a multilinear generalization of the best rank-lapproximation problem for matrices obtained by truncating the SVD. In fact, the HOSVD enables us to perform a low dimensional approximation of tensors as the SVD allows to approximate bivariate data. Also, a recent variant of the HOSVD is the Truncated Sequential HOSVD (ST-HOSVD), introduced in [26], that consists in sequentially performing a SVD analysis of the successive modes. The fact that each step for the construction of the ST-HOSVD separates two variables by the SVD procedure allows to exactly compute the truncation error in terms of the remaining singular values. This paper focuses on the study of the low tensor representation of the solution of the reaction-diffusion partial differential equation. The Recursive POD (RPOD) was proposed in [5], as a procedure to build low dimensional representation of trivariate functions. It was applied to the diffusion-reaction equation, for which an exponential rate of convergence, depending only on the diffusion rate, was proved. In the present paper we study the extension of the RPOD to approximate the solution of the advectiondiffusion-reaction equations. To do this, we start first by proving a more general result of convergence of the RPOD expansion to approximate parametric analytic function in L2 space. Then the analysis of the RPOD expansion of solution of reaction-diffusion partial differential equation turns to be a particular case where the parameters are advection velocity –assumed constant–, diffusion, reaction rate and space-time coordinates. We prove a spectral rate of convergence, that depends only on the diffusion rate much as in [5] for the convection-reaction equation. However here this rate is affected by a factor that exponentially grows with the P´eclet number, and is independent on the reaction rate if this lives in a bounded set. We also present some additional numerical tests that show the suitability of the RPOD to approximate analytic functions depending on a moderate number of parameters, with an exponential convergence rate. The outline of the paper is as follows. Section 2 focuses on the presentation of the RPOD expansion to multivariate functions with an arbitrary number of parameters. In section 3 we study the application of the RPOD expansion to approximate parametric 1558 M. Aza¨ıez et al. / Commun. Comput. Phys., 24 (2018), pp. 1556-1578 analytic functions valued in a L2space. We study the application of the RPOD to the approximation of the solution to the advection-reaction-diffusion equation in section 4. The paper ends with section 6 where we test the convergence rates for the solution of the advection-diffusion-reaction equations. Notation. Let X⊂Rdbe a given Lipschitz domain. We denote by L2(X)the space of measurable and square integrable functions on X. The scale of fractional Sobolev spaces Hτ(X),τ>0 is defined as in [1]. Finally, given the sets Xj⊂Rdj,j=1,···,nwe denote Xn i=Xi×Xi+1×···×Xnand xn i=(xi,xi+1,···,xn)∈Xi,i=1,···,n. 2 Low dimensional tensor representation via RPOD This section focuses on the presentation of the RPOD expansion to multivariate functions with an arbitrary number of parameters. 2.1 RPOD representation of multivariate functions In order to illustrate the RPOD representation of n-variate functions, at first we briefly review the basic concepts for bivariate functions, investigated in [5]. Let X1⊂Rd1and X2⊂ Rd2two bounded domains. The goal consists in approximating, in a low-dimensional variety, a function fin the Lebesgue space L2(X1×X2). Given a function f∈L2(X1×X2)≡L2(X1,L2(X2)), we define the integral operator Y:L2(X1)→L2(X2), with kernel f, as: (Yϕ)(x2)=ZX1 f(x1,x2)ϕ(x1)dx1, where ϕ∈L2(X1). The bounded operator Yhas an adjoint Y∗:L2(X2)→L2(X1), defined as: (Y∗v)(x1)=ZX2 f(x1,x2)v(x2)dx2, with v∈L2(X2). Once we define the integral operator Yand its adjoint, we consider the POD integral operator defined as R=Y∗Y. It follows that Ris a linear, bounded, nonnegative, self-adjoint and compact operator. The latter property is a consequence of the Kolmogorov compactness criterion in L2(X1)which states that Y∗is compact. Then, since Yis bounded, it implies that Ris compact as well [23]. Taking into account the properties of the operator R, the following spectral decomposition for the POD operator holds. Property 2.1. There exists a complete orthonormal basis of L2(X1)consisting of eigenvectors (ϕi)i∈Nand associated non-negative real eigenvalues (λi)i∈Nof Rthat we assume ordered in decreasing value. Each nonzero eigenvalue of Rhas finite multiplicity and 0 is the only possible accumulation point of the spectrum of R. Moreover, limi→∞λi=0. M. Aza¨ıez et al. / Commun. Comput. Phys., 24 (2018), pp. 1556-1578 1559 Concerning the computation of such basis, the sequence (vi)i∈N, given by: vi=1 σi Yϕi, (2.1) where σi=√λiare the so-called singular values and ϕi=1 σiY∗vi, is an orthogonal basis of L2(X2). Since Ris a positive operator, by virtue of the Mercer’s theorem [22], the following decomposition holds: Corollary 2.1. f(x1,x2)= ∑ i∈N σiϕi(x1)vi(x2). (2.2) The expansion (2.2) is convergent in L2(X1×X2). As already mentioned, the POD decomposition achieves the best approximation property, stated below. Property 2.2. Let Vl=span(v1,···,vl)⊂L2(X2)and let Wlbe any subspace of L2(X2)of dimension l. Then: ZX1 dL2(X2)(f(x1,·),Vl)2dx1≤ZX1 dL2(X2)(f(x1,·),Wl)2dx1, where dL2(X2)(ϕ,Wl)= inf ψ∈Wl||ϕ−ψ||L2(X2),ϕ∈L2(X2) denotes the distance from the element ϕ∈L2(X2)to the subspace Wl. In order to formulate a recursive POD representation for multivariate functions, let f be in the Lebesgue space L2(X1×X2×···×Xn), where X1⊂Rd1,X2⊂Rd2,···,Xn⊂Rdn are bounded domains. From Corollary 2.1 and taking into account that L2(Xn 1)and L2(Xn 2,L2(X1)) are isometric, we know that there exist two orthonormal sets (ϕi1)i1∈N and (vi1)i1∈N, which are respectively complete in L2(X1)and L2(Xn 2), such that: f(x1,x2,···,xn)= ∑ i1∈N σi1vi1(x2,···,xn)ϕi1(x1), (2.3) where the sum is convergent in L2(Xn 2,L2(X1)). The singular values are convergent to zero and non-negative. Now, we continue recursively applying the POD expansion to each mode vi1(x2,···,xn). Thus, since L2(Xn 2)and L2(Xn 3,L2(X2)) are isometric, there exist two orthonormal sets (ϕ(i1) i2)i2∈Nand (v(i1) i2)i2∈N, which are respectively complete in L2(X2)and L2(Xn 3), such that: vi1(x2,···,xn)= ∑ i2∈N σ(i1) i2v(i1) i2(x3,···,xn)ϕ(i1) i2(x2), 1560 M. Aza¨ıez et al. / Commun. Comput. Phys., 24 (2018), pp. 1556-1578 converges in L2(Xn 3,L2(X2)). Continuing in this fashion way we find the POD expansion of v(i1) i2,···,v(i1,i2,···,in−3) in−2, where the latter is found by considering two orthonormal sets (ϕ(i1,i2,···,in−2) in−1)in−1∈Nand (v(i1,i2,···,in−2) in−1)in−1∈N, respectively complete in L2(Xn−1)and L2(Xn)and such that: v(i1,i2,···,in−3) in−2(xn−1,xn)= ∑ in−1∈N σ(i1,i2,···,in−2) in−1v(i1,i2,···,in−2) in−1(xn)ϕ(i1,i2,···,in−2) in−1(xn−1), where the sum is convergent in L2(Xn,L2(Xn−1)) and the singular values (σ(i1,i2,···,in−2) in−1)in−1∈Nconverge to zero and are non-negative. Lemma 2.1. The function f ∈L2(Xn 1)admits the expansion: f=∑ i1∈N ∑ i2∈N··· ∑ in−1∈N σi1σ(i1) i2···σ(i1,i2,···,in−2) in−1ϕi1⊗ϕ(i1) i2⊗···⊗ϕ(i1,i2,···,in−2) in−1⊗v(i1,i2,···,in−2) in−1 =∑ i1∈N σi1∑ i2∈N σ(i1) i2···∑ in−1∈N σ(i1,i2,···,in−2) in−1ϕi1⊗ϕ(i1) i2⊗···⊗ϕ(i1,i2,···,in−2) in−1 ⊗v(i1,i2,···,in−2) in−1, (2.4) where sums are convergent in L2(Xn 1). Proof. It is enough to prove that the series are absolutely convergent. Thus we consider the following partial absolute sum: SN= N ∑ i1=0 N ∑ i2=0··· N ∑ in−1=0σi1σ(i1) i2···σ(i1,i2,···,in−2) in−1ϕi1⊗ϕ(i1) i2⊗··· ⊗ϕ(i1,i2,···,in−2) in−1⊗v(i1,i2,···,in−2) in−1 2 L2(Xn 1). Taking into account that eigenfunctions are orthonormal we have: SN= N ∑ i1=0 N ∑ i2=0··· N ∑ in−1=0σi1 2σ(i1) i2 2···σ(i1,i2,···,in−2) in−1 2 ≤ N ∑ i1=0 ∑ i2∈N··· ∑ in−1∈Nσi1 2σ(i1) i2 2···σ(i1,i2,···,in−2) in−1 2≤ N ∑ i1=0σi1 2≤f 2 L2(Xn 1), where the last inequality holds, since: ∑ i2∈Nσ(i1) i2 2=···=∑ in−1∈Nσ(i1,i2,···,in−2) in−1 2=1. (2.5) Then indeed the series are absolutely convergent. M. Aza¨ıez et al. / Commun. Comput. Phys., 24 (2018), pp. 1556-1578 1561 A similar result to the one outlined in Property 2.2 holds also for the RPOD: Property 2.3. Let f∈L2(Xn 1), then: f−fI 2 L2(Xn 1)≤f−˜ fI 2 L2(Xn 1), (2.6) where fI= I ∑ i=0 σivi(xn 2)ϕi(x1), (2.7) and ˜ fIis any n-variate approximation of fwith Imodes of the form: ˜ fI= I ∑ i=0 ˜ X(1) i(x1)˜ X(2) i(x2)··· ˜ X(n) i(xn). (2.8) Proof. We start pointing out that fIis the orthogonal projection in L2(Xn 2)of fon VI, where VI=span(ϕ1,···,ϕI). Let WIbe any other subspace of dimension Iof L2(Xn 2), then from Property 2.2 it follows: ZX1(f−fI)(x1) 2 L2(Xn 2)dx1≤ZX1(f−ˆ fI)(x1) 2 L2(Xn 2)dx1, for any ˆ fIin WI. Moreover since L2(Xn 1)and L2(X1,L2(Xn 2)) are isometric, (2.6) follows taking ˆ fI=˜ fI. Remark 2.1. The algebraic version of the RPOD is a particular case of the so called hierarchical tensor representation (see [28]). The link between the two approaches can be established in the same way as for the POD and SVD. See Volkwein [27] for the mathematical analysis. 3 RPOD applied to parametric analytic functions In this section we study the application of the RPOD expansion to approximate parametric analytic functions valued in a L2space. We intend to apply the results of this section to the approximation of the solution of the advection-diffusion-reaction equation, among other equations. Let us consider an open set Q ⊂ Rm1, a tensor parameter set G=∏n i=2Gi⊂Rm2×···×Rmn−1×Cand an analytic function f:G 7→ L2(Q). We shall denote by x1the elements of Qand by xithe elements of Gi,i=2,3,···,n. We further assume that all the Giare compact, and that Gncontains some ellipse Eρ=ζ∈C:| |ζ−1|+|ζ+1|≤ρ+ρ−1 1562 M. Aza¨ıez et al. / Commun. Comput. Phys., 24 (2018), pp. 1556-1578 contained in the half right complex plane, that is Eρ⊂C+={ζ∈C|Reζ>0}. Our purpose is to build an approximation of f, named fP, by truncating the POD expansion, defined as: fP= M1 ∑ i1=0 σi1 M2 ∑ i2=0 σ(i1) i2··· Mn−1 ∑ in−1=0 σ(i1,i2,···,in−2) in−1ϕi1⊗ϕ(i1) i2⊗···⊗ϕ(i1,i2,···,in−2) in−1⊗v(i1,i2,···,in−2) in−1. (3.1) The notation Pis a multi-index notation and stands for P=M1,M2,···,Mn−1. The estimation of rate of convergence of the truncated expansion fPtowards fis stated as follows Theorem 3.1. The truncated RPOD series expansion fPgiven by (3.1)satisfies the error estimate kf−fPkL2(G×Q)≤CρKρ−M1+pM1ρ−M2+···+pM1M2···Mn−1ρ−Mn−1, (3.2) for any 1<ρ<ρ∗, where Cρ=2 ρ−1|G|1/2, and K =max(x2,···,xn)∈Gkf(x2,···,xn)kL2(Q). Therefore, the RPOD expansion converges with spectral accuracy in terms of the number of truncation modes in the main, secondary and tertiary expansions. Our analysis is based upon the estimates for rate of convergence of polynomial interpolation of vector-valued functions. We shall adapt a result by S. Bernstein in 1912, improved since then in many works (see for instance [20]). For some ρ>1, let the ellipse Eρin the complex plane be defined as Eρ=ζ∈C:| |ζ−1|+|ζ+1|≤ρ+ρ−1. Consider a function F:Eρ→Hwhere His a Hilbert space. For a given integer number M≥0 let FMbe the truncated Chebyshev polynomial series expansion of Fof degree M with coefficients in H. Then we have (cf. [24, Theorem 8.2]), Lemma 3.1. Assume that F is analytic and bounded in Eρ. There holds that max ξ∈[−1,1]kF(ξ)−FM(ξ)kH≤Cρρ−M, where Cρ=2 ρ−1kFkL∞(Eρ,H). Proof of Theorem 3.1.Let us consider the truncated primary expansion fM1(x1,x2,x3,···,xn)= M1 ∑ i1=0 σi1vi1(x2,x3,···,xn)ϕi1(x1), for some integer M1≥0. We adapt the proof of Lemma 0.10 in [5]. By hypothesis, for any given (x2,x3,···,xn−1)∈G2×G3×···×Gn−1, the vector-valued function ζ∈G47→ f(·,x2,x3,···,xn−1,ζ)∈L2(Q) M. Aza¨ıez et al. / Commun. Comput. Phys., 24 (2018), pp. 1556-1578 1563 is analytic in the ellipse Eρ. For brevity of notation we denote this function by f(x2,x3,···,xn−1,ζ). Consider the change of coordinates ζ=τ(ˆ ζ):=¯ γ−γ 2ˆ ζ+¯ γ+γ 2,ˆ ζ∈Eρ. It is affine and bijective from Eρinto Eρand transforms the reference interval [−1,1] into [γ,¯ γ]. Let ˆ S(x2,x3,···,xn−1) M1(ˆ ζ)be the truncated Chebyshev series expansion of the (transformed) function ˆ f(x2,x3,···,xn−1)(ˆ ζ) = f(x2,x3,···,xn−1,ζ). Let S(x2,x3,···,xn−1) M1(ζ) = ˆ S(x2,x3,···,xn−1) M1(ˆ ζ). From Lemma 3.1 we obtain max γ∈[γ,¯ γ]kf(x2,x3,···,xn)−S(x2,x3,···,xn−1) M1(x4)kL2(Q) ≤2 ρ−1kf(x2,x3,···,xn−1,·)kL∞(Eρ,L2(Q))ρ−M1≤2K ρ−1ρ−M1. (3.3) In view of the best approximation property stated in Lemma 2.2, it holds kf−fM1kL2(G×Q)≤kT−SM1kL2(G×Q) ≤|G|1/2 max (x2,x3,···,xn)∈Gkf(x2,x3,···,xn)−S(x2,x3,···,xn−1) M1(xn)kL2(Q). (3.4) By estimate (3.3) it follows that kf−fM1kL2(G×Q)≤CρKρ−M1, where Cρ=2|G|1/2 ρ−1. (3.5) Let us now consider the secondary expansion f(M1,M2)(x1,x2,x3,···,xn)= M1 ∑ i1=0 M2 ∑ i2=0 σi1σ(i1) i2ϕi1(x1)ϕ(i1) i2(x2)v(i1) i2(x3,···,xn). As the sequence (ϕi1)m≥0is orthonormal in L2(Q), then kfM1−f(M1,M2)k2 L2(G×Q)≤ M1 ∑ i1=0 σ2 i1kvi1−v(M2) i1k2 L2(G), where v(M2) i1(x2,x3,···,xn)= M2 ∑ i2=0 σ(i1) i2ϕ(i1) i2(x2)v(i1) i2(x3,···,xn) is the truncated POD expansion of vi1to M2+1 terms. 1564 M. Aza¨ıez et al. / Commun. Comput. Phys., 24 (2018), pp. 1556-1578 As vi1(x2,x3,···,xn)= 1 σi1ZQf(x1,x2,x3,···,xn)ϕi1(x1)dx1dt (3.6) then vi1is an analytic function from Ginto R. By an argument similar to the one that has lead to estimate (3.3), we prove that for any (x2,x3,···,xn−1)∈G2×G3×·×Gn−1there exists a polynomial in xn,S(i1,x2,x3,···,xn−1) M2(xn), of degree less or equal than M2such that max xn∈Gn|vi1(x2,x3,···,xn)−S(i1,x2,x3,···,xn−1) M2(xn)| ≤2 1−ρkvi1(x2,x3,···,xn−1,·)kL∞(E⊂)ρ−M2. (3.7) From (3.6), we deduce |σi1||vi1(x2,x3,···,xn)| ≤ kf(·,x2,x3,···,xn)kL2(Q)for all (x2,x3,···,xn)∈G, and then |σi1|max (x2,x3,···,xn)∈Gvi1(x2,x3,···,xn)−S(i1,x2,x3,···,xn−1) M2(xn)≤2 1−ρKρ−M2, Consequently, by the best approximation property, |σi1|kvi1−v(M2) i1kL2(G)≤|σi1|kvi1−S(i1) M2kL2(G) ≤|G|1/2 |σi1| kvi1−S(i1) M2kL∞(G)≤CρKρ−M2, (3.8) where S(i1) M2(x2,x3,···,xn)=S(i1,x2,x3,···,xn−1) M2(xn). Therefore kfM1−f(M1,M2)kL2(G×Q)≤CρpM1Kρ−M2. (3.9) Let us consider next the tertiary expansion, f(M1,M2,M3)= M1 ∑ i1=0 M2 ∑ i2=0 M3 ∑ i3=0 σi1σ(i1) i2σ(i1,i2) i3ϕi1(x1)ϕ(i1) i2(x2)ϕ(i1,i2) i3(x3)v(i1,i2) i3(x4,···,xn). As the functions ϕi1(x1)ϕ(i1) i2(x2)are orthogonal in L2(G1×G2), it holds kf(M1,M2)−f(M1,M2,M3)k2 L2(G×Q)≤ M1 ∑ i1=0 M2 ∑ i2=0 σ2 i1|σ(i1) i2|2kvi1 i2−v(i1,M3) i2k2 L2(G3×···×Gn), (3.10) where v(i1,M3) i2(x3,···,xn)= M3 ∑ i3=0 σ(i1,i2) i3ϕ(i1,i2) i3(x3)v(i1,i2) i3(x4,···,xn) M. Aza¨ıez et al. / Commun. Comput. Phys., 24 (2018), pp. 1556-1578 1571 5.2 Practical implementation In this subsection we analyze the determination of the number of modes required to meet a prescribed accuracy by the truncated RPOD expansion. This is based upon the following lemma. Lemma 5.1. Assume that the coefficients of the RPOD expansion (2.4) satisfy the estimates ∑ i1>I1σi1 2≤Σ1,∑ ik>I(i1,···,ik−1) kσ(i1,···,ik−1) ik 2≤Σk,k=2,···,n−1. (5.18) Then, the following error estimate in L2norm holds: f−fP 2 L2(Xn 1)≤Σ1+CPf 2 L2(Xn 1), (5.19) where CP=Σ2+···+Σn−2+Σn−1. We omit the proof of this lemma for brevity. Estimate (5.19) allows to derive a strategy for the practical computation of the RPOD expansion (2.4), if estimates (5.18) for the singular values are asymptotically sharp. Indeed assume that we want to determine the number of modes for truncation, M,Iand Jfor the solution of the advection-diffusion- reaction equation, such that ku−u(M1,M2,M3)kL2(G×Q)≤ε for some targeted error ε>0. By estimate (5.19), it holds ku−u(M1,M2,M3)kL2(G×Q)≤Σ1+(Σ2+Σ3)kuk2 L2(G×Q), and Σ1=∑ m≥1|σm|2,Σ2=sup m≥1 ∑ i≥I+1|σ(m) i|2,Σ3=sup m,i≥1 ∑ j≥J+1|σ(m,i) j|2. If estimates (5.10), (5.13) and (5.14) are sharp, then good estimators for Σ1,Σ2and Σ3are provided by Σ1≃S1=|σM+1|2,Σ2≃S2=σ(M) I+1 2,Σ3≃S3=σ(M1,M2) J+1 2. This suggests to recursively determine M,Iand Jin such a way that S1<ε2/3, S2<ε2/(3kuk2 L2(G×Q)),S3<ε2/(3kuk2 L2(G×Q)). This is possible if the truncated primary, secondary and tertiary expansions are constructed term by term, for instance by the power iteration algorithm (see [10]). 1572 M. Aza¨ıez et al. / Commun. Comput. Phys., 24 (2018), pp. 1556-1578 6 Numerical experiments This section is devoted to determine the effective convergence rate of the RPOD approximation of some analytic functions, so as of the solutions to the advection-diffusion- reaction equation (4.1) when parameterized by the diffusivity, advection velocity and reaction coefficients. 6.1 Numerical validation of RPOD expansion with toy models To assess the efficiency of the RPOD to approximate multidimensional data let us consider tensors sampled from the following test functions: f1(x1,x2,x3,x4,x5)=sin(x1x4x2 5)+p9+x3+4x2, (6.1) f2(x1,x2,x3,x4,x5)= 1 1+x1+x2+x3+x4+x5 , (6.2) f3(x1,x2,x3,x4,x5)= x2 1{sin[5x2π+3log(x3 1+x2 2+x2 4+x3+π2)]−1}2 +(x1+x3−1)(2x2−x3)(4x5−x4)cos[30(x1 +x3+x4+x5)]log(6+x2 1x2 2+x3 3) −4x2 1x2x3 5(1−x3)3/2. (6.3) The spatial domain is the unit hypercube, and the toy model tensors are constructed by discretization these functions in a structured equispaced 10×10×10×10×10-mesh. Thus, we have three 5-order tensors whose size is 105. Note that f1and f2are fairly simple, but f3is more complex since it is constructed through algebraic operations plus composition of various transcendental functions. Also note that because of that term proportional to (1−x3)3/2, the function is not analytic. Fig. 1 shows the L2(X1)errors of the RPOD expansion (5.8) in semi logarithmic scale for the three functions f1,f2and f3. All three curves are consistent with the exponential convergence, modified by the factor Ln−2 2(n−1), that yields a concave curve for Llarge enough. The concavity of the error curves is larger for smaller convergence rate, a characteristic of the function that yields the error estimate, e(L)= D′ nLn−2 2(n−1)ρ−L1 n−1. (6.4) In order to test (6.4), which should hold for enough large values of L, we fit the two parameters D′ nand ρ. Specifically, in the fitting process, we consider for functions f1,f2 and f3the modes for which L1/4 (in this case n=5) is larger than 2, 3 and 6.5, respectively. The graphical result of this fit, shown in figure 1, confirms our expectations. In fact, in all the considered cases, the experimental values approach the ones found out by fitting the function (6.4). The values obtained for ρrespectively are 4.54×104, 1.06×102 and 1.26×104. These values are quite large, so the exponential factor in (6.4) is largely M. Aza¨ıez et al. / Commun. Comput. Phys., 24 (2018), pp. 1556-1578 1573 Figure 1: Convergence history of RPOD expansion error in logarithmic scale. In blue are plotted the experimental values and in red the theoretical ones. The x-axis corresponds to L1/4, where Lis the number of modes Lin the re-arranged expansion. dominant and the fitting curves are quite close to a straight line. These results are in good agreement with the error estimate (5.8), suggesting that it could be asymptotically sharp. 6.2 Advection-diffusion-reaction equations In this sub-section we study the effective convergence rate of the RPOD approximation of the solution to de advection-diffusion-reaction equation. We assess the exponential convergence rate and investigate the variation of this rate with respect to the set G= BR×[γ,¯ γ]×[α,¯ α]. We consider the time-dependent advection-diffusion-reaction equation in the domain Q=(0,1)×(0,1)and the source term and initial condition, given by f(t,x)=q|x−t−0.3|,T0(x)=|x−0.4|. These data have mild singularities, so the temperature solutions of (4.1), have a reduced regularity with respect to xand t, in particular for t=0. The data we are considering to assess the efficiency of the RPOD approach are provided as a solution of the advectiondiffusion-reaction problem which is discretized by an implicit Euler second order scheme for the time derive. The advection term is approximated using an explicit second order Adams-Bashforth [15] scheme. Gauss-Lobatto-Legendre spectral method is used for the space approximation see [9]. Various integrals (with respect to either γ,α,~aor x) are computed using Gauss-Lobatto quadrature formulas with high resolution in the corresponding intervals. The numerical experiments are achieved using time step δt=10−3 and the polynomial degree N=24. 1574 M. Aza¨ıez et al. / Commun. Comput. Phys., 24 (2018), pp. 1556-1578 Test 1: Exponential convergence rate. Fig. 2 shows the convergence history of the RPOD expansion for the advection-diffusion-reaction equation in terms of the total number of modes in the expansion. We have considered the sets of diffusivities γ∈[1,51], advection velocity a∈[0,30]and reaction rates α∈[0,100]. The error is measured in L2(G×Q)norm. The numbers of principal and secondary modes have been determined to be σ(Im+1) m≤ε=10−5. The modes have been re-arranged in decreasing order of the effective singular values. We observe that the singular values indeed are good error estimators for this re-arranged expansion, as was argued in Subsection 5.2. 050 100 150 200 L 10-5 10-4 10-3 10-1 100 Singular Values L2-error 12 3 4 56 M 10-5 10-4 10-3 10-1 100 Singular Values L2-error Figure 2: Convergence history for RPOD expansion of the solution of the advection-diffusion-reaction equation. Left : with respect to L, right with respect to M=L1/3 Test 2: Dependence of the convergence rate with respect to the diffusion parameters range. The dependence with respect to the ratio of diffusivities r=¯ γ/γof the exponential convergence rate, stated by Theorem 3.1, is illustrated in Fig. 3. We depict the convergence history for our, computed for r=2, 10, 25, 100, 400 in all cases with a fixed interval of velocity advection [a,¯ a] = [0,5]and reaction rates [α,¯ α] = [0,10], with respect 11,5 22,5 33,5 M 10-5 10-4 10-3 10-1 r = 2 r = 10 r = 25 r = 100 r = 400 Figure 3: Variation of the RPOD errors (in logarithmic scale) with respect to the ratio r=¯ γ/γ, for fixed [am,¯ a]=[0,5]and α=0,¯ α=10. The variable Mstands for L1/3 where Lis the number of modes. M. Aza¨ıez et al. / Commun. Comput. Phys., 24 (2018), pp. 1556-1578 1575 to M=L1/3. We can point out that the convergence rate degrades as rincreases. This is consistent with the dependence upon rof estimate (4.7). Test 3: Dependence of the convergence rate with respect to the reaction parameters range. We show in Fig. 4 the convergence rates history corresponding to α=0, ¯ α= 1,25,50,75,100 for fixed γ=1, ¯ γ=51 and a=0, ¯ a=5. We observe a decrease of the rate as ¯ αincreases, that however appears to be uniformly bounded, in agreement with estimate (4.7). 11,5 22,5 33,5 44,5 M 10-5 10-4 10-3 10-2 10-1 (0,10) (0,100) (0,400) (0,700) (0,1000) Figure 4: Variation of the RPOD errors with respect to the reaction rate α. The curves correspond to α=0, ¯ α=1,25,50,100,200 with γ=1,¯ γ=51 and a=0,¯ a=20 in all cases. The variable Mstands for L1/3 where Lis the number of modes. Test 4: Dependence of the convergence rate with respect to the advection parameters range. We continue our study by testing the dependence of the convergence rate with respect to the range of the advection velocity [a,¯ a]. We show in Fig. 5 the convergence rates history corresponding to a=0, ¯ a=1,25,50,100,200 for fixed γ=1, ¯ γ=51 and α=0, 12 3 4 5 M 10-5 10-4 10-3 10-1 100 R = 1 R = 10 R = 20 R = 30 R = 40 Figure 5: Variation of the RPOD errors with respect to the advection velocity variation R=¯ a. The curves correspond to a=0,R=1,10,20,30,40 and 50, with γ=50,¯ γ=101 and α=0,¯ α=10 in all cases. The variable Mstands for L1/3 where Lis the number of modes. 1576 M. Aza¨ıez et al. / Commun. Comput. Phys., 24 (2018), pp. 1556-1578 ¯ α=50. We can point out that the convergence rate degrades as ¯ aincreases, but the slope of the convergence curves remains fixed. This also is strongly consistent with estimate (4.7). The last numerical experiment studies whether the dependence of the exponential convergence on the diffusivities range [γ,¯ γ]indeed takes place in terms of the ratio r=¯ γ/γ. This is confirmed by the result plotted in Fig. 6, where we consider the couples (γ,¯ γ) = (1,2)and (4,8), corresponding to r=2, and (γ,¯ γ) = (1,25)and (4,100), corresponding to r=25, with fixed a=0, ¯ a=5 and α=0, ¯ α=10. We observe that the slope of the convergence curve for r=2 is larger than for r=25 (for large values of M), in agreement with estimate (4.7). 0,9 1,2 1,5 1,8 2,1 2,4 2,7 M 10-6 10-4 10-2 100 (1,2) (4-8) (1-25) (4-100) Figure 6: Analysis of dependence of the RPOD errors with respect to the ratio r=¯ γ/γ. The curves correspond to the indicated pairs (γ,¯ γ). The variable Mstands for L1/3 where Lis the number of modes. 7 Conclusion We have studied in this paper the RPOD expansion of the solutions of the advectionreaction-diffusion problem. We have demonstrated the exponential rate of convergence of the RPOD expansion for its solution. The proof is based upon the analyticity of its solution with respect to the parameters appearing in the equation. We have also applied the RPOD to approximate analytic functions, that also show an exponential convergence rate. We conclude the ability of the RPOD expansion to construct low-rank approximations of smooth multi-parametric functions, depending on a moderate number of parameters. Acknowledgments The research of Tom´as Chac´on has been partially supported by the Spanish Government – EU FEDER Fund project MTM2015-64577-C2-1-R, and by the IDEX program of the Universit´e de Bordeaux. The research of Mejdi Azaiez has been partially supported by the M. Aza¨ıez et al. / Commun. Comput. Phys., 24 (2018), pp. 1556-1578 1577 Spanish Government EUFEDER Fund project MTM2015-64577-C2-1-R by NSFC/ANR joint program 51661135011. The research of Emma Perracchione is accomplished within Rete ITaliana di Approssimazione (RITA), GNCS-INdAM and supported by the research project Data Assimilation Methods for Machine Learning with application to environmental data, founded by UE ERA-PLANET GA n. 689443. References [1] R.A. ADAMS, J.J.F FOURNIER,Sobolev Spaces, Academic Press, 2003. [2] D. V. ANOSOV, S. KH. ARANSON, V. I. ARNOLD, I. U. BRONSHTEIN, V. Z. GRINES, YU. S. IL’YASHENKO,Ordinari Differential Equations and Smooth Dynamical Systems. Springer (1991). [3] AZA¨ IEZ, M. AND BEN BELGACEM, F., Karhunen-Lo`eve’s truncation error for bivariate functions, Computer Methods in Applied Mechanics and Engineering, Vol 290, pp 57-72, 2015 [4] M. AZA¨ IEZ, F. BEN BELGACEM AND T. CHAC ´ ON REBOLLO., Error Bounds for POD expansions of parameterized transient temperatures CMAME. 305 (2016) 501-511 [5] M. AZA¨ IEZ, F. BEN BELGACEM, T. CHAC ´ ON REBOLLO,Recursive POD expansion for reactiondiffusion equation, Adv.Model. and Simul. in Eng. Sci. (2016) 3:3. DOI 10.1186/s40323-016-0060- 1 [6] M. BALAJEWICZ,A new approach to model order reduction of the Navier-Stokes Equations, PhD thesis, Duke University, Durham, 2012. [7] D. BELTON,Improving and extending the information on principal component analysis for local neighborhoods in 3D point clouds, The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 37 (2008), B5: 477 ff. [8] G. BERKOOZ, P. HOLMES, J.L. LUMLEY,The proper orthogonal decomposition in the analysis of turbulent flows, Annu. Rev. Fluid Mech. 25 (1993), pp. 539–575. [9] Bernardi, C. and Maday, Y., Approximations spectrales de probl`emes aux limites elliptiques, Math´ematiques et applications, Springer, Paris, Berlin, Heidelberg, 1992. [10] F. CHINESTA, A. AMMAR, A. LEYGUE, R. KEUNINGS,An overview of the proper generalized decomposition with applications in computational rheology, J. Non-Newtonian Fluid Mech. 166 (2011), pp. 578–592. [11] A. COHEN AND R. DEVORE,Kolmogorov widths under holomorphic mappings, IMAJ. Numer. Anal., 36 (2016), pp. 112. [12] A. COHEN, R. DEVORE AND CH. SCHWAB,Analytic regularity and polynomial approximation of parametric and stochastic elliptic PDEs, Analysis and Applications. Vol. 09, Issue 01, January 2011 [13] L. DELATHAWER, B. DEMOOR, J. VANDEWALLE,On the best rank-1and rank-(R1,R2,···,RN) approximation of higher-order tensors, SIAM J. Matrix. Anal. Appl. 21 (2000), pp. 1324–1342. [14] L. DELATHAWER, B. DEMOOR, J. VANDEWALLE,A multilinear singular value decomposition, SIAM J. Matrix. Anal. Appl. 21 (2000), pp. 1253–1278. [15] C. A. J. FLETCHER,Computational Techniques for Fluid Dynamics 1 Fundamental and General Techniques, Springer-Verlag. 1988. [16] G.H. GOLUB, C.F. VAN LOAN,Matrix Computations, Johns Hopkins University Press, Baltimore, MD, 1996. [17] C. HEYBERGER, P.A. BOUCARD, D. N´ ERON,Multiparametric analysis within the proper generalized decomposition framework, Comput. Mech. 49 (2012), pp. 277-289. [18] I.T. JOLLIFFE,Principal Component Analysis, New York, Springer Verlag, 1986. 1578 M. Aza¨ıez et al. / Commun. Comput. Phys., 24 (2018), pp. 1556-1578 [19] M. KIRBY, L. SIROVICH,Application of the Karhunen-Lo`eve procedure for the characterization of human faces, IEEE Transactions on Pattern Analysis and Machine Intelligence, 12 (1990), pp. 103–108. [20] G. LITTLE, J. B. READE,Eigenvalues of analytic kernels. SIAM J Math Anal. 15 (1984), pp. 133–6. [21] M.M. LO` EVE,Probability Theory, Van Nostrand, Princeton, NJ, 1988. [22] J. MERCER,Functions of positive and negative type and their connection with the theory of integral equations, Phil. Trans. Royal Society, 209 (1909), pp. 415–446. [23] M. M ¨ ULLER,On the POD Method: An abstract investigation with applications to reduced-order modeling and suboptimal control, PhD thesis, Georg-August Universit¨at, G¨ottingen, 2008. [24] Trefethen, L. N., Approximation theory and approximation practice. Software, Environments, and Tools., Society for Industrial and Applied Mathematics (SIAM), Philadel-phia, PA, 2013. [25] S. ULLMANN, J. LANG,POD-Galerkin modeling and sparse-grid collocation for a natural convection problem with stochastic boundary conditions, Lecture Notes in Computational Science and Engineering 97, pp. 295–315. [26] N. VANNIEUWENHOVEN, R. VADEBRIL, K. MEERBERGEN A new truncation strategy for the higher-order singular value decomposition, SIAM J. SCi. Comput, Vol. 34, No. 2, pp. A1027A1052, 2012. [27] S. VOLKWEIN,Proper Orthogonal Decomposition and Singular Value Decomposition. Rapport technique 153, Institut fr Mathematik, universit de Graz, 1999. [28] H. WOLFGANG,Tensor spaces and numerical tensor calculus. Springer series in computational mathematics, V 42, 2012