Analysis of an approximation to a fractional extension problem

TL;DR

This paper investigates an approximation method for a fractional extension problem related to Bessel-type operators, demonstrating its consistency, stability, and convergence through theoretical analysis and numerical experiments.

Contribution

It introduces a novel approximation approach for fractional extension problems, extending the framework of time-stepping methods to abstract Bessel-type problems.

Findings

The method is consistent, stable, and convergent.

Numerical experiments confirm theoretical predictions.

The approach generalizes existing techniques for fractional operators.

Abstract

The purpose of this work is to study an approximation to an abstract Bessel-type problem, which is a generalization of the extension problem associated with fractional powers of the Laplace operator. Motivated by the success of such approaches in the analysis of time-stepping methods for abstract Cauchy problems, we adopt a similar framework, herein. The proposed method differs from many standard techniques, as we approximate the true solution to the abstract problem, rather than solve an associated discrete problem. The numerical method is shown to be consistent, stable, and convergent in an appropriate Banach space. These results are built upon well understood results from semigroup theory. Numerical experiments are provided to demonstrate the theoretical results.

Figures

Equations

$$(-\Delta {)}^{s}v=f,x\in {\mathbb{R}}^d,$$

$$(-\Delta {)}^{s}v(x):=C_{d,s}{\int }_{{\mathbb{R}}^d}\frac{v(x)-v(z)}{|x-z{|}^{d+2s}}dz$$

