Reduced Basis Approximations of the Solutions to Spectral Fractional Diffusion Problems

TL;DR

This paper introduces a reduced basis method to efficiently approximate solutions to spectral fractional diffusion problems, significantly reducing computational costs while maintaining accuracy across different fractional powers.

Contribution

The work develops a fractional power-independent reduced basis strategy for spectral fractional diffusion, enabling fast online evaluations with proven exponential convergence.

Findings

Reduced basis method accelerates reaction-diffusion problem solutions.

The approach is effective uniformly over fractional powers in [s_min, s_max].

Numerical experiments confirm exponential convergence.

Abstract

We consider the numerical approximation of the spectral fractional diffusion problem based on the so called Balakrishnan representation. The latter consists of an improper integral approximated via quadratures. At each quadrature point, a reaction-diffusion problem must be approximated and is the method bottle neck. In this work, we propose to reduce the computational cost using a reduced basis strategy allowing for a fast evaluation of the reaction-diffusion problems. The reduced basis does not depend on the fractional power $s$ for $0<s_{\min}\leq s \leq s_{\max}<1$. It is built offline once for all and used online irrespectively of the fractional power. We analyze the reduced basis strategy and show its exponential convergence. The analytical results are illustrated with insightful numerical experiments.

Equations

$$u(s)\text{mathchar}58=(-\Delta {)}^{-s}f\text{mathchar}58=\sum _{k=1}^{\infty }{\lambda }_k^{-s}{f}_k{\psi }_k,$$

$$u(s)=\frac{\sin (s\pi )}{\pi }{\int }_{-\infty }^{\infty }{e}^{(1-s)y}w(y)dy,$$

$$(e^{y}I-\Delta )w(y)=f.$$

$$u(s)\approx u_{h,k}(s)=\frac{k\sin (s\pi )}{\pi }\sum _{l=-M_s}^{N_s}{e}^{(1-s)y_l}{w}_h(y_l)$$

$$M_s\text{mathchar}58=\lceil \frac{{\pi }^2}{(1-s)k^2}\rceil ,N_s\text{mathchar}58=\lceil \frac{{\pi }^2}{s{k}^2}\rceil$$

$$u_{h,k}^n(s)\text{mathchar}58=\frac{k\sin (s\pi )}{\pi }\sum _{l=-M_s}^{N_s}{e}^{(1-s)y_l}{w}_h^n(y_l).$$