$$\left\{\begin{array}[]{ll}\partial_{t}^{2}u(x,t)+\frac{1-2s}{t}\partial_{t}u(x,t)=-\Delta_{x}u(x,t),&\quad(x,t)\in{\mathbb{R}}^{d}\times{\mathbb{R}}_{+},\\ u(x,0)=v(x),&\quad x\in{\mathbb{R}}^{d},\end{array}\right.$$

$$c_{s}f(x)=c_s(-\Delta {)}^{s}v(x)=-\underset{t\to 0^{+}}{\lim }{t}^{1-2s}{\partial }_{t}u(x,t),x\in {\mathbb{R}}^d,$$

$$\left\{\begin{array}[]{ll}u^{\prime\prime}(t)+\frac{\alpha}{t}u^{\prime}(t)=-Au(t),&\quad t\in(0,\infty),\\ u(0)=u_{0}\in X,\end{array}\right.$$

$$\left\{\begin{array}[]{ll}u^{\prime\prime}(t)+\frac{\alpha}{t}u^{\prime}(t)=-Au(t),&\quad t\in(0,\infty),\\ u(0)=u_{0}\in X,&\end{array}\right.$$

$$\Vert (I-tA)u\Vert >\Vert u\Vert ,\text{mbox}forallu\in D(A),t\ge 0,$$

$$R(I-tA)=X,\text{mbox}forallt\ge 0.$$

$$[u,v]:=\Vert v\Vert \underset{\epsilon \to 0^{-}}{\lim }\frac{1}{\epsilon }\left(\Vert v+\epsilon u\Vert -\Vert v\Vert \right),u,v\in X.$$

$$\text{mbox}Re[Au,u]<0,\text{mbox}forallu\in D(A).$$

$$-\Vert u\Vert \Vert v\Vert \le [u,v]\le \Vert u\Vert \Vert v\Vert ,$$

$$[u,v]={\int }_{\Omega }uvd\Omega \text{mbox}and\Vert u{\Vert }^2=[u,u].$$

$$[Au,u]=-[\nabla u,\nabla u]\le -C[u,u],$$

$$\mu (A):=\underset{u\in D(A)}{\sup }\frac{[Au,u]}{\Vert u{\Vert }^2}.$$

$$\Vert T(t)\Vert \le e^{t\mu (A)},t\ge 0.$$

$$\left\{\begin{array}[]{ll}x^{\prime}(t)=Ax(t),&\quad t>0,\\ x(0)=x_{0}.\end{array}\right.$$

$$D_t^{+}x(t):=\underset{\epsilon \to 0^{+}}{\mathrm{limsup}}\frac{x(t+\epsilon )-x(t)}{\epsilon }.$$

$$D_t^{+}\Vert x(t)\Vert =\frac{[Ax,x]}{\Vert x{\Vert }^2}\Vert x\Vert$$

$$\rho (A):=\left\{\lambda \in \mathbb{C}:(\lambda -A)\text{mbox}isinjectiveand(\lambda -A{)}^{-1}\in L(X)\right\}$$

$$J_{-A}^{s}u:=\frac{1}{\Gamma (s)\Gamma (1-s)}{\int }_0^{\infty }{z}^{s-1}(z-A{)}^{-1}(-A)udz.$$

$$D(A)=\{u\in \overline{D(A)}:\exists {\delta }_0>0,\forall \delta \in (0,{\delta }_0)$$

$$\text{mbox}s.t.u\in D((-A+\delta {)}^s),\underset{\delta \to 0^{+}}{\lim }(-A+\delta {)}^{s}u\text{mbox}exists\}$$

$$\overline{J_{-A}^s}=(-A{)}^s.$$

$$u(t)=\frac{1}{\Gamma (s)}{\int }_0^{\infty }{z}^{s-1}{e}^{-t^2/4z}T(z)(-A{)}^s{u}_{0}dz$$

$$\underset{t\to 0^{+}}{\lim }{t}^{1-2s}{u}^{\prime }(t)=c_s(-A{)}^s{u}_0,$$

$$u^{\prime }(t)=\frac{1}{\Gamma (s)}{\int }_0^{\infty }\left[\frac{-t}{2z}\right]z^{s-1}{e}^{-t^2/4z}T(z)(-A{)}^s{u}_{0}dz$$

$$u^{\prime \prime }(t)=\frac{1}{\Gamma (s)}{\int }_0^{\infty }\left[\frac{-1}{2z}+\frac{t^2}{4z^2}\right]z^{s-1}{e}^{-t^2/4z}T(z)(-A{)}^s{u}_{0}dz.$$

$$u^{\prime \prime }(t)+\frac{\alpha }{t}{u}^{\prime }(t)$$

$$t^{1-2s}{u}^{\prime }(t)$$

$${\int }_0^{\infty }{z}^s\left(\frac{d}{dz}{e}^{-t^2/4z}\right)T(z)(-A{)}^s{u}_{0}dz$$

Full text

1 Introduction

Due to its wide array of applications, the computation of fractional powers of the Laplacian, and other elliptic operators, has become a problem of great interest [2, 9, 25, 26]. Fractional powers of operators have received this attention due to their applicability to the accurate modeling of real-world problems with varying scales. Examples of such problems are found in porous media flow, peridynamics, nonlocal continuum field theory, finance, and many others (see [8, 10, 13, 29], and the references therein). However, due to the nonlocal nature of fractional problems, the development of accurate and efficient computational algorithms has been greatly hindered. In response to this hurdle, the revolutionary work by Caffarelli and Silvestre demonstrated that the nonlocal fractional Laplacian problem may be recast into an equivalent local problem which is amenable to certain standard numerical techniques [7].

In order to demonstrate, fix $\displaystyle 0<s<1$ and consider the function $\displaystyle v\,:\,\mathbb{R}^{d}\to\mathbb{R}$ which solves

[TABLE]

where $\displaystyle f$ is a given function of appropriate regularity and the fractional power of the Laplacian is defined via the hypersingular integral

[TABLE]

with $\displaystyle C_{d,s}$ being some normalization constant [7]. Caffarelli and Silvestre then showed that for $\displaystyle 0<s<1,$ the fractional Laplacian in $\displaystyle(\ref{aaa1})$ can be realized via the Dirichlet-to-Neumann map for a function $\displaystyle u\,:\,{\mathbb{R}}^{d}\times{\mathbb{R}}_{+}\to{\mathbb{R}}.$ The function $\displaystyle u$ then satisfies the following Bessel-type equation

[TABLE]

where $\displaystyle v$ is the solution to $\displaystyle(\ref{aaa1}).$ Moreover, the function $\displaystyle u$ satisfying $\displaystyle(\ref{e1})$ may be viewed as the harmonic extension of $\displaystyle v$ to the fractional dimension $\displaystyle 2-2s.$ One can then calculate $\displaystyle(-\Delta)^{s}v$ as

[TABLE]

where $\displaystyle c_{s}$ is a constant depending upon $\displaystyle s,$ but not on the dimension $\displaystyle d,$ and $\displaystyle u$ solves $\displaystyle(\ref{e1}).$ Moreover, the final equality is $\displaystyle(\ref{e2})$ is known as the conormal derivative. We have emphasized the relationship between $\displaystyle v$ and $\displaystyle f$ in $\displaystyle(\ref{e2})$ in order to clarify the ensuing exposition. Further details and applications of this method may be found in [3, 7, 14, 31].

Several numerical algorithms for solving problems involving fractional powers of the Laplacian have been developed based upon the Caffarelli-Silvestre extension technique [12, 18, 25, 26]. These methods have also been extended to problems involving general elliptic operators and various other physically relevant nonlocal problems. Despite the increased interest in such methods, these explorations are still in their infancy and merit further consideration. Moreover, there is the need to improve the efficiency and accuracy of existing numerical methods.

In order to develop more robust approximations to $\displaystyle(\ref{e1}),$ we consider the following abstract Bessel-type problem

[TABLE]

where $\displaystyle X$ is a real separable Banach space, $\displaystyle A\,:\,D(A)\subseteq X\to X$ is a densely defined linear operator, and $\displaystyle\alpha\mathrel{\mathop{:}}=1-2s,\ s\in(0,1).$ There have been numerous studies concerned with approximating particular forms of $\displaystyle(\ref{e3}),$ but herein we present a novel method and demonstrate its properties in an arbitrary Banach space for appropriate abstract operators $\displaystyle A.$ In particular, such analysis allows for mesh independent results, which are quite useful in application.

Analysis and representations of the true solution of $\displaystyle(\ref{e3})$ have been studied in multiple settings [3, 14, 31]. However, to the author’s knowledge, approximations to $\displaystyle(\ref{e3})$ in this abstract setting have yet to be explored. In particular, the true solution to $\displaystyle(\ref{e3})$ involves the abstract semigroup generated by the operator $\displaystyle A$ and there have been few methods employing direct approximations of this operator. These approximations preserve numerous desirable properties inherent to the original semigroup operator, often resulting in superior approximations, resulting in qualitatively superior approximations [6, 15, 27, 28].

This article is organized as follows. In Section 2 relevant mathematical preliminaries are outlined for the reader’s convenience. Section 3 is concerned with the development and stability of the proposed approximation to $\displaystyle(\ref{e1}).$ Section 4 provides the necessary results to guarantee convergence of the proposed method. Section 5 presents some concrete example of our abstract problem and provides numerical experiments verifying the analytic result. Finally, concluding remarks and an outline of future endeavors are given in Section 6.

2 Preliminaries

Let $\displaystyle X$ be a real separable Banach space with norm denoted $\displaystyle\|\cdot\|$ and let $\displaystyle B(X)$ denote the Banach algebra of bounded operators on $\displaystyle(X,\|\cdot\|).$ We wish to consider the following Bessel-type problem

[TABLE]

where $\displaystyle\alpha\mathrel{\mathop{:}}=1-2s,\ s\in(0,1),$ and $\displaystyle A\,:\,D(A)\subseteq X\to X$ is linear. The operator norm on $\displaystyle X$ will also be denoted by $\displaystyle\|\cdot\|.$ Before providing details regarding a solution to $\displaystyle(\ref{b1})$ we detail a necessary assumption and some relevant background.

Assumption 1**.**

Let the operator $\displaystyle A$ be closed and densely defined in $\displaystyle X.$ Moreover, let $\displaystyle A$ be strictly $\displaystyle m-$dissipative in $\displaystyle X.$

We recall that an operator $\displaystyle A$ is called strictly dissipative if and only if

[TABLE]

and it is called strictly $\displaystyle m-$dissipative if and only if it satisfies the additional range condition

[TABLE]

The notion of a dissipative operator may also be formalized in the following manner. Define the (left) semi-inner-product $\displaystyle[\cdot,\cdot]\,:\,X\times X\to{\mathbb{R}}$ on the Banach space $\displaystyle X$ as

[TABLE]

Then $\displaystyle A$ is strictly dissipative if and only if

[TABLE]

Note that the limit in $\displaystyle(\ref{semi-inner})$ exists as $\displaystyle[\cdot,\cdot]$ is a G-differential of the norm $\displaystyle\|\cdot\|.$ If $\displaystyle X$ is a Hilbert space then we have that $\displaystyle(\ref{semi-inner})$ coincides with the true inner product on $\displaystyle X,$ that is $\displaystyle[\cdot,\cdot]=\langle\cdot,\cdot\rangle.$ Moreover, we have that $\displaystyle[u,u]=\|u\|^{2},$ for any $\displaystyle u\in X,$ and

[TABLE]

for any $\displaystyle u,v\in X.$ More details regarding properties of $\displaystyle(\ref{semi-inner})$ may be found in [11].

Remark 2.1*.*

It is worth noting that many existing works in the direction of abstract numerical analysis often require that an operator simply be m-dissipative. The current work requires a slightly stricter result in order to guarantee the appropriate decay of the solution (a requirement similar to that in [14, 25]).

From this, it readily follows that the resolvent operator $\displaystyle R(t,A)\mathrel{\mathop{:}}=(I-tA)^{-1}\,:\,X\to X$ of a strictly $\displaystyle m-$dissipative operator is well-defined and a contraction in $\displaystyle X.$ It also follows that strictly $\displaystyle m-$dissipative operators generate strongly continuous semigroups of contractions [16, 21]. We will denote the semigroup generated by $\displaystyle A$ by $\displaystyle(T(t))_{t\geq 0}.$ Note that if $\displaystyle X$ is reflexive then every strictly $\displaystyle m-$dissipative operator is densely defined. More details regarding $\displaystyle m-$dissipative operators and their properties can be found in [4, 21].

The following example demonstrates that standard elliptic operators are examples of the proposed abstract formulation.

Example 2.1**.**

Let $\displaystyle\Omega\subset{\mathbb{R}}^{d}$ be a bounded domain with smooth boundary. We may then consider the operator $\displaystyle A=\Delta$ with homogeneous Dirichlet boundary conditions, on $\displaystyle\partial\Omega.$ Let $\displaystyle X=H_{0}^{1}(\Omega)\cap H^{2}(\Omega),$ and for $\displaystyle u,v\in X,$ we define

[TABLE]

We then have, by employing integration by parts and standard Sobolev-type inequalities, that

[TABLE]

where $\displaystyle C$ is some positive constant. Thus, we have that $\displaystyle\mbox{Re}[Au,u]<0$ for all $\displaystyle u\in X.$ The range condition follows from standard results, as well.

We also define the logarithmic norm of an operator $\displaystyle A$ with respect to the semi-inner-product $\displaystyle[\cdot,\cdot]$ given by $\displaystyle(\ref{semi-inner})$ to be

[TABLE]

Logarithmic norms have been well studied and more details regarding their properties may be found in [30].

Example 2.2**.**

If we consider the situation in Example 2.1, again, then it is easy to see that $\displaystyle\mu(A)<0,$ where $\displaystyle A=\Delta.$

Remark 2.2*.*

It is worth noting that one could easily have considered (the negation) of any strongly elliptic operator in Examples 2.1 and 2.2, not just the Laplacian. In fact, one could equivalently define strongly elliptic operators as any operator $\displaystyle L$ such that $\displaystyle\mu(L)>0.$

Remark 2.3*.*

While not the focus of this work, one may easily use $\displaystyle(\ref{lognorm})$ to define the notion of ellipticity of nonlinear operators on $\displaystyle X.$ This approach can be quite useful and allows for the consideration of more generalized notions of numerical stability for stiff systems of equations.

The above formulations allow us to develop the following, very useful, result.

Lemma 2.1**.**

Let $\displaystyle A\,:\,D(A)\subseteq X\to X$ be the generator of a strongly continuous semigroup $\displaystyle(T(t))_{t\geq 0}$ in $\displaystyle D(A).$ Then,

[TABLE]

Proof.

This result has been proven elsewhere, but we include a proof for completeness. Let $\displaystyle A$ be the generator of the semigroup $\displaystyle(T(t))_{t\geq 0}.$ Then it follows that $\displaystyle x(t)=T(t)x_{0}$ is the solution to

[TABLE]

Recall that the upper right Dini derivative of a function $\displaystyle x(t)$ is defined as

[TABLE]

By direct calculation, we have

[TABLE]

(see, for instance, [11]). Employing $\displaystyle(\ref{lognorm})$ and solving $\displaystyle(\ref{log2})$ directly yields the desired result.

The remainder of this section will serve as a brief introduction to fractional powers of operators in Banach spaces and provide some existence results for problem $\displaystyle(\ref{b1})$. These introductory notes are based primarily on the reference [22], with some influence from [23]. Let $\displaystyle L(X)$ denote the set of all bounded linear operators mapping $\displaystyle X$ into $\displaystyle X$ and

[TABLE]

denote the resolvent set of the linear operator $\displaystyle A.$ For operators $\displaystyle A\,:\,D(A)\subseteq X\to X$ satisfying Assumption 2.1 and $\displaystyle u\in D(A),$ we define the Balakrishnan operator with power $\displaystyle s\in(0,1)$ and base $\displaystyle-A$ as

[TABLE]

This operator will be used to define fractional powers of the operator $\displaystyle A.$ For more details see [22].

Remark 2.4*.*

The operator $\displaystyle J_{-A}^{s}$ inherits numerous desirable properties from the underlying operator $\displaystyle A.$ For instance, if $\displaystyle A$ is bounded or injective, then so is $\displaystyle J_{-A}^{s}$ [22].

With $\displaystyle(\ref{bala}),$ we are now able to define fractional powers of abstract operators.

Definition 2.1**.**

Let $\displaystyle A$ be as in Assumption 2.1 and fix $\displaystyle s\in(0,1).$ Then we define $\displaystyle(-A)^{s}$ as

• i.

$\displaystyle(-A)^{s}\mathrel{\mathop{:}}=J_{-A}^{s},$ for $\displaystyle A\in L(X);$

• ii.

$\displaystyle(-A)^{s}\mathrel{\mathop{:}}=\left(J_{(-A)^{-1}}^{s}\right)^{-1},$ for $\displaystyle A$ being unbounded.

• iii.

$\displaystyle(-A)^{s}\mathrel{\mathop{:}}=\lim_{\delta\to 0^{+}}(-A+\delta)^{s}u,$ for $\displaystyle A$ being unbounded, $\displaystyle 0\in\sigma(A)$ (where $\displaystyle\sigma(A)$ is the spectrum of $\displaystyle A$), and $\displaystyle D((-A)^{s})$ given by

[TABLE]

It is worth noting that the above definition yields a well-defined closed operator which extends the original Balakrishnan operator. Moreover, in [23], it was shown that for a densely defined linear operator $\displaystyle A,$ we have

[TABLE]

There is still a wealth of information that one could present to provide background on fractional powers of operators, but these ideas will not be pertinent to the study at hand. Interested readers should see [22] for more details.

Based on Assumption 2.1, we are in fact considering a special case of that studied in [14]. Though the following result was included in [14], the proof was omitted for this special case, hence, we include a proof for completeness. The methods employed below are similar to those in [23].

Theorem 2.1**.**

Let $\displaystyle u_{0}\in D((-A)^{s}).$ Then a solution to $\displaystyle(\ref{b1})$ is given by

[TABLE]

and also satisfies

[TABLE]

where $\displaystyle c_{s}\mathrel{\mathop{:}}=2^{1-2s}\Gamma(1-s)/\Gamma(s).$

Proof.

We first show $\displaystyle(\ref{usol1})$ by demonstrating directly that $\displaystyle(\ref{usol1})$ satisfies $\displaystyle(\ref{b1})$ for $\displaystyle t>0.$ To that end we have

[TABLE]

and

[TABLE]

Thus, for $\displaystyle t>0$ we have

[TABLE]

By noting that $\displaystyle\textstyle\frac{d}{dz}T(z)(-A)^{s}u_{0}=AT(z)(-A)^{s}u_{0}$ and employing integration by parts, we have

[TABLE]

which shows that $\displaystyle(\ref{usol1})$ is a solution to $\displaystyle(\ref{b1}).$

Next we show $\displaystyle(\ref{uder})$ holds for $\displaystyle u_{0}\in D((-A)^{s}).$ Note that

[TABLE]

By employing integration by parts in $\displaystyle(\ref{lim1b})$ we obtain

[TABLE]

where the boundary terms in $\displaystyle(\ref{lim1c})$ go to zero due to Assumption 2.1. After substituting $\displaystyle(\ref{lim1c})$ into $\displaystyle(\ref{lim1b})$ and employing integration by parts yet again, we obtain

[TABLE]

Note that the first integral in $\displaystyle(\ref{lim2})$ is Bochner integrable and is in fact the desired result, as

[TABLE]

by Definition 2.1 and $\displaystyle(\ref{bala1}).$ Thus, $\displaystyle(\ref{uder})$ follows if the second integral in $\displaystyle(\ref{lim2})$ goes to zero. By Taylor expanding T(z) about $\displaystyle z=0$ and employing Assumption 2.1, we have

[TABLE]

Since $\displaystyle D(A)\subseteq D((-A)^{s}),$ it follows that $\displaystyle\|Au_{0}\|$ is bounded by assumption and we have

[TABLE]

Combining $\displaystyle(\ref{lim3})$ and $\displaystyle(\ref{lim5})$ yields the desired result.

Lemma 2.2**.**

For all $\displaystyle u_{0}\in X,$ a solution to $\displaystyle(\ref{b1})$ is given by $\displaystyle(\ref{usol1}).$

Proof.

This result follows from [23] and the fact that $\displaystyle A$ is densely defined in $\displaystyle X.$

It is worth noting that the solution to $\displaystyle(\ref{b1})$ is not known to be unique in the current setting. However, there are numerous results regarding situations in which uniqueness is known, for instance, when $\displaystyle A$ is the Laplacian with Dirichlet boundary conditions, when $\displaystyle A$ is a sectorial operator and $\displaystyle X$ is a Hilbert space, and for all operators $\displaystyle A$ when $\displaystyle s=1/2$ [3, 14, 23, 31]. In order to avoid these issues in our ensuing analysis, we introduce the following assumption.

Assumption 2**.**

The solution to $\displaystyle(\ref{b1})$ is unique in $\displaystyle X$ for all $\displaystyle u_{0}\in D((-A)^{s}),$ and hence, given by $\displaystyle(\ref{usol1}).$

Remark 2.5*.*

One could replace Assumption 2.2 with the condition that the solution to $\displaystyle(\ref{b1})$ decays to zero at infinity, in some weak sense, and obtained the same result.

Remark 2.6*.*

It is worth noting that $\displaystyle(\ref{usol1})$ may seem a bit strange as it involves the fractional power of the operator $\displaystyle-A.$ However, as discussed in the Introduction, when considering $\displaystyle(\ref{e1})$ or $\displaystyle(\ref{b1})$ in applications the function $\displaystyle(-A)^{s}u_{0}$ will be given.

We close this section by outlining some basic properties regarding the semigroup generated by $\displaystyle A$ and the operator $\displaystyle(-A)^{s}.$ For more details, see [17].

Lemma 2.3**.**

Let $\displaystyle\alpha\geq 0$ and $\displaystyle 0\leq\gamma\leq 1.$ Then there exists a constant $\displaystyle C>0$ such that

• i.

$\displaystyle\|(-A)^{\alpha}T(t)\|\leq Ct^{-\alpha},$* for $\displaystyle t>0;$*

• ii.

$\displaystyle(-A)^{\alpha}T(t)=T(t)(-A)^{\alpha},$* on $\displaystyle D((-A)^{\alpha});$*

• iii.

If $\displaystyle\alpha\geq\gamma,$ then $\displaystyle D((-A)^{\alpha}\subseteq D((-A)^{\gamma}).$

3 Derivation of Approximation

The properties of various forms of $\displaystyle(\ref{usol1})$ have been studied by several authors, see [3, 14, 23, 31] and the references therein. However, in this paper, we are concerned with the development and analysis of an approximation to $\displaystyle(\ref{usol1}).$ To that end, fix $\displaystyle 0\ll M<N<\infty$ with $\displaystyle M,N\in{\mathbb{N}}.$ Then we define $\displaystyle h\mathrel{\mathop{:}}=M/(N+1)$ and set $\displaystyle z_{k}\mathrel{\mathop{:}}=kh,\ k=0,\ldots,N+1.$ For the sake of the ensuing analysis, we define

[TABLE]

We then have the following result.

Lemma 3.1**.**

Let $\displaystyle u$ and $\displaystyle u^{(M)}$ be defined as in $\displaystyle(\ref{usol1})$ and $\displaystyle(\ref{a1}),$ respectively. Then, if $\displaystyle u_{0}\in D((-A)^{s}),$ we have

[TABLE]

where $\displaystyle C$ is a constant independent of $\displaystyle s$ and $\displaystyle M.$

Proof.

We proceed by direct calculation. We have

[TABLE]

where $\displaystyle\mu(A)<0$ by Assumption 2.1 and is defined in $\displaystyle(\ref{lognorm}).$ It now remains to bound the integral in $\displaystyle(\ref{lemma1a}).$ To that end, let $\displaystyle w=-\mu(A)z,$ to obtain

[TABLE]

Employing the substitution $\displaystyle y=w+\mu(A)M$ in $\displaystyle(\ref{1int1})$ then yields

[TABLE]

where

[TABLE]

Since $\displaystyle s\in(0,1),$ we have

[TABLE]

By substituting $\displaystyle(\ref{1int1})$-$\displaystyle(\ref{1int3})$ into $\displaystyle(\ref{lemma1a}),$ we obtain

[TABLE]

Remark 3.1*.*

While the constant, $\displaystyle C,$ in Lemma 3.1 is independent of $\displaystyle s,$ $\displaystyle t,$ and $\displaystyle M,$ it does depend on $\displaystyle\mu(A).$ Hence, in practice, if $\displaystyle\mu(A)$ is close to zero we will need to increase the value of $\displaystyle M$ to obtain the desired accuracy.

Now that it has been established that the truncation of the integral will not have drastic effects on the approximation of $\displaystyle(\ref{usol1}),$ for large $\displaystyle M,$ we present the proposed approximation to $\displaystyle u(t)$ in $\displaystyle X.$ For $\displaystyle u_{0}\in D((-A)^{s}),$ we define the following continuous in time approximation to $\displaystyle(\ref{usol1})$

[TABLE]

where

[TABLE]

with $\displaystyle z_{k+1/2}\mathrel{\mathop{:}}=z_{k}+h/2.$

Assumption 3**.**

Let the operator $\displaystyle S(w)$ be nonexpansive in $\displaystyle X.$

Remark 3.2*.*

Note that if $\displaystyle X$ is in fact a Hilbert space, then the assumption that $\displaystyle S(w)$ is nonexpansive will always be valid due to Assumption 2.1. That is, let $\displaystyle v\in D(A),$ then

[TABLE]

where $\displaystyle[\cdot,\cdot]$ is the inner product on the Hilbert space $\displaystyle X.$ However, in general Assumption 3.1 does not hold in general Banach spaces, $\displaystyle X,$ with $\displaystyle A$ satisfying Assumption 2.1. For more details see [4].

We begin by demonstrating that $\displaystyle(\ref{approx})$ is stable in $\displaystyle X.$

Lemma 3.2**.**

Let $\displaystyle u_{0}\in D((-A)^{s})$ and $\displaystyle v(t)$ be as defined in $\displaystyle(\ref{approx}).$ Then we have

[TABLE]

where $\displaystyle C$ is a constant independent of $\displaystyle h$ and $\displaystyle k=0,\ldots,N.$

Proof.

Direct calculation yields

[TABLE]

by Assumption 3.1. Thus it we need only consider the $\displaystyle\gamma_{k}(t).$ We first note that since $\displaystyle\gamma_{k}(t)\geq 0$ for $\displaystyle t\geq 0,$ we have

[TABLE]

for some $\displaystyle\xi_{k}\in(z_{k},z_{k+1}),\ k=0,...,N.$ Now let $\displaystyle\beta\in{\mathbb{N}}$ be the index such that $\displaystyle z_{\beta}\leq 1<z_{\beta+1}.$ We then have, by $\displaystyle(\ref{new1}),$

[TABLE]

By noting that

[TABLE]

we then obtain

[TABLE]

where $\displaystyle C$ is a constant independent of $\displaystyle h$ and $\displaystyle k=0,\ldots,N.$ Substituting $\displaystyle(\ref{st1e})$ into $\displaystyle(\ref{st1})$ completes the proof.

4 Convergence

We now present several convergence estimates that will be necessary for the main result of this section.

Lemma 4.1**.**

Let $\displaystyle u_{0}\in D((-A)^{\ell+s}),\ \ell=1,2.$ Then we have

[TABLE]

for $\displaystyle t\in[0,\infty),$ where $\displaystyle C$ is a constant independent of $\displaystyle h$ and $\displaystyle k=0,...,N.$

Proof.

We begin by defining

[TABLE]

for all $\displaystyle t\in[0,\infty).$ Thus, by employing $\displaystyle(\ref{ffff}),$ we have by direct calculation

[TABLE]

for some $\displaystyle\xi_{k}\in(z_{k},z_{k+1}).$ By noting that

[TABLE]

it is clear from $\displaystyle(\ref{ffff1})$ that it only remains to show

[TABLE]

for $\displaystyle\ell=1,2.$ We first note that, by Hille’s theorem, we have

[TABLE]

From $\displaystyle(\ref{ffff2}),$ Assumption 2.1, and Lemma 2.3, we have

[TABLE]

since $\displaystyle z^{-\eta}e^{-t^{2}/4z}$ is bounded for $\displaystyle z\in(0,\infty)$ and $\displaystyle\eta=2,3,4.$

We next note that

[TABLE]

Combining $\displaystyle(\ref{fffff})$ and $\displaystyle(\ref{ffff3})$ into $\displaystyle(\ref{ffff1})$ then yields

[TABLE]

where $\displaystyle\ell=1,2$ for $\displaystyle u_{0}\in D((-A)^{\ell+s}),$ respectively, and $\displaystyle C$ is a constant independent of $\displaystyle h$ and $\displaystyle k=0,...,N.$

We now define a family of operators which will be employed in the proof of the following lemma. Let $\displaystyle\Lambda_{k}\,:\,[0,\infty)\to B(X)$ be defined as

[TABLE]

with $\displaystyle\Lambda_{0}(z)\mathrel{\mathop{:}}=T(z).$ The $\displaystyle\Lambda_{k}$ satisfy the following recurrence relation

[TABLE]

Lemma 4.2**.**

Let $\displaystyle u_{0}\in D((-A)^{\ell+s}),\ \ell=1,2.$ Then we have

[TABLE]

where $\displaystyle C$ is a constant independent of $\displaystyle h$ and $\displaystyle k=0,\ldots,N.$

Proof.

By employing $\displaystyle(\ref{lam1})$ and $\displaystyle(\ref{lam2}),$ we readily obtain

[TABLE]

for $\displaystyle u_{0}\in D((-A)^{2+s}).$ Further, consideration of $\displaystyle(\ref{semi1})$ gives

[TABLE]

Combining $\displaystyle(\ref{semi1})$ and $\displaystyle(\ref{semi1a})$ finally gives

[TABLE]

The result for $\displaystyle\ell=1$ then follows by taking the norm of both sides of $\displaystyle(\ref{semi1c}).$ Since $\displaystyle A(I-A)^{-1}$ is a bounded operator on $\displaystyle X,$ it follows that $\displaystyle(\ref{semi1c})$ holds for $\displaystyle u_{0}\in D(A^{1+s}),$ and taking the norm of both sides yields the result for $\displaystyle\ell=2.$

Combining the above results, we have the following theorem demonstrating the desired convergence rate of the proposed scheme.

Theorem 4.1**.**

Let $\displaystyle u_{0}\in D((-A)^{\ell+s}),\ \ell=1,2.$ Then, we have

[TABLE]

where $\displaystyle C$ is a constant independent of $\displaystyle h$ and $\displaystyle k=0,\ldots,N.$

Proof.

We begin by writing the difference as follows

[TABLE]

We will now apply the previous lemmas to $\displaystyle I_{1},\ I_{2},$ and $\displaystyle I_{3}.$ First, by Lemma 3.1, it is immediately clear that we have

[TABLE]

and thus $\displaystyle(\ref{est1})$ may be made arbitrarily small for large enough $\displaystyle M.$ Next, by Lemma 4.1 we have

[TABLE]

where the value of $\displaystyle C$ changes throughout the calculations, but remains independent of $\displaystyle h$ and $\displaystyle k=0,...,N,$ with $\displaystyle\ell=1$ or $\displaystyle\ell=2,$ depending on the regularity of the initial data. Finally, by Lemma 4.2, Assumption 2.1, and Assumption 3.1, we have

[TABLE]

where, once again, $\displaystyle C$ is independent of $\displaystyle h$ and $\displaystyle k=0,...,N,$ with $\displaystyle\ell$ as before. Combining the obtained bounds in $\displaystyle(\ref{est1})$-$\displaystyle(\ref{est3})$ yields the desired result.

Remark 4.1*.*

The first term in our estimate

[TABLE]

is an expected and standard term (a similar term appears in, for example, [5, 9]). However, we intend to remove this estimate in future work by relating the choice of $\displaystyle M$ to the parameter $\displaystyle h.$

5 Numerical Experiments

In order to further verify the sharpness of our analysis and demonstrate the efficacy of the proposed method, we present some numerical experiments. The first and second examples consider a one-dimensional and two-dimensional problem, respectively, where the given operator is the Laplace operator. In both cases, we see that the expected theoretical convergence results are recovered. The third example outlines how our theory easily includes general strongly elliptic differential operators which still satisfy the theoretical results from Theorem 4.1.

It is the case that the nonlocality of fractional problems can increase memory requirements to the point that some standard methods become quite inefficient. To that end, we briefly explore the memory requirements needed to implement the proposed method $\displaystyle(\ref{approx}).$ If we fix a particular numerical approximation of the operator $\displaystyle A,$ then our method can be shown to be $\displaystyle\mathcal{O}(N).$ Initially, it may seem that the computations are more involved, due to the evaluation of the operator $\displaystyle S^{k}(h)$ (seeming to require $\displaystyle N$ additions and $\displaystyle(N^{2}+N)/2$ multiplications), but these evaluations can be be completed with only $\displaystyle N$ additions and $\displaystyle N$ multiplications by employing nested operations similar to the classical Horner algorithm for polynomials. To determine an exact efficiency, one must consider the precise discretization method employed to approximate the operator $\displaystyle A.$ Due to the inversion needed to compute $\displaystyle S,$ it is clear that discretizations resulting in a banded matrix structure will prove the most efficient. We leave more detailed studies of the efficiency, as well as computational run time trials, for future work.

5.1 A One-Dimensional Example

Let $\displaystyle\Omega\mathrel{\mathop{:}}=(-1,1)\subset{\mathbb{R}}$ and define $\displaystyle X=L^{2}(\Omega)$ with its usual inner product and norm. Consider the problem

[TABLE]

where

[TABLE]

and the fractional Laplacian is defined via the spectral definition (which can be done due to the smoothness of $\displaystyle f$). That is,

[TABLE]

where the set $\displaystyle\{\varphi_{k}\}_{k\in{\mathbb{N}}}\subset H_{0}^{1}(\Omega)$ forms an orthonormal basis of $\displaystyle X$ and satisfy

[TABLE]

for each $\displaystyle k\in{\mathbb{N}},$ and

[TABLE]

where $\displaystyle[\cdot,\cdot]$ is the usual inner product on $\displaystyle X.$ The true solution to $\displaystyle(\ref{ex1})$ can be found explicitly and is given by

[TABLE]

This true solution may be obtained by using $\displaystyle(\ref{spec_frac}).$ We will first demonstrate how our method can approximate problems such as $\displaystyle(\ref{ex1})$ by taking the trace of $\displaystyle(\ref{approx}),$ and then consider the convergence properties of the approximation of the solution to the extension problem given by $\displaystyle(\ref{b1}).$

As has been demonstrated in this work, it is the case that $\displaystyle(\ref{ex1})$ can be recovered from the following extension problem:

[TABLE]

where $\displaystyle\alpha=1-2s$ and $\displaystyle A\mathrel{\mathop{:}}=\Delta\,:\,H_{0}^{1}(\Omega)\cap H^{2}(\Omega)\subset X\to L^{2}(\Omega).$ In order to implement the problem, we replace the operator $\displaystyle A$ with some finite dimensional discretization. For simplicity, we employ the standard second-order finite difference method on uniform grids [19]. It is worth noting that the use of any other consistent spatial disretization method is acceptable, and as seen below, the convergence rate in the spatial direction will not be affected. The freedom to choose the discretization method for the operator $\displaystyle A$ is one of the primary benefits of this approach.

The tables in Fig. 5.1 demonstrate the numerical convergence rates of the scheme given by $\displaystyle(\ref{approx}).$ For these experiments, we chose a modest value of $\displaystyle M=100,$ in order to demonstrate the efficacy of the method with respect to the truncation. It is clear from this table that the numerical solution converges at the expected rates given by Theorem 4.1. For the sake of clarity, we present a plot of the numerical solution to $\displaystyle(\ref{ex2})$ in Fig. 5.2.

Finally, we finish this example by presenting convergence rates for several other choices of the parameter $\displaystyle s.$ These results are presented in Fig. 5.3 and further verify our theoretical results. For brevity, we omit the convergence rates for the parameter $\displaystyle h,$ as they are not affected by these choices and match those of Fig. 5.1.

Remark 5.1*.*

We wish to note that our experiments have shown that more modest choices of $\displaystyle M$ may still provide accurate numerical solutions, as well. However, we leave a more detailed study of these choices for future endeavors.

Remark 5.2*.*

We also wish to note that one could potentially increase efficiency of the proposed numerical method by allowing it to be run in parallel. For instance, one may rewrite the method as

[TABLE]

where

[TABLE]

While the implementation of $\displaystyle(\ref{para1})$ requires more direct evaluations than the original algorithm, one can easily compute $\displaystyle(\ref{para1})$ in parallel by noting that

[TABLE]

While there is a slight decrease in accuracy due to employing $\displaystyle(\ref{para2}),$ it is the case that for complicated operators $\displaystyle A$ this could potentially be a favorable trade-off.

5.2 A Two-Dimensional Example

Let $\displaystyle\Omega=\{(x,y)\in{\mathbb{R}}^{2}\ :\ x^{2}+y^{2}<1\}.$ We then once again consider the problem $\displaystyle(\ref{ex1}),$ but with

[TABLE]

where $\displaystyle(\lambda_{1,1},\varphi_{1,1})$ are the first eigenpair of the operator $\displaystyle\Delta$ on $\displaystyle\Omega,$ and are given by

[TABLE]

and

[TABLE]

where $\displaystyle J_{1}$ is the first Bessel function of the first kind, $\displaystyle j_{1,1}$ is the first zero of $\displaystyle J_{1},$ and $\displaystyle A_{1,1}$ and $\displaystyle B_{1,1}$ are real-valued constants that ensure $\displaystyle\|\varphi_{1,1}\|_{X}=1.$

By using $\displaystyle(\ref{spec_frac})$ one may determine the true solution of the extension problem to be given by

[TABLE]

where $\displaystyle K_{s}$ is the modified Bessel function of the second kind (see [1] for more information regarding the aforementioned Bessel functions and their properties). Moreover, it can be shown that $\displaystyle j_{1,1}\approx 3.8317.$

Just as before, we discretize the operator $\displaystyle A$ using standard finite differences (see, for example, [20] for more details on doing so for $\displaystyle\Omega$ as above). However, in this example we will only consider convergence aspects of the method in the extended variable direction. A table with convergence rates for various choices of the parameter $\displaystyle s$ are given in Fig. 5.4. As before, the results are in agreement with the rates given in Theorem 4.1.

5.3 A General Second-Order Elliptic Operator

As should be clear from the outlined theory, the proposed numerical algorithm will still be applicable when the operator $\displaystyle A$ is a general strongly elliptic operator on some domain $\displaystyle\Omega\subset{\mathbb{R}}^{d}.$ Following the construction developed in [9], let $\displaystyle\mathscr{L}$ be given by

[TABLE]

where $\displaystyle x\in\Omega$ (here we emphasize the difference between the variable $\displaystyle x$ and the extended variable, $\displaystyle t.$). We assume that the function $\displaystyle c\in L^{\infty}(\Omega)$ with $\displaystyle c\geq 0$ almost everywhere, $\displaystyle A\in C^{0,1}(\Omega,\mbox{GL}(d,{\mathbb{R}}))$ is symmetric and positive definite, and that $\displaystyle\Omega$ is sufficiently smooth. Then, it follows that given $\displaystyle f\in X\mathrel{\mathop{:}}=L^{2}(\Omega),$ there exists a unique $\displaystyle u\in D(\mathscr{L})\mathrel{\mathop{:}}=H_{0}^{1}(\Omega)\cap H^{2}(\Omega)$ such that

[TABLE]

Moreover, we have that $\displaystyle\mathscr{L}$ satisfies Assumption 2.1. It is clear that the operator’s domain is densely defined in $\displaystyle X.$ If we employ the same framework as that in Example 2.1, then we have

[TABLE]

where $\displaystyle[\cdot,\cdot]$ is the standard inner product on $\displaystyle X,$ $\displaystyle C$ is a positive constant, and the final inequality follows from the fact that $\displaystyle A$ is positive definite and $\displaystyle c$ is positive almost everywhere. Thus, we have $\displaystyle\mbox{Re}[-\mathscr{L}u,u]<0$ for all $\displaystyle u\in D(\mathscr{L}),$ as desired. Further, the range condition of a strictly $\displaystyle m-$dissipative operator is satisfied, as well. Thus, it follows that our approximation, given by $\displaystyle(\ref{approx}),$ of the Bessel-type problem

[TABLE]

is well-defined, stable, and convergent. Moreover, it follows that the initial datum, $\displaystyle u,$ will be the solution of the following nonlocal elliptic problem

[TABLE]

It should be noted that we have omitted some more technical details, such as the fact that the operator $\displaystyle\mathscr{L}^{s}$ may be extended to the fractional space $\displaystyle\mathbb{H}^{s}(\Omega),$ since such things are not relevant to the current discussion (this is because we are considering functions with higher regularity to guarantee the theoretical rates of convergence). However, interested readers may see [14], and the reference therein, for more details.

In order to provide an explicit example, we let $\displaystyle d=1,$ $\displaystyle\Omega\mathrel{\mathop{:}}=(-1,1)\subset{\mathbb{R}},$ $\displaystyle A\equiv 1,$ $\displaystyle c\equiv 1,$ and $\displaystyle f(x)=\sin(\pi x).$ Since we do not have an explicit solution to $\displaystyle(\ref{elip_frac})$ on hand, we will approximate the numerical rate of convergence. Let $\displaystyle v_{h}$ denote the numerical solution obtained by $\displaystyle(\ref{approx})$ on a grid defined by the parameter $\displaystyle h.$ Then we define $\displaystyle v_{h/j}$ to be the numerical solution obtained by $\displaystyle(\ref{approx})$ on a grid with spacing $\displaystyle h/j,\ j=2,4.$ We then have the following Milne device for determining the numerical rate of convergence, denoted $\displaystyle p_{h},$ given by

[TABLE]

where the differences are defined on the grid associated to $\displaystyle h$ (the least refined grid) and the norm $\displaystyle\|\cdot\|$ is an appropriately chosen norm (in our case, the discrete $\displaystyle X-$norm).

Remark 5.3*.*

The convergence of the numerical method given by $\displaystyle(\ref{approx})$ is also influenced by the error of truncating domain of the Bessel-type problem. However, by assuming this error is roughly constant, which is reasonable for a fixed value of $\displaystyle M$ and $\displaystyle s,$ the approximation given by $\displaystyle(\ref{num_rate})$ will be accurate.

The table in Fig. 5.5 demonstrates the obtained numerical convergence rates for various choices of $\displaystyle s,$ with $\displaystyle M=100.$ We see from the table that the computed values still reflect the expected theoretical values obtained in Theorem 4.1. Moreover, we present the convergence rate for the case when $\displaystyle s=1,$ to verify the expected second-order convergence of our method.

6 Concluding Remarks

In this article, we consider an approximation to an abstract extension problem, that is the generalization of the method demonstrated by Caffarelli and Silvestre in [7]. This generalized setting allows for the consideration of a wider class of problems, while also allowing for the use of robust nonlinear functional analytic techniques. By taking this novel approach, we develop methods which differ from those presented in [9, 12, 18, 25], yet are quite effective. In particular, this approach allows for the consideration of semigroup theory within the numerical framework. As such, the method is able to exhibit favorable qualities.

Herein, we provide an approach where we approximate the true solution to the extended problem, as compared to the standard approach of discretizing the underlying system first. The approximation is shown to be stable, consistent, and convergent for appropriately smooth initial data. Moreoever, the employed Padé approximation allows for the efficient construction of the numerical solution via factorization techniques. However, one could easily replace this approximation with any desired second-order (or higher) approximation to the associated abstract Cauchy problem.

Moving forward, we will consider improved approximations methods within the same abstract framework. In particular, we will consider abstract analysis of higher-order Padé approximations and quadrature techniques, including Runge-Kutta and exponential integration methods, with the goal of constructing methods of arbitrary order. This will also allow for the consideration of the underlying algebraic structures of these methods, as has been done with Butcher and Lie-Butcher series for classical problems [24]. We also intend to extend these results to abstract problems posed on abstract manifold structures, allowing for the consideration of novel problems in differential geometry. Moreover, we hope to better understand how well this extension technique preserves and approximates the Dirichlet-to-Neumann map that is of fundamental importance in the underlying theory. Preliminary work indicates that our proposed numerical method faithfully represents this operator after applying the conormal derivative given by $\displaystyle(\ref{e2})$, but there is still a need to finish developing a rigorous analysis of this application. However, more immediate concerns will be the deeper comparison of our newly proposed methods with popular existing methods for the fractional Laplacian problem.

Acknowledgements

This work was supported by the NSF grant number 1903450. The author would like to thank Akif Ibraguimov, of Texas Tech University, for introducing this problem to them. The author is also thankful for Akif’s time and insightful suggestions that greatly improved the quality of the work herein. Finally, the author is thankful to the reviewers who provided numerous useful suggestions that greatly improved the presentation of the numerical experiments in Section 5.