$$\mathbb{H}^{r}(\Omega)\mathrel{\mathop{\mathchar 58\relax}}=\left\{\begin{array}[]{ll}\left(L^{2}(\Omega),H_{0}^{1}(\Omega)\right)_{r}&\mbox{for }r\in[0,1]\\ H_{0}^{1}(\Omega)\cap H^{r}(\Omega)&\mbox{for }r\in(1,2],\end{array}\right.$$

$$\Vert v{\Vert }_{H_0^1(\Omega )}\text{mathchar}58=\Vert v{\Vert }_{{\mathbb{H}}^1(\Omega )}=\Vert \nabla v{\Vert }_{L^2(\Omega )}\forall v\in H_0^1(\Omega ),$$

$$\Vert v{\Vert }_{L^2(\Omega )}\le C_P\Vert \nabla v{\Vert }_{L^2(\Omega )}\forall v\in H_0^1(\Omega ).$$

$$u(s)=\frac{1}{2\pi i}{\int }_{\mathcal{C}}{z}^{-s}(zI+\Delta {)}^{-1}fdz,$$

$$u(s)=\frac{\sin (s\pi )}{\pi }{\int }_0^{\infty }{\mu }^{-s}(\mu I-\Delta {)}^{-1}fd\mu .$$

$$u(s)=\frac{\sin (s\pi )}{\pi }{\int }_{-\infty }^{\infty }{e}^{(1-s)y}(e^{y}I-\Delta {)}^{-1}fdy.$$

$$u_h(s)\text{mathchar}58=\frac{\sin (s\pi )}{\pi }{\int }_{-\infty }^{\infty }{e}^{(1-s)y}{w}_h(y)dy,$$

$$a(w_h(y),v_h;y)=F(v_h),\forall v_h\in {\mathbb{V}}_h.$$

$$a(w,v;y)\text{mathchar}58=a_0(w,v)+e^y{a}_1(w,v)\text{mathchar}58={\int }_{\Omega }\nabla w\cdot \nabla v+e^y{\int }_{\Omega }wv$$

$$F(v)\text{mathchar}58={\int }_{\Omega }fv,$$

$$\Vert v{\Vert }_{H_0^1(\Omega )}^2\le a(v,v;y)\text{and}a(v,w;y)\le (1+C_P^2{e}^y)\Vert v{\Vert }_{H_0^1(\Omega )}\Vert w{\Vert }_{H_0^1(\Omega )},$$

$$\Vert \nabla w_h(y)\Vert \le C_P\Vert f\Vert ,\Vert w_h(y)\Vert \le e^{-y}\Vert f\Vert ,$$

$$\Vert \nabla (w_h(y)-w_h(\bar{y}))\Vert \le C_P^3|e^y-e^{\bar{y}}|\Vert f\Vert ,$$

$$\Vert w_h(y)-w_h(\bar{y})\Vert \le e^{-y}|e^{y-\bar{y}}-1|\Vert f\Vert ,$$

$$\Vert \nabla (w_h(y)-w_h(\bar{y}))\Vert \le \frac{1}{2}{e}^{-\frac{y}{2}}|e^{y-\bar{y}}-1|\Vert f\Vert .$$

$$\Vert \nabla w_h(y){\Vert }^2+e^y\Vert w_h(y){\Vert }^2={\int }_{\Omega }f{w}_h(y)\le \Vert f\Vert \Vert w_h(y)\Vert ,$$

$${\int }_{\Omega }\nabla (w_h(y)-w_h(\bar{y}))\cdot \nabla v_h+e^y{\int }_{\Omega }(w_h(y)-w_h(\bar{y}))v_h=(e^{\bar{y}}-e^y){\int }_{\Omega }{w}_h(\bar{y})v_h\forall v_h\in {\mathbb{V}}_h.$$

$$\Vert \nabla (w_h(y)-w_h(\bar{y})){\Vert }^2+e^y\Vert w_h(y)-w_h(\bar{y}){\Vert }^2\le |e^{\bar{y}}-e^y|\Vert w_h(\bar{y})\Vert \Vert w_h(y)-w_h(\bar{y})\Vert .$$

$$\Vert \nabla (w_h(y)-w_h(\bar{y}))\Vert \le \frac{1}{2}{e}^{-\frac{y}{2}}|e^{\bar{y}}-e^y|\Vert w_h(\bar{y})\Vert .$$

$$\Vert u-u_h{\Vert }_{{\mathbb{H}}^r(\Omega )}\lesssim \ln (h^{-1})h^{2{\alpha }^{\ast }}\Vert f{\Vert }_{{\mathbb{H}}^{r+2{\alpha }^{\ast }-2s}(\Omega )}.$$

$$\Vert u-u_h{\Vert }_{{\mathbb{H}}^r(\Omega )}\lesssim h^{2{\alpha }^{\ast }}\Vert f{\Vert }_{{\mathbb{H}}^{r+2{\alpha }^{\ast }-2s+2\epsilon }(\Omega )}.$$

$$\Vert u-u_h{\Vert }_{{\mathbb{H}}^r(\Omega )}\lesssim h^{2{\alpha }^{\ast }}\Vert f{\Vert }_{L^2(\Omega )}.$$

$$\Vert u_h(s)-u_{h,k}(s){\Vert }_{{\mathbb{H}}^r(\Omega )}\lesssim \left(e^{-\frac{{\pi }^2}{k}}+e^{-(1-s)M_{s}k}+e^{-s{N}_{s}k}\right)\Vert f{\Vert }_{L^2(\Omega )}.$$

$$\Vert u_h(s)-u_{h,k}(s){\Vert }_{{\mathbb{H}}^r(\Omega )}\lesssim \left(e^{-\frac{{\pi }^2}{k}}+e^{-(1-s)M_{s}k}+e^{-s{N}_{s}k}\right)\Vert f{\Vert }_{{\mathbb{H}}^{r-2s+\epsilon }(\Omega )}.$$

Full text

This research was supported by the NSF grant DMS-1817691 (AB-DG-AZ); DG was supported by the Swiss National Science Foundation grant P2ELP2-175056 and IAMCS at TAMU

1. Introduction

Nonlocal models have recently received a great attention due to their apparent ability to capture novel effects such as in mechanics [27] and in particular in peridynamics [26, 19], turbulence [10], biophysics [7] and image denoising [17], to mention a few.

In most of the applications, the type of nonlocal interactions are different and their scaling laws are unknown. Initiated by the work in [2], an algorithm is proposed and analyzed in [1] to identify the fractional power $s\in[s_{\min},s_{\max}]$ governing the state equation in an optimization framework. As expected, the algorithm exploits the smoothness of the map $s\mapsto(-\Delta)^{-s}f$ and requires many (costly) evaluations of $(-\Delta)^{-s}f$ for different $s\in(0,1)$. Several numerical methods are available to approximate $(-\Delta)^{-s}f$ and we refer to the surveys [3, 20] for the description of different fractional Laplacians along with their numerical approximations. Here, for $f\in L^{2}(\Omega)$ and $\Omega$ a Lipschitz domain of $\mathbb{R}^{d}$, $d=1,2,3$, we set

[TABLE]

where $\{\lambda_{k},\psi_{k}\}_{k\in\mathbb{N}}\subset\mathbb{R}^{+}\times H^{1}_{0}(\Omega)$ are the eigenpairs of $(-\Delta)$ and $f_{k}\mathrel{\mathop{\mathchar 58\relax}}=\int_{\Omega}f\psi_{k}$. The eigenfunctions $\{\psi_{k}\}_{k\in\mathbb{N}}$ are chosen orthogonal in $H_{0}^{1}(\Omega)$ and orthonormal in $L^{2}(\Omega)$. In (1) the fractional operator is referred to as the spectral fractional Laplacian and is the one considered in this work. It is worth mentioning that the methodology proposed here is not limited to the Laplacian operator and can be easily extended to regularly accretive operators as in [6].

In this work, we follow the approach proposed in [5] to approximate (1), see also [6], which is based on the Dunford-Taylor-Balakrishnan representation

[TABLE]

where $w(y)\in H^{1}_{0}(\Omega)$ solves

[TABLE]

Originally introduced in [5] and later improved in [4], a sinc quadrature coupled with a standard finite element method is used for the approximation of the integration in the variable $y$. For $k>0$, it reads

[TABLE]

with $y_{l}\mathrel{\mathop{\mathchar 58\relax}}=lk$,

[TABLE]

and where $w_{h}(y_{l})\in\mathbb{V}_{h}$ are standard finite element approximations of $w(y_{l})$.

The numerical approximation of $(-\Delta)^{-s}f$ requires $M_{s}+N_{s}+1$ finite element solves to determine $w_{h}(y_{l})$, $y_{l}\in[-M_{s}k,N_{s}k]$. This can become prohibitive when the computation of $(-\Delta)^{-s}f$ is needed for many values of $s$ such as within an optimization loop as mentioned above. The reduced basis method seems to be a natural approach to reduce the computational cost when approximating the parametrized reaction-diffusion problems (2). In fact, reduced basis method for this type of one dimensional parametric elliptic partial differential equation has already been partially analyzed in [22, 21] and recently in [13] from which part of our analysis is inspired.

A (weak) greedy strategy is advocated (offline stage) to iteratively select snapshots $w_{h}(y^{l})$, $l=1,...,n<\!\!<\dim(\mathbb{V}_{h})$, defining the $s$-independent reduced basis space $\mathbb{V}_{h}^{n}\mathrel{\mathop{\mathchar 58\relax}}=\operatorname{span}\{w_{h}(y^{1}),\ldots,w_{h}(y^{n})\}$. Galerkin approximations $w_{h}^{n}(y_{l})\in\mathbb{V}_{h}^{n}$ of $w_{h}(y_{l})$ can then be easily computed (online stage) to produce a reduced basis approximation of $u_{h,k}(s)$

[TABLE]

We point out that one of the difficulty faced in this study is that the approximation of the parametric elliptic partial differential equation (2) is required for $y$ in the parametric domain $[-M_{s}k,N_{s}k]$ whose length increases as the sinc quadrature parameter $k$ decreases (improving the precision of the algorithm).

The proposed algorithm provides an approximation of the entire map $s\mapsto u(s)$, $s\in[s_{\min},s_{\max}]$ using the same reduced basis space $\mathbb{V}_{h}^{n}$. Our main result is Theorem 4 which guarantees an exponential convergence of the reduced basis approximation $u_{h,k}^{n}(s)$ toward $u_{h,k}(s)$ in a wide range of Sobolev norms, uniformly in the fractional power $s\in[s_{\min},s_{\max}]$.

We end this introduction by noting that the idea of using the reduced basis method for fractional problems has been recently proposed for instance in [13] and [15]. In [13], the reduced basis is used for the approximation of interpolation norms, as well as evaluations of both types $s\mapsto(-\Delta)^{s}u$ with $u$ fixed and variable $s\in(0,1)$ and $u\mapsto(-\Delta)^{s}u$. A reduced basis space based on best rational approximations for a reaction diffusion problem similar to the one satisfied by $w_{h}$ is proposed and exponential convergence of the approximation with respect to the dimension of the reduced basis space is obtained. Worth mentioning, the numerical method is based on the extension method [23] but seemingly apply to other approximation techniques. Actually, this method boils down to the approximation of several reaction diffusion problems as in [5]. We take advantage of the technology developed in [13] to derive an exponential decay in the approximation of (3) by (5). In [15], a similar approximation $u_{h,k}(s)$ is proposed for a different quadrature. Exponential decay of the reduced basis error is observed numerically but without analysis. In some sense, this work provides a mathematical justifications of the experimental observations in [15]. Finally, we mention that the reduced basis method has also been used in [9] to approximate the parametric PDEs $(-\Delta)^{s}u=f$, where $s$ is the parameter, in the case of the integral fractional Laplacian.

The rest of the paper is organized as follows. In Section 2, we describe the numerical approximation of $u(s)$ by $u_{h,k}(s)$. Section 3 describes the construction of the reduced basis space and its corresponding error analysis - the main result of this work. Section 4 provides numerical experiments to illustrate the performance of the proposed methodology.

2. Spectral Fractional Laplacian and its Numerical Approximations

We start with some notations. Let $\mathbb{H}^{r}(\Omega)$ be the interpolation space defined by

[TABLE]

where $(\cdot,\cdot)_{r}$ denotes interpolation using the real method.

Notice that for the particular case $r=1$, we have

[TABLE]

which is equivalent to the $H^{1}(\Omega)$ norm thanks to the Poincaré inequality

[TABLE]

To simplify the notation, we write when $r=0$, $\|\cdot\|\mathrel{\mathop{\mathchar 58\relax}}=\|\cdot\|_{L^{2}(\Omega)}=\|.\|_{\mathbb{H}^{0}}$. Moreover, $a\lesssim b$ means that $a\leq Cb$ for a constant $C$ that does not depend on $a$, $b$ and the discretization parameters and whose value might change at each occurrence. Also, $a\approx b$ indicates $a\lesssim b$ and $b\lesssim a$.

2.1. Dunford-Taylor Representation

The function $u(s)\in L^{2}(\Omega)$ in (1) has the following representation [29]

[TABLE]

where $\mathcal{C}$ is a Jordan curve oriented to have the spectrum of $-\Delta$ to its right. Deforming the contour $\mathcal{C}$ to the negative real axis, we obtain the Balakrishnan formula, valid for $s\in(0,1)$,

[TABLE]

The numerical integration of the above improper integral relies on a sinc quadrature method after the change of variable $y=\ln(\mu)$, leading to

[TABLE]

2.2. Finite Element Approximation

We assume that $\Omega$ is a polyhedral domain and we consider a sequence $\{\mathcal{T}_{h}\}_{h>0}$ of conforming and shape-regular partitions of $\Omega$ into $d$-simplices with maximal mesh size $h<1$. Let $\mathbb{V}_{h}$ be the space of continuous and piecewise linear finite element functions associated with $\mathcal{T}_{h}$. The finite element approximation of (10) is then defined by

[TABLE]

where $w_{h}(y)\in\mathbb{V}_{h}$ is the solution to

[TABLE]

Here we used the notation

[TABLE]

for $w,v\in H_{0}^{1}(\Omega)$ and

[TABLE]

for $v\in H_{0}^{1}(\Omega)$. The Poincaré inequality (7) implies that for $v,w\in H^{1}_{0}(\Omega)$,

[TABLE]

which guarantees that (12) has a unique solution for any parameter $y\in\mathbb{R}$ by the Lax-Milgram lemma.

We now collect some estimates for $w_{h}(y)$, which will be used in the analysis later.

Lemma 2.1**.**

Let $C_{P}$ be the Poincaré constant in (7). For any $y,\bar{y}\in\mathbb{R}$, we have

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

Proof.

Choosing $v_{h}=w_{h}(y)$ in (12) yields

[TABLE]

from which the two relations in (16) can be easily deduced. From (12) we get

[TABLE]

We now choose $v_{h}=w_{h}(y)-w_{h}(\bar{y})$ to get

[TABLE]

This, the Poincaré inequality (7) and (16) with $y=\bar{y}$ yield (17).

The estimate (18) follows from (20) together with (16) with $y=\bar{y}$. For (19), we invoke Young’s inequality to estimate the right hand side of (20) and get

[TABLE]

It remains to invoke (16) with $y=\bar{y}$ to derive the desired result and ends the proof. ∎

We mention that both results in (16) are standard, while the estimates (17), (18) and (19) are less common yet useful in the error analysis below. Moreover, note that (16)-left and (17) are favorable for negative $y$, while (16)-right, (18) and (19) for positive $y$. This plays a role in our analysis below and is observed in the numerical experiments, see Section 4.

We end this section by stating the error in the finite element method derived and analyzed in [6]. Before doing this, we define $\alpha\in(0,1]$ to be the elliptic pick-up regularity index, i.e. $\alpha$ is the largest number in $(0,1]$ such that $(-\Delta)$ is an isomorphism from $\mathbb{H}^{r}(\Omega)$ to $\mathbb{H}^{r+1}(\Omega)$ for all $r\in[0,\alpha]$. Notice that $\alpha>0$ for Lipschitz domains and $\alpha=1$ when $\Omega$ is convex.

Theorem 1**.**

Let $f\in L^{2}(\Omega)$, $\alpha>0$ denote the elliptic regularity pick-up and $\alpha^{*}\mathrel{\mathop{\mathchar 58\relax}}=\frac{\alpha+\min(\alpha,1-r)}{2}$. Then for any $r\in[0,1]$, we have

If $r+2\alpha^{*}-2s\geq 0$ and $f\in\mathbb{H}^{r+2\alpha^{*}-2s}(\Omega)$ then

[TABLE] 2. 2.

If $r+2\alpha^{*}-2s\geq 0$ and $f\in\mathbb{H}^{r+2\alpha^{*}-2s+2\varepsilon}(\Omega)$ with $r+2\alpha^{*}-2s+2\varepsilon\leq 1+\alpha$ then

[TABLE] 3. 3.

If $r+2\alpha^{*}-2s<0$ then

[TABLE]

2.3. Sinc quadrature approximation

We discuss the sinc quadrature approximation leading to the fully discrete approximation $u_{h,k}$ given by (3). Recall that $k>0$ is the sinc quadrature parameter and that $M_{s}$, $N_{s}$ are given by (4). This choice is dictated from the analysis of the sinc quadrature error, which is the subject of the following theorem; we refer to [4] for its proof.

Theorem 2**.**

Let $f\in L^{2}(\Omega)$ and $r\in[0,1]$.

If $s>r/2$ then

[TABLE] 2. 2.

If $s\leq r/2$ and $f\in\mathbb{H}^{r-2s+\varepsilon}(\Omega)$ with $r-2s+\varepsilon\in[0,1+\alpha]$ then

[TABLE]

3. Reduced Basis Approximation

The computation of $u_{h,k}$ in (3) involves the finite element solution $w_{h}(y)$ to the reaction diffusion problem (12) for $M_{s}+N_{s}+1$ different values of the parameter $y\in\mathcal{D}_{s}\mathrel{\mathop{\mathchar 58\relax}}=[-M_{s}k,N_{s}k]$. We propose to use the reduced basis method to approximate the entire map $y\mapsto w_{h}(y)$. Notice that the bilinear form $a(\cdot,\cdot;y)$ defining $w_{h}(y)$ in (12) is not affine in $y$. However, it becomes affine for $\mu=e^{y}$.

3.1. Construction for a fixed $s$

The reduced basis space

[TABLE]

is constructed using a greedy strategy [8]. Starting with $y^{1}=0$, $y^{m+1}\in\mathcal{D}_{s}$ is selected iteratively to maximize the error

[TABLE]

on $\mathcal{D}_{s}$, i.e.,

[TABLE]

Here $P_{\mathbb{V}_{h}^{m}}w_{h}(y)\in\mathbb{V}_{h}^{m}$ is the unique solution (from Lax-Milgram theory) to

[TABLE]

Notice that in view of the definition (12) of $w_{h}(y)$, the above relation is equivalent to

[TABLE]

The enrichment of the reduced basis space ends when

[TABLE]

for a prescribed accuracy $\varepsilon>0$ or when a maximum number of basis functions $N_{\max}$ is reached. Note that the relations (17) and (19) guarantee that (23) can always be achieved by a uniform selection of points $y$ in $\mathcal{D}_{s}$. The aim of the Greedy algorithm is to provide an alternate selection performing as well but of less cardinality.

The error $W_{m}(y)$ defined in (21) is not a computable quantity and is usually replaced by an equivalent computable quantity leading to the so-called weak greedy algorithm [14]. In this work, we use the residual based a posteriori error estimate [24]

[TABLE]

where

[TABLE]

We have the following equivalence relation between the error $W_{m}(y)$ and its surrogate

[TABLE]

for all $y\in\mathbb{R}$. This follows from

[TABLE]

where $w_{h}(y)$ satisfies (12), and the coercivity and continuity of the bilinear form $a$ (15).

Hence, selecting the samples using $\|r_{m}(\cdot;y)\|_{\mathbb{V}_{h}^{\prime}}$ as surrogate for the error $W_{m}(y)$ yields

[TABLE]

where

[TABLE]

The parameter $\gamma_{s}$ corresponds to the constant in the weak greedy algorithm, see [14] for more details, and will appear in the analysis below.

The dual norm can be computed using the Riesz representation theorem, see for instance [24, 16] for more details. However, evaluating $\|r_{m}(\cdot;y)\|_{\mathbb{V}_{h}^{\prime}}$ for every $y$ in $\mathcal{D}_{s}$ remains unfeasible. In practice, the maximization is performed over a finite dimensional training set $\Theta_{s}\subset\mathcal{D}_{s}$, either chosen sufficiently fine to retain the performance of the algorithm, see for instance [12], or based on a random selection of moderate size [11].

The reduced basis space is constructed offline and gives the following online approximation of $u_{h,k}(s)$

[TABLE]

Remark 1**.**

The approximate solution (3) obtained without using the reduced basis method requires to solve $M_{s}+N_{s}+1$ sparse finite element systems of dimension $N_{h}$. In comparison, the solution of $M_{s}+N_{s}+1$ systems are needed for the approximation (27), each requiring $\mathcal{O}(n^{3})$ operations, where $n$ stands for the dimension of the reduced basis space. The latter is built once for all (offline stage) with a computational cost dominated by the solve of $n$ sparse finite element systems. The value of $n$ depends on the Kolmogorov $n$-width of the solution manifold $\{w_{h}(y)\mathrel{\mathop{\mathchar 58\relax}}\,y\in\mathcal{D}_{s}\}$. We refer for instance to [25, 12] for a detailed complexity analysis of the reduced basis method but note that typically for elliptic problems we have $n<\!\!<N_{h}$. Finally, we anticipate that the proposed reduced basis space is independent on $s$, see Section 3.3.

3.2. Error analysis for a fixed $s$

We now analyze the distortion between $u(s)=(-\Delta)^{-s}f$ and its reduced basis approximation $u_{h,k}^{n}(s)$ given by (27) in the $\mathbb{H}^{r}(\Omega)$ norm. In order to avoid unnecessary technicalities, we assume from now on that $r\in[0,1]$ is chosen such that $r<2s$, which includes the natural choice $r=s$ leading to the energy error and $r=0$ for the $L^{2}(\Omega)$ error. The discussion below can be readily extended to the case $r\geq 2s$ by accounting for the log factor $\ln(h^{-1})$ in the finite element approximation (see Theorem 1). For $u(s)\in\mathbb{H}^{r}(\Omega)$, we decompose the error into three parts

[TABLE]

corresponding to the finite element error, the sinc quadrature error and the reduced basis error, respectively. Given a target tolerance $\varepsilon>0$, we construct a reduced basis space such that (23) holds. In view of Theorems 1 and 2, we select the space discretization and sinc quadrature parameters $h$ and $k$ to balance the finite element and sinc quadrature errors, i.e.,

[TABLE]

for some absolute constants $C_{\textrm{FEM}}$ and $C_{\textrm{SINC}}$.

We now assess the error in the reduced basis modeling by analyzing the behavior of $\|u_{h,k}(s)-u_{h,k}^{n}(s)\|_{\mathbb{H}^{r}(\Omega)}$ as $n$ increases. From the definitions (3) and (27) of $u_{h,k}(s)$ and $u_{h,k}^{n}(s)$, respectively, we have

[TABLE]

Key ingredients in our analysis are estimates for the reduced basis errors $\|w_{h}(y_{l})-w_{h}^{n}(y_{l})\|_{\mathbb{H}^{r}(\Omega)}$ in approximating the inner problems. We discuss this now. Recall that the reduced basis error for the inner problem is given by

[TABLE]

The Kolmogorov $n$-width

[TABLE]

is the benchmark for the best achievable decay. By convention, we set

[TABLE]

It is quite remarkable that the linear space $\mathbb{V}_{h}^{n}$ constructed by the (weak) greedy selection discussed in Section 3.1 leads to an error (31) equivalent to $d_{n}$ [14]. In particular, an exponential decay of the Kolmogorov $n$-width guarantees an exponential decay of the (weak) greedy error. In order to prove that the error in (31) decays exponentially with $n$, see Lemma 3.1 below, we will thus show that the Kolmogorov $n$-width exhibits an exponential decay.

To facilitate the analysis of $d_{n}$, we use the notations in (13) and provide a representation of the finite element functions $w_{h}(y)$ in term of the eigenpairs $\{\mu_{i},\varphi_{i}\}_{i=1}^{N_{h}}\subset\mathbb{R}^{+}\times\mathbb{V}_{h}$, $N_{h}\mathrel{\mathop{\mathchar 58\relax}}=\textrm{dim}(\mathbb{V}_{h})$, of the generalized eigenvalue problem

[TABLE]

Without loss of generality, we assume that the $\varphi_{i}$ are $H_{0}^{1}$-orthonormal, i.e.

[TABLE]

The inverse inequality

[TABLE]

together with the Poincaré inequality (7), yield

[TABLE]

With these notations, we can rewrite $w_{h}(y)$ in (12) as

[TABLE]

We are now in position to assess the reduced basis approximation property.

Lemma 3.1**.**

For any $n\geq 1$ we have

[TABLE]

where $\gamma_{s}$ is given by (26), $C_{1}$ is a constant only depending on $C_{P}$ and

[TABLE]

Proof.

We follow [13] to construct a linear space $\mathbb{W}_{h}^{n}\subset\mathbb{V}_{h}$ with $\dim(\mathbb{W}_{h}^{n})\leq n$ such that for some constants $c_{1}$, $c_{2}$ and $n\geq 1$ we have

[TABLE]

Let $\mathbb{W}_{h}^{n}\mathrel{\mathop{\mathchar 58\relax}}=\operatorname{span}\{w_{h}(y_{1}),\ldots,w_{h}(y_{n})\}$, where the $y_{j}$ are chosen such that the $e^{y_{j}}$ are the transformed Zolotarëv points $[C_{P}^{-2},C_{I}^{2}h^{-2}]$ as in [13], see also [18, 28]. Notice that this interval is dictated by the lower and upper bound of the eigenvalues $\mu_{i}$, see (34). Now, given $y\in\mathcal{D}_{s}$, we define the approximation

[TABLE]

where the coefficients $\alpha_{j}(y)$ are such that

[TABLE]

The above system is a particular rational interpolation problem and has a unique solution according to Lemma 5.13 in [13]. Furthermore, thanks to Lemma 5.17 in [13], we have

[TABLE]

where

[TABLE]

Hence, the error between $w_{h}(y)$ in (35) and $v_{h}^{n}({y})$ in (39) satisfies

[TABLE]

where we have used the $H_{0}^{1}$-orthonormality of the $\{\varphi_{i}\}_{i=1}^{N_{h}}$. With the help of (16), this implies

[TABLE]

The above estimate is (38) with $c_{2}=C^{*}$ and $c_{1}$ only depending on $C_{P}$ and the hidden constant in (40). Moreover, thanks to (16) we also have $d_{0}\leq C_{P}\|f\|$, where $d_{0}$ is defined in (33). Therefore, we have shown that the Kolmogorov $n$-width (see (32) and (33)) satisfies

[TABLE]

To conclude, it remains to relate the error decay of the reduced basis generated by the weak greedy algorithm with parameter $\gamma_{s}$ (see $\eqref{e:gamma_s}$) and the Kolmogorov $n$-width $d_{n}$. Corollary 8.4 in [12], see also Corollary 3.3 in [14], guarantees that (41) implies (36) with $C_{2}=c_{2}/6=C^{*}/6$ and $C_{1}=c_{1}\max(\sqrt{2},\gamma_{s}e^{C_{2}})\lesssim c_{1}$. ∎

Remark 2**.**

We mention that an exponential decay for the reduced basis error for one dimensional parametric problem of the form (12) has already been obtained in [22, 21]. However, the exponential decay is guaranteed for $n\geq n_{\textrm{crit}}$ for some integer $n_{\textrm{crit}}$ depending on the length of the parameter interval $[-M_{s}k,N_{s}k]$. We did not pursue this route as the latter restriction seems prohibitive to take full advantage of the performances of the reduced basis method.

We now use the exponential decay of the reduced basis error for $w_{h}$ obtained in Lemma 3.1 to estimate the error for $u_{h,k}$ defined in (30).

Lemma 3.2**.**

Let $h$ and $k$ be the finite element and sinc quadrature parameters. Let $r\in[0,1]$. For $n\geq 1$, we have

[TABLE]

where $C_{1}$ and $C_{2}(h)$ are the constants in (36).

Proof.

Because $\mathbb{H}^{r}(\Omega)$ are interpolation spaces between $L^{2}(\Omega)$ and $H^{1}_{0}(\Omega)$, see (6), the Poincaré inequality (7) yields

[TABLE]

Therefore, using the estimate (30) for the error and invoking Lemma 3.1, we get

[TABLE]

which is the claimed estimate. ∎

Remark 3**.**

From (37), we see that the constant $C_{2}(h)$ that appears in (36) and (42) tends to [math] as $h$ tends to [math]. In other words, the reduced basis performances deteriorate as the target accuracy $\varepsilon$ tends to [math]. This phenomenon is observed in the numerical experiments reported in Figure 4 of Section 4.

Using Lemma 3.2, we directly derive the following result providing a sufficient condition on the dimension $n$ of the reduced space to achieve an error under a specified tolerance $\delta>0$.

Theorem 3**.**

(offline construction of the reduced space) Let $\varepsilon>0$ be a given tolerance. Assume that $h$ and $k$ are chosen so that (29) holds and that the reduced basis space is constructed such that (23) holds. Then, for any $\delta\geq\varepsilon$ we have

[TABLE]

provided

[TABLE]

where $C$ is a constant only depending on $s$ and $r$ and $\alpha^{*}$ is as in Theorem 1. In particular

[TABLE]

when $\delta=\varepsilon$.

Proof.

The claims directly follow from Lemma 3.2 together with (29). ∎

3.3. Universal Reduced Basis Space

In the previous section we constructed a reduced basis space $\mathbb{V}_{h}^{n}$ to approximate $u_{h,k}(s)$ for a fixed $s\in(0,1)$. We now show that it is possible to take real advantage of the offline work and construct reduced basis spaces approximating the map $s\mapsto u_{h,k}(s)$ for $s\in[s_{\min},s_{\max}]$, with $0<s_{\min}\leq s_{\max}<1$ fixed.

To see this, it suffices to adjust the constant depending on $s$ as follows. First, we let

[TABLE]

Then, we define the domain $\mathcal{D}\mathrel{\mathop{\mathchar 58\relax}}=[-Mk,Nk]$ containing $\mathcal{D}_{s}$ for all $s\in[s_{\min},s_{\max}]$ and, similarly to (26), we introduce the parameter

[TABLE]

Finally, for all $y\in\mathcal{D}$ we approximate $w_{h}(y)$ by $w_{h}^{n}(y)=P_{\mathbb{V}_{h}^{n}}w_{h}(y)$, where the reduced basis space $\mathbb{V}_{h}^{n}$ is constructed as detailed in Section 3.1 upon replacing $M_{s}$, $N_{s}$, $\mathcal{D}_{s}$ and $\gamma_{s}$ by $M$, $N$, $\mathcal{D}$ and $\gamma$, respectively.

With this uniform construction, we directly obtain the universal version of Lemma 3.2 and Theorem 3.

Lemma 3.3**.**

Let $h$ and $k$ be the finite element and sinc quadrature parameters. Let $r\in[0,1]$. For $n\geq 1$ and any $s\in[s_{\min},s_{\max}]$ we have

[TABLE]

where $C_{1}$ and $C_{2}(h)$ are the constants in (36).

Theorem 4**.**

(offline construction of the universal reduced space) Let $\varepsilon>0$ be a given tolerance. Assume that $h$ and $k$ are chosen so that (29) holds and that the reduced basis space is constructed such that (23) holds. Then, for any $\delta\geq\varepsilon$ we have

[TABLE]

provided

[TABLE]

where $C$ is a constant only depending on $s_{\min}$, $s_{\max}$ and $r$, and $\alpha^{*}$ is as in Theorem 1. In particular

[TABLE]

when $\delta=\varepsilon$.

4. Numerical Experiments

We present numerical results to illustrate the performances of the reduced basis approach analyzed in the previous section. Since the focus of this paper is on the reduced basis approximation, the finite element meshsize $h$ and the sinc quadrature parameter $k$ are chosen sufficiently small not to influence the total error, unless otherwise specified. We refer to [5, 4] for an extensive numerical study on the influence of the discretization parameters $h$ and $k$. The space $\mathbb{V}_{h}^{n}$ is built using a weak greedy algorithm on $[-Mk,Nk]$, starting with the snapshot $w_{h}(0)$. Moreover, the training set $\Theta$ consists of $10000$ uniformly distributed points in $[-Mk,Nk]$. Finally, we set $f=1$ in all the numerical examples.

4.1. 1D example

We consider the case $\Omega=(0,1)$. The subdivision $\mathcal{T}_{h}$ of $\Omega$ consists of a uniform partition of $[0,1]$ with subintervals of length $h=2^{-12}$. The sinc quadrature parameter is fixed to $k=0.5$ and the fractional power $s$ varies from $s_{\min}=0.1$ to $s_{\max}=0.9$. In this setting, we have $N=M=395$ and $[-Mk,Nk]=[-197.5,197.5]$.

4.1.1. Reduced basis error for $w_{h}$

We provide in Figure 1-left the evolution of

[TABLE]

versus $n\geq 1$ as indicator of the error

[TABLE]

The observed exponential decay matches the estimate of Lemma 3.1.

Moreover, Figure 1-right reports the values of the selected parameters $y^{n}$ by the weak greedy procedure. We observe that except for $y^{2}$, they are all located in the interval $[0,20]$. This behavior can be in part explained by the estimates provided in Lemma 2.1, indicating the robustness of $w_{h}(y)$ for small values of $y$ and the smallness of $\|w_{h}(y)\|$ for $y$ large. The fact that no negative $y$ is selected is also attributable to the choice of the initial snapshot, namely $w_{h}(0)$, which already provides a good approximation of $w_{h}(y)$ for $y\leq 0$. We mention that similar results are obtained when changing the range for $s$, for instance setting $s_{\min}=0.01$ and $s_{\max}=0.99$.

We comment on the use of an a posteriori error estimate (weak greedy) in place of the true error (greedy), see (25). For this, we compare in Figure 2 the performance of both algorithms and we can conclude that very little efficiency is lost in using the computable a posteriori error estimator.

4.1.2. Reduced basis error for $u_{h,k}$

We now turn our attention to the approximation of $u_{h,k}(s)$ by $u_{h,k}^{n}(s)$. Figure 3 depicts the evolution of

[TABLE]

for various values of fractional power. In agreement with Theorem 4, exponential decay is observed in all cases when using the universal reduced basis space. Notice that in this experiment the sinc quadrature requires $440$ points for $s=0.1,0.9$, $190$ points for $s=0.3,0.7$ and $159$ points for $s=0.5$ to guarantee a sinc quadrature error of the order $e^{-\pi^{2}/k}\approx 2.7\times 10^{-9}$ for $k=0.5$. A comparable reduced basis accuracy in the $L^{2}(\Omega)$ norm is achieved for only $n=20$ for $0.5\leq s\leq 0.9$.

Finally, we study numerically the behavior of the constant $C_{2}(h)$ given by (37). We set $s=0.1$ and consider a sequence of uniform partitions of $[0,1]$ with subintervals of length $h=2^{-k}$ for $k=6,8,10,12$. The reduced basis error $e_{u(0.1)}(n)$ versus $n$ is reported in Figure 4 for each finite element discretization. As predicted by Theorem 4, the exponential decay encoded in $C_{2}=C_{2}(h)$ deteriorates as $h\to 0$. The $L^{2}(\Omega)$ norm of the error behaves like $e^{-1.4n}$ for $h=2^{-6}$ and $e^{-0.7n}$ for $h=2^{-12}$.

4.2. 2D examples

We now consider two dimensional domains: a square domain $\Omega=(0,1)^{2}$ and an L-shaped domain $\Omega=(0,1)^{2}\setminus([0,0.5]\times[0.5,1])$. The space discretization consists of a Delaunay triangulation with $22968$ elements for $\Omega=(0,1)^{2}$ and $17190$ elements for $\Omega=(0,1)^{2}\setminus([0,0.5]\times[0.5,1])$. In both cases, the elements in the triangulation have diameters between $0.005$ and $0.01$. All the other parameters are the same as in Section 4.1.

The evolution of the reduced basis error $e_{u(s)}(n)$ for different values of $s$ is reported in Figure 5. As for the one dimensional case, exponential decays are observed for all values of $s\in[0.1,0.9]$ and irrespectively of the shape of the domain.