Numerical Approximation of Fractional Powers of Elliptic Operators
Beiping Duan Raytcho Lazarov Joeseph Pasciak





 
	In this paper, we develop and study algorithms  for approximately solving the linear algebraic systems: 
	𝒜_h^α u_h = f_h, 0< α <1, for u_h, f_h ∈
	V_h with  V_h a finite element approximation space.
	Such problems arise in finite element or finite difference
	approximations of the  
	problem ^α u=f with , for example, coming from   
	a second order elliptic operator  with homogeneous boundary
	conditions. 
	The algorithms are motivated by  the method of Vabishchevich
	<cit.> that relates the algebraic problem to a
	solution of a time-dependent initial value problem on the interval
	[0,1].   
	Here we develop and study two time stepping schemes 
	based on diagonal Padé approximation to (1+x)^-α.
	The first one uses geometrically graded meshes in order to compensate
	for the singular behavior of the solution for t close to 0. 
	The second algorithm uses uniform time stepping  
	but requires smoothness of the data f_h in discrete norms. 
	For both methods,  we estimate the error in terms of the number of time
	steps, with the regularity of f_h playing a major role for the second method. 
	Finally, we present numerical experiments for _h coming from the
	finite element approximations of  second order elliptic boundary value problems in 
	one and two spatial dimensions. 












[2010]Primary 35S15, 65R20, 65N12, 65N50, 65N30.

Two-field Warm Inflation and Its Scalar Perturbations on Large Scales
    Jian-Yang Zhu
    December 30, 2023
=====================================================================










§ INTRODUCTION



 §.§ Motivation and problem formulation









Nonlocal operators arise in a wide variety  of mathematical
models such as modes of long-range interaction in elastic deformations
<cit.>, nonlocal electromagnetic fluid flows <cit.>, 
image processing <cit.> and many more.
A recent discussion about the properties of such models and their applications to chemistry, 
geosciences, and engineering can be found in <cit.>.


The nonlocal operators considered in this 
paper involve fractional powers of operators  associated with
second order elliptic equations  in bounded domains with homogeneous 
Dirichlet boundary conditions. The fractional power of   
is defined through the Dunford-Taylor integral, <cit.>, 
which is equivalent to the definition by the spectrum of . For a detailed
discussion about this setting and other possible ways to define fractional powers of
the Laplacian (and more general elliptic operators) we refer to <cit.>.  
We focus on issue of solving the corresponding algebraic system 
that arises in approximating such operators by the finite element method, 
e.g. <cit.>.


We begin with the definition of the fractional power of a second order elliptic operator 

in a bounded domain Ω⊂^d, d=1,2,… with a Lipschitz continuous boundary.
On  V × V, with V=H^1_0(Ω),  we consider the bilinear form:

A(w,ϕ)= ∫_Ω(a(x) ∇w ·∇ϕ+q(x) w ϕ)  dx, 
and assume that the coefficients are such that  the bilinear form is  coercive  and bounded on V. 
 Further, we define :L^2(Ω) → V, where v=w ∈ V is the unique solution to
A(w,ϕ)=(v,ϕ),    ϕ∈V.

Here (·,·) denotes the L^2(Ω)-inner product.

Following <cit.>, we define an unbounded operator  with
domain of definition D() being the image of  on
L^2(Ω) and set u= ^-1 u for u∈ D().
This is well defined as  is injective on L^2(Ω).
Negative fractional powers can be defined by  Dunford-Taylor integrals,
i.e., for α>0 and v∈ L^2(Ω),

    ^-αv = 1/2 π i∫ _ z^-α
    R_z()v
      dz
 
where R_z()=(-z )^-1 is the resolvent operator 
and  is an appropriate contour in the complex plane (see, e.g., <cit.>).  


Equivalently, fractional powers for the above example
can be defined by eigenvector expansions.
As  is a
compact, symmetric and positive  definite operator on L^2(Ω), its eigenpairs
{ψ_j,μ_j}, for j=1,2,…,∞, with suitably
normalized eigenvectors, provide an orthonormal basis
for L^2(Ω).  We also set λ_j=μ_j^-1.   For α≥ 0
and v∈ L^2(Ω), 

    ^-α v= ∑_j=1^∞μ_j^α (v,ψ_j)  ψ_j.

Positive fractional powers of  are also given by similar series.
For α≥ 0, define

    D(^α):= { v ∈ L^2  :  λ_j^2α
    |(v,ψ_j)|^2<∞}

and

    ^α  v:= λ_j^α (v,ψ_j)  ψ_j     v ∈ D(^α).


We consider the  fractional order elliptic equation:  Find u∈ D(^α) satisfying

^αu = f      f∈L^2(Ω)
and note that its solution is given by

    u= ^α  f:=μ_j^α (f,ψ_j) ψ_j.
 

Our goal is to approximate u by using finite element or finite differences.
The finite element approximation  on V_h ⊂ V is based on the 
discrete solution operator _h:V_h→ V_h defined  by 
_hv_h:=w_h where w_h is the unique function in
V_h satisfying

    A(w_h,v_h)=(v_h,v_h), v_h∈ V_h.

The inverse of _h is denoted by _h and satisfies

    (_h w_h,v_h)=A(w_h,v_h), v_h∈ V_h.


We obtain a “semi-discrete” approximation to u of equation
(<ref>)
by defining  

u_h =_h^απ_h f := _h^- α π_h f,
where 
π_h is the L^2(Ω)-orthogonal projection into V_h.  Note that u_h
can be expanded in the L^2(Ω)-orthogonal  eigenfunctions of _h, i.e.,  if
{ψ_h,j,μ_h,j}, for j=1,…, M, denotes the eigenpairs and
M is the dimension of V_h then

    π_h f =∑_j=1^M (f,ψ_h,j) ψ_h,j

and 
u_h=∑_j=1^M μ_h,j^α(f,ψ_h,j)  ψ_h,j.

In this paper, we shall study a technique for approximating the solution
to (<ref>) which avoids computing the eigenvectors and
eigenvalues of _h.

We note that the technique to be developed can be applied to finite
difference approximations as well.   
In this case, the discrete space is a finite dimensional space
of grid point values and π_h f is replaced by the interpolant of f
at the grid nodes.  
The matrix _h comes from applying finite difference
approximations to the derivatives in the strong form (see,
(<ref>)) of problem (<ref>).  Its scaling is not arbitrary   if we
expect _h^απ_h f to converge to the grid point values of f.
This, in turn, implies the proper scaling for the discrete norms.

It has been shown in  <cit.> that, if the operator  
satisfies elliptic regularity pickup with index 
s ∈  (0,1] (see Assumption <ref>)
then  for appropriate f,

    u - u_h  =^αf- _h^απ_h f  =O(h^2s).
  
See Section <ref> for more details. 
Here · denotes the L^2(Ω)-norm.

Obviously, the discrete operator _h is symmetric and 
positive definite and the corresponding  matrix is full. 
We note that problems involving
finding approximations of  _h^1/2, cf. <cit.>,
evaluating the sign function of _h, cf. <cit.>, 
and other related functions of matrices have a long history in numerical linear algebra. 



 §.§ The idea of the method of Vabishchevich


In this paper, we develop and study a method for approximating the
solution (<ref>).  Our proposed 
method is related to an  idea of P. Vabishchevich,
<cit.>, which exhibits  u_h
as a solution of a special time dependent problem. Now we 
briefly explain the idea of his method.

We start with the
observation that the unique solution (t) to the ordinary differential
equation initial value problem,
_t+ α(λ-δ)/δ+t (λ-δ) = 0,    
(0) = δ^-α v̂

with 0<δ<λ and α≥ 0, is given by
(t) = (δ+ t(λ-δ))^-α  v̂

and hence
(1)= λ^-α v̂.

Note also that for k>0,
(t+k) = [1+k(λ-δ)/(δ+t (λ-δ))]^-α (t).



Now suppose that the spectrum of _h is contained
in the interval [λ_1, λ_M  ] 
with λ_1>0.   
Let 0<δ<λ_1 and set =_h-δ. We
consider the vector  valued ODE:  Find U(t):[0,1]→ V_h satisfying 
U_t + α(δ+t )^-1 U     = 0,

U(0)      =  δ^-α f_h 

where f_h:=π_h f.  Expanding the solution to (<ref>) as 

    U(t)= ∑_j=1^M c_j(t) ψ_j,h,
 
we find that c_j(t) solves (<ref>) with λ=λ_j,h. 
Moreover, it follows from (<ref>) and (<ref>) that  
U(t):=(δ+ t )^-αf_h 
and 
U(t+k) = (+k (δ+t)^-1)^-α U(t).

As proposed by
Vabischchevich <cit.>, it is then natural to consider numerical approximations to
(<ref>) based on a time stepping method.


In <cit.>, Vabishchevich proposed a time stepping scheme based on
the backward Euler method and applied it to approximate fractional
powers of a discrete approximation _h 
of the Laplace
operator with homogeneous Dirichlet boundary conditions. The results of
numerical computations illustrating the accuracy, convergence, and
some theoretical aspects of the method were provided.










In this paper, we  take a different but related approach.   Instead
of approximating the solution of (<ref>), we simply
approximate the function U(t_i) given by (<ref>) 
on an increasing sequence of nodes 0=t_0<t_1<⋯< t_=1.   
We start from the recurrence (<ref>).
The “time stepping” methods that we
shall study  are based on
diagonal Padé approximation to (1+x)^-α, i.e.,
(1+x)^-α≈r_m(x) := P_m(x)/Q_m(x)

with P_m and Q_m being polynomials of degree m and Q_m(0)=1.
The polynomials P_m and Q_m are then uniquely defined by requiring
that the first 2m+1 terms of the Maclauren expansion of 

    (1+x)^-α-r_m(x)

vanish.   The method that we
study and analyze is  given by setting  U_0:=u(0)=δ^-αf_h and applying the recurrence
U_l=r_m(k_l (δ+t_l-1 )^-1) U_l-1,  l=1,2,…,
with k_l=t_l-t_l-1.  Here U_l is our approximation of U(t_l)
so that U_K approximates U(1)=u_h.
We shall see that these methods are unconditionally stable for
m=1,2,…
and α∈ (0,1).


Even though we take a different point of view, we are still solving (<ref>), 
as suggested by P. Vabishchevich, <cit.>.
It is important to note that although Problem (<ref>) appears
harmless, it behaves considerably different than, for example, the classical parabolic
problem:
w_t + _h w = 0.

For example,  if w(0)= ψ_j,h then w(1)= e^-λ_j,hψ_j,h
while if f_h=ψ_j,h, the solution of (<ref>)  is
u(1)=_h^-αψ_j,h  =λ_j,h^-αψ_j,h.  This means that initial
time step errors in the high frequency components for our problem have much
stronger effect on the accuracy of the final solution. This is especially important
for problems whose solutions have minimal regularity. 



 §.§ Our contributions



In this paper, we consider two time stepping schemes, one involving mesh
refinement near t=0 and the other using a fixed time step.
In both cases, we shall be using (<ref>)  to define our
solution but on different meshes in time.

The refinement scheme starts with an initial basic mesh with t_0=0, t_i=2^i-1-L, for
i=1,2,… L+1 with L chosen so that
2^-L<(λ_M)^-1. Subsequent finer meshes are defined by
partitioning each of the above intervals into N equally spaced subdivisions.   The
refinement scheme leads to an error estimate

    _h^-αf_h -U_N(L+1)≤ C N^-2mf_h
 
for the Padé scheme  based on r_m(x).    

The second scheme that we study is the simpler one using a fixed step
size k_N=1/N.  In this case, we obtain the error estimate
_h^-αf_h -U_N  
 ≤C k_N^α+γA_h^γf_h,   for 0≤γ≤2m-α.



It is clear that L in the first method grows like the logarithm of 
λ_M so that more steps are required by the refinement scheme
when the same N is used in both.   However, in all of our numerical
examples, if one adjusts the values of N in both schemes to obtain the
same absolute convergence, the refinement scheme requires less steps
overall. 

The question of when the norm on the right hand side of (<ref>)
can be controlled by natural norms on the data f is open.  Although,
such a result for γ≤ 1 was provided in <cit.>,
the result for larger γ is not known even in the finite element
case. We discuss this in more detail in Section <ref>. In fact, our numerical results 
in  Section <ref>  suggest that the result is not true in general.




§ PADÉ APPROXIMATIONS
   


In this section, we develop diagonal Padé
approximations to  (1+x)^-α for α∈ (0,1) based on the
classical theory of Padé approximations given by Baker <cit.>.


Our approximations are of the form of (<ref>) with m=1,2,…
and we shall write down explicit formula for the polynomials P_m(x)
and Q_m(x). The starting point is the formula  <cit.>
or <cit.>:


    (1+x)^-α =  (α,1;1;-x).

Here (a,b;c;x) denotes the hypergeometric function defined by

    (a,b;c;x)=∑_j=0^∞(a)_j (b)_j/j! (c)_j x^j.

Here (a)_0=1 and (a)_j=a(a+1)⋯(a+j-1) for j>0.
This series converges for |x|<1 provide that c is not in
{0,-1,-2,…}.  

Baker, <cit.>, also gives an explicit expression 

for the denominator:
Q_m(x)   = (-m,-α-m; -2m;-x)

   :=∑_j=0^m (-m)_j  (-α-m)_j/j!  (-2m)_j (-x)^j 
= 1 +∑_j=1^m  a_j b_j(α) x^j.

Here b_0(α)=a_0=1,

    b_j(α)= (m+α)((m-1)+α)⋯((m+1-j)+α)

and

    a_j=m(m-1)⋯(m+1-j)/j! (2m(2m-1)⋯(2m+1-j)    for j=1,2,…,m.


Now, Theorem 9.2 of <cit.> implies that Q_m(x)/P_m(x) is the
diagonal Padé approximation to (1+x)^α and again applying
(5.12) of <cit.>, we find that 
P_m(x)   =_2F_1(-m,α-m; -2m;-x)

   :=∑_j=0^m (-m)_j
	(α-m)_j/j! (-2m)_j (-x)^j 
= 1 +∑_j=1^m  a_j b_j(-α) x^j.

Using the above formulas, we find, for example:

    r_1(x)= 1 + [(1-α)/2]x / 1 + [ (1+α)/2 ]x

and 

    r_2(x)= 1 + [ (2-α)/2]x +[(2-α)(1-α)/12]x^2 / 1 +
    	[(2+α)/2]x+ [(2+α)(1+α)/12]x^2 .


For our further considerations, we need to discuss a relation
between the  denominators appearing in Padé approximations and
orthogonal polynomials with respect to an appropriate weight w.  
We first note the series expansion of 

    (1-x)^-α = ∑_i=0^∞(α)_i  x^i/i!:=
    ∑_i=0^∞ c_i x^i.

The coefficients above satisfy
c_i = 1/Γ(α)Γ(1-α) ∫_0^1 x^i
((1-x)^-α x^α-1 )   dx 
:= 1/Γ(α)Γ(1-α) ∫_0^1 x^i
w(x)    dx, 

where w(x)=(1-x)^-αx^α-1. By (7.7) of
<cit.>, utilizing the fact that P_m(-x)/Q_m(-x) is  the
Padé approximation of (1-x)^-α, we obtain that the denominator Q_m(x) can be expressed by 

    Q_m(x)=(-x)^m q_m(-1/x)

where q_m is the monic polynomial of order m which is orthogonal to the set of polynomials 
of degree less than m with the weight xw(x). As the roots of q_m(x) are all in the interval (0,1), those of
Q_m(x) are in the interval (-∞,-1).

The expressions for the numerator and denominator in r_m(x) imply the following
proposition.

  
	Let α be in (0,1) and m be a positive integer.
	Then, there are positive constants ρ_m 
	satisfying
	ρ_m≤r_m(x)≤1,    for all x≥0.
	
	


	
	As α is in (0,1), 
	
    0<b_j(-α)<b_j(α).

	This means for x≥ 0, 
	a_j b_j(-α)x^j≤ a_j b_j(α) x^j 
	with strict inequality when x>0.  The second inequality of  (<ref>) follows by summation. 
	
	For the first inequality in (<ref>), we note that both Q_m(x) and P_m(x) are positive for
	x≥ 0 and 
	
    lim_x→∞ r_m(x) = b_m(-α)/b_m(α)>0.

	The first inequality in (<ref>) immediately follows from the fact
	that r_m(x) is continuous on [-1,∞) since all of the roots of
	Q_m(x) are in (-∞,-1) and  takes values in (0,1] for x≥ 0. 



To clarify further the convergence of these approximations, we include
the following proposition.

   
	For 0≤ s ≤ 2m+1,
	|(1+x)^-α -r_m(x)| ≤c_m,s   x^s   for x∈[0,∞)
	where 
	
    c_m,s=max{ Q_m(-1) 2^s-2m,2^1+s}.





	We use Theorem 5 of <cit.> which gives that 
	
    (1+x)^-α-r_m(x)    = Q_m(-1)/Q_m(x)∑_n=2m+1^∞(α)_n(n-2m)_m/n! (n+α-m)_m (-x)^n
       :=Q_m(-1)/Q_m(x) (-x)^2m+1∑_i=0^∞ d_2m+1+i (-x)^i.

	A simple computation shows that 0<d_n<1 for n≥ 2m+1 and so for |x|<τ<1,
	
    |∑_i=0^∞ d_2m+1+i (-x)^i  |≤ (1-τ)^-1.

	Thus,  for x∈ [0,τ] and s∈ [0,2m+1],
	
    |r_m(x)-(1+x)^-α |    ≤ |Q_m(-1)| (1-τ)^-1 x^2m+1
       ≤ |Q_m(-1)| τ^2m+1-s (1-τ)^-1 x^s.

	Finally, for the case x ∈[τ, ∞)  by  Proposition <ref>,
	
    |r_m(x)-(1+x)^-α|≤ 2≤ 2x^s/τ^s,   
    		for  x∈ [τ,∞)

	and (<ref>) follows by  taking τ=1/2.





§ THE TIME STEPPING SCHEMES AND THEIR ANALYSIS


In this section, we define and analyze both equally spaced time stepping 
schemes as well as schemes employing refinement near the origin.    
We shall restrict ourselves to the approximating the solution to finite dimensional 
problem (<ref>) described in the introduction even though generalizations to hermitian
and non-hermitian bounded operators on infinite dimensional spaces are possible.  
Recall that =_h-δ with δ∈
(0,λ_1) and that · and (·,·) denote, respectively, 
the norm and inner product in L^2(Ω).



 §.§ Time-stepping method on geometrically refined meshes (GRM)

We  first consider the geometrically refined mesh that is constructed in two steps.

First, we take 

    =⌈log(λ_M)/log(2) ⌉

and set t_i=2^i-1-L for i=1,…,+1 and t_0=0. Note that t_+1=1 and 
2^-≤ (λ_M)^-1.
Next, we define 

    k_n= { t_n/N   :    for n=1,…,,
    
    t_1/N   :    when n=0. .

Note that k_0=k_1 and 
k_0 λ_M = 2^-
	λ_M/N≤N^-1.

The grid that we use in our computations is obtained by partitioning 
each subinterval 
I_n=[t_n,t_n+1],    n=0,1,…,,
into N subintervals with end points

    t_n,j := t_n + jk_n     for  j=0,…,N.

It can be seen that t_n,0=t_n and t_n,N=t_n+1, n=0,1,…,,
so that the mesh has totally (+1)N intervals.

Based on this partitioning, the refined time stepping method for
approximating A_h^-αv for v∈ V_h is given by:



After executing the above algorithm, U_+1 is the approximation to
_h^-α.  Note that the notation differs slightly from
that used in the introduction.  The computation of U_+1 requires
=( +1)N
time steps.

The discrete eigenvalues and eigenvectors will play a major role in our
analysis so, for notational simplicity, we denote them by
{ψ_j,λ_j}  (instead of {ψ_j,h,λ_j,h} as in
the introduction).
We then have 

    =∑_j=1^ (,ψ_j) ψ_j.

Moreover,

    _h^-α v =∑_j=1^λ_j^-α (,ψ_j) ψ_j.

The expansion for U_L+1 is given by

    U_L+1=∑_j=1^μ(λ_j)  (,ψ_j) ψ_j

where the coefficient μ(λ_j) is given by the following algorithm:





	Let N be a positive integer.   Then
	
    |λ^-α - μ(λ)| ≤c N^-2m, λ∈ [λ_1,λ_M]

	and 
	
    _h^-α -U_ +1≤c N^-2m, ∈ V_h.

	Here c is a
	constant depending only on δ, m and α.


  Fix λ in  [λ_1,λ_M] and let
	{μ_n,μ_n,j} be as in Algorithm <ref>. Further,
	for n=0,…, and j=1,…,N let 
	
	
    v_n,j=(δ +t_n,j(λ-δ))^-α,   
    	e_n,j=v_n,j-μ_n,j,    θ_n,j= k_n (λ-δ)
    	/(δ+t_n,j-1(λ-δ)).
 
	We note that as in (<ref>), 
	
    v_n,j=(1+θ_n,j)^-αv_n,j-1.

	Thus,
	
    e_n,j = r_m(θ_n,j) 
    	e_n,j-1 + [(1+θ_n,j)^-α-r_m(θ_n,j)] v_n,j-1

	and  by Proposition <ref> and <ref>,
	
    |e_n,j| ≤ 
    		|e_n,j-1| +
    		|(1+θ_n,j)^-α-r_m(θ_n,j)| |v_n,j-1|,

	for n=0,1,…,L and j=1,…,N.
	
	Using  e_0=0, (<ref>),  Proposition <ref> and the
	definition of v_n,j-1,
	we have 
	
    |λ^-α-μ(λ)|    = |e_+1|=∑_n=0^∑_j=1^N(
    	|e_n,j|-|e_n,j-1|)
       ≤ c_m,2m+1∑_n=0^∑_j=1^Nθ_n,j^2m+1
    	|v_n,j-1|
       = c_m,2m+1∑_n=0^∑_j=1^N(k_n
    		(λ-δ))^2m+1/(δ+t_n,j-1(λ-δ))^2m+α+1.

	
	We note that by (<ref>),
	∑_j=1^N (k_0
		(λ-δ))^2m+1/(δ+t_0,j-1(λ-δ))^2m+α+1
	    ≤δ^-2m-α-1∑_j=1^N N^-2m-1 =δ^-2m-α-1 N^-2m.
	
	The remaining terms in (<ref>) will be
	bounded by integration. Applying mean value theorem for integration it
	follows that for some t_θ∈ [t_n,j-1,t_n,j],
	
    ∫_t_n,j-1^t_n,jt^2m (λ-δ)^2m/(δ+t (λ-δ))^2m+α+1  dt= t_θ^2m (λ-δ)^2m/(δ+t_θ (λ-δ))^2m+α+1k_n
       =(t_θ/t_n,j-1)^2m(δ+t_n,j-1 (λ-δ)/δ+t_θ (λ-δ))^2m+α+1t_n,j-1^2m (λ-δ)^2m/(δ+t_n,j-1 (λ-δ))^2m+α+1k_n
       ≥(t_n,j-1/t_θ)^α+1t_n,j-1^2m (λ-δ)^2m/(δ+t_n,j-1 (λ-δ))^2m+α+1k_n
       ≥(N/N+1)^α+1t_n,j-1^2m (λ-δ)^2m/(δ+t_n,j-1 (λ-δ))^2m+α+1k_n.

	Now  (N+1)/N ≤ 2 and t_n≤ t_n,j-1 for j=1,2,…,N so
	that
	(k_n
		(λ-δ))^2m+1 /(δ+t_n,j-1(λ-δ))^2m+α+1
	   =N^-2m t_n^2m k_n  (λ-δ)^2m+1/(δ+t_n,j-1 (λ-δ))^2m+α+1

	   -80pt≤2^α+1 N^-2m   (λ-δ)∫_t_n,j-1^t_n,j  t^2m (λ-δ)^2m/(δ+t (λ-δ))^2m+α+1  dt.
	
	Applying (<ref>) to the remaining terms in (<ref>) we have
	
    ∑_n=0^∑_j=1^N(
    	|e_n,j|-|e_n,j-1|) 
       ≤δ^-2m-α-1 N^-2m+ 2^α+1N^-2m 
    	(λ-δ)∫_0^1  t^2m (λ-δ)^2m/(δ+t (λ-δ))^2m+α+1  dt
       ≤δ^-2m-α-1 N^-2m+2^α+1 N^-2m∫_0^∞z^2m/(δ+z)^2m+α+1  dz.

	As the last integral on the right converges,  (<ref>) follows
	combining  (<ref>), (<ref>) and (<ref>).
	
	The inequality (<ref>) follows immediately from (<ref>) and  the Parseval's identities:
	
	
	 _h^-α -U_+1^2   =∑_j=1^|λ_j^-α-μ(λ_j)|^2 |( ,ψ_j)|^2 

	   ≤c^2 N^-4m
	∑_j=1^|(v,ψ_j)|^2=c^2 N^-4m   ^2.
	
	


   For any s∈ and v∈ V_h, let
	
    _s,h:=(∑_j=1^λ_j^s
    		|( ,ψ_j)|^2)^1/2=_h^s/2.

	Then,  it follows from the proof of the above theorem that for any
	s∈, 
	
    _h^-α   -  U_+1_s,h≤c N^-2m_s,h.








 §.§ Time-stepping method on uniform meshes (UM)

We next consider uniform time stepping.  In this case, given a positive
integer N, we set k_N=1/N and t_n=nk_N. The approximation is
obtained from the recurrence
U_0   =δ^-α ∈V_h,

U_n   = r_m(k_N (δ+t_n-1 )^-1) U_n-1,   for 
n=1,2,…,N.

In this case, U_N is our approximation to _h^-α.  The
analysis of the error requires the following proposition.

  For λ≥δ,
	set
	θ_n= k_N (λ-δ)/δ+t_n-1(λ-δ).
	
	Then for q≥ p>1,  
	
    ∏_n=p^q r_m(θ_n)≤ c ∏_n=p^q (1+θ_n)^-α,
 
	with c depending only on m and α.



	We note that for n>1, θ_n≤ (n-1)^-1∈ [0,1].
	By Proposition  <ref>,  
	
		
		r_m(θ_n)   ≤(1+θ_n)^-α+c_m,2m+1|θ_n|^2m+1

		   ≤(1+θ_n)^-α (1+2^αc_m,2m+1 (n-1)^-2m-1)
		
        
	and hence
	
		 r_m(θ_n)/(1+θ_n)^-α≤(1+2^αc_m,2m+1
		(n-1)^-2m-1) .
	
	Thus,  for q≥ p>1, 
	
	∏_n=p^q r_m(θ_n)/(1+θ_n)^-α≤∏_n=p^q  (1+2^αc_m,2m+1(n-1)^-2m-1)  ≤c
	
	with
	
		c=  ∏_j=1^∞  (1+2^αc_m,2m+1
		j^-2m-1).
	



Theorem 7.2 of <cit.> provides error estimates for single
step approximations for the standard parabolic problem (<ref>) with non-smooth
initial data.     
The next theorem
has the same flavor however differs
significantly as the solutions of our problem exhibit  less regularity.


	Let N>1, ∈ V_h and U_N be defined by (<ref>).  Then,
	for γ≥0 and  α+γ≤ 2m,
	 _h^-α -  U_N  ≤c k_N^α+γ _h^γ
	
	with c depending only on α, m,  γ, and δ.
	As in Remark <ref>, the left hand norm above can be replaced by
	·_h,r provided that the right is replaced by
	_h^γ_h,r.




   In this proof, c denotes a generic positive constant
	only depending on α,
	γ, m and δ.    We fix λ≥δ and define, 
	for j≥ l≥ 1, 
	
    r_m^j,l(λ):=r_m^j,l=r_m(θ_j) r_m(θ_j-1) ⋯ r_m(θ_l)

	and 
	
    w^j,l(λ):=w^j,l=(1+θ_j)^-α (1+θ_j-1)^-α ⋯ (1+θ_l)^-α.

	Finally, we set
	
    e_N(λ):=e_N=(w^N,1-r_m^N,1) δ^-α.

	
	We note that it is a consequence of (<ref>) and  (<ref>) that 
	δ^-α  w^N,1 =
	λ^-α. 
	
	For any j>1,
	
    w^j,1-r_m^j,1=  [(1+θ_j)^-α-r_m(θ_j)] w^j-1,1
    	+r_m(θ_j) [w^j-1,1-r_m^j-1,1].

	Repeated application of this identity leads to 
	e_N      = δ^-α [w^N,1-r_m^N,1] 

	    =   δ^-α[(1+θ_N)^-α - r_m(θ_N)] w^N-1,1 

	       +
	δ^-α r_m(θ_N) [(1+θ_N-1)^-α - r_m(θ_N-1)] w^N-2,1

	       + ⋯

	       + δ^-α  r_m^N,2 [(1+θ_1)^-α -  r_m(θ_1)].  
	
	Proposition <ref>  implies that for j≥ 2,
	|r_m^N,j|≤c w^N,j.
	We first bound the last term of  (<ref>) by applying this
	and (<ref>) to obtain
	T_1:=δ^-α|r_m^N,2 [(1+θ_1)^-α -
	r_m(θ_1)]|     ≤c δ^-α w^N,2| (1+θ_1)^-α -
	r_m(θ_1)|

	   = cλ^-α (1+θ_1)^α|(1+θ_1)^-α -
	r_m(θ_1)|.
	
	When θ_1≤ 1,  Proposition  <ref> with s=α+γ gives,
	
    T_1 ≤ c θ_1^α+γλ^-α.

	When θ_1>1, since γ≥ 0,
	
    T_1 ≤ c (1+θ_1)^αλ^-α≤ c θ_1^α+γλ^-α.

	Thus, in either case, since θ_1 = k_N (λ-δ)/δ≤
	k_N λ/δ,
	T_1≤c k_N^α+γ λ^γ.
	
	The absolute value of the other terms in (<ref>) are given by
	
    T_j:=| r_m^N,j+1  [(1+θ_j)^-α - r_m(θ_j)]
    	w^j-1,1|,     for  j=2,3,… N,

	where we have defined r_m^N,N+1=1 for convenience of notation.
	Similar to (<ref>), we have 
	
    T_j ≤ c λ^-α (1+θ_j)^α |(1+θ_j)^-α -
    	r_m(θ_j)|.

	In this case, θ_j is in [0,1]
	and we apply Proposition <ref> with
	s=1+α+γ to obtain
	
    T_j≤ c λ^-αθ_j^1+α+γ≤
    	c k_N (λ-δ) k_N^α+γλ^γ/(δ+t_j-1 (λ-δ))^1+α+γ.

	Thus,  
	∑_j=2^N  T_j ≤c (λ-δ) k_N^α+γ
	λ^γ
	∫_0^1(δ+t
	(λ-δ))^-1-α-γ  dt
	≤c  k_N^α+γ
	λ^γ.
	
	Combining (<ref>) and (<ref>) gives
	
    |e_N(λ)| ≤ c  k_N^α+γλ^γ.

	
	We note that by (<ref>),
	
    _h^-α = δ^-α∑_j=1^ w^N,1(λ_j)
    	( ,ψ_j)ψ_j
    	       
    	U_N= δ^-α∑_j=1^ r_m^N,1(λ_j)
    	( ,ψ_j)ψ_j .

	Thus,
	
    _h^-α -U_N^2    = ∑_j=1^ |e_N(λ_j)|^2
    	|( ,ψ_j)|^2
       ≤  c k_N^2α+2γ∑_j=1^λ_j^2γ |(v,ψ_j)|^2 = c k_N^2α+2γ_h^γ^2

	and this completes the proof.


As seen in the above theorem, the discrete regularity of the 
solution determines the rate of convergence for the uniform step size  time stepping
method.  To some extent, the regularity of the discrete solution is related to the
regularity properties of the continuous problem which is being
approximated. This will be discussed in the next section.  




§ FINITE ELEMENT APPROXIMATION TO FRACTIONAL POWERS OF SECOND ORDER ELLIPTIC OPERATORS



We start with the second order elliptic problem  associated with the
bilinear form (<ref>) of the introduction, namely 
the  boundary value problem:  

-∇·(a(x)∇w) + q(x) w    =f,    for x∈Ω,

w(x)   =0,   for x∈∂Ω.
Here q(x),a(x) and Ω are as in the introduction.  The bilinear
form (<ref>) results from (<ref>) in the usual way, i.e.,
integration against a test function and integration by parts.

We
start by providing some results for the error between the
semi-discrete approximation u_h given by (<ref>) and the
solution u of (<ref>).    These results depend on the following
regularity condition:


  
	 satisfies elliptic regularity pickup with index s ∈ (0,1],
	that is 
	
		
  * For f∈ H^-1+s(Ω), f is in H^1+s(Ω) and
		there is a constant c not depending on f satisfying
		
    f _H^1+s(Ω)≤ c f_H^-1+s(Ω).

		
  * =^-1 is a bounded map of H^1+s(Ω) into
		H^-1+s(Ω).
	


  
	In <cit.>, it has been shown that this implies
	
    D(^t/2) =H^1_0(Ω)∩ H^t(Ω),    for   t ∈[1,1+s]

	with equivalent norms.


The above remark shows that D(^t/2) coincides with a Sobolev space
of index t.   Accordingly, we introduce the notation 

    ^t=D(^t/2),    for  t≥ 0.


A detailed estimation of the error u - u_h= ^αf- _h^απ_h f
can be found  in <cit.> and is  summarized below
(see also, <cit.> for the case when s=1).

 
	(see, <cit.>) Let Assumption <ref>
	hold and assume that the mesh is globally quasi-uniform so that the 
	inverse inequality holds, e.g. <cit.>
	or <cit.>.
	Set β =s -α when s > α and β=0 when s ≤α.
	For γ≥β, there is a constant C uniform in h such that
	
	^αf- _h^α π_h f    ≤C_h,γh^2s  f_^2γ
	f ∈^2 γ,
	where
	
    C_h,γ = {
    	C log(1/h),        when  γ = β   s ≥α,
    
    	C,        when  γ > β   s ≥α ,
    
    	C,       when  α > s. .



As we see from this theorem, the rate of convergence in the
L^2(Ω)-norm is 
the result of 
an interplay between the fractional order α, the regularity 
pick up s of the solution of problem  (<ref>), and the 
regularity of the right hand side f.  The bottom line is that one recovers optimal 
convergence rate O(h^2s) for α > s when f ∈ L^2(Ω).  
However, if α < s, the solution is not in H^2s(Ω)
without  extra regularity from f so this additional smoothness is
needed  to get the same rate.

Now if we approximate the problem (<ref>) using the method of Vabishchevich 
on the geometrically refined mesh as described in  Algorithm <ref>, we get the
following bound for the total error (approximation by finite elements and time-stepping):
  
	Assume the conditions of Theorem <ref> hold.
	Then  U_L+1, obtained by  Algorithm <ref> with v=π_h f,
	=_h-δ and
	performing K=(L+1)N steps, satisfies 
	
    ^α f - U_L+1≤ C (h^2sf _^2γ + N^-2mf ).




	Assume the conditions of Theorem <ref> hold. Then
	U_N, obtained by applying (<ref>) with v=π_h f and
	=_h - δ, 
	satisfies 
	
    ^α f - U_N≤ C (h^2sf _^2γ
    	+ N^-α-β_h^βπ_h f),   for
    	0≤β≤ 2m-α.

	


We next consider the question of bounding the norm _h^βπ_h f in 
terms of the regularity of f.   For β∈ [0,1/2], this reduces to showing 
that the L^2(Ω)-projector into V_h is a bounded operator on H^1_0(Ω) 
with bound independent of h. For globally quasi-uniform meshes, this result is given in
<cit.> while the case of certain refined meshes is given in 
<cit.>. When the H^1 bound holds, by interpolation, there is a constant c
depending only on β∈ [0,1/2] satisfying 
_h^βπ_h v  ≤c ^βv, v∈D(^β).



We extend the above inequality to β∈ [1/2,(1+s)/2] in the next lemma whose
proof is included for completeness as it was already observed in
<cit.>.

 Assume that Assumption <ref> holds and that the
	mesh is globally quasi-uniform.   Then there is a constant c  depending on
	β∈ [1/2,(1+s)/2] such that (<ref>) holds.


 
	Let P_h:H^1_0(Ω)→ V_h denote the elliptic
	projector, i.e., P_hw=w_h∈ V_h is the unique solution of 
	
    A(w_h,θ_h)=A(w,θ_h),θ_h∈ V_h.

	Without loss of generality, we can take the norm on H^1_0(Ω) to
	be 
	
    v_H^1(Ω) = A(v,v)^1/2=^1/2 v.

	We then have  
	
    _h^1/2P_h w=  P_h w_H^1(Ω)≤w_H^1(Ω)=^1/2w,

	for all w∈ H^1_0(Ω)=D(^1/2), while the identity 
	_h P_h v=π_h  v for v∈ D() implies that
	
    _h P_h v = π_h  v ≤ v, v∈ D().

	It follows by interpolation that for r∈ [1/2,1],
	
    ^r_hP_h v ≤^rv,  v∈ D(^r).

	Now for t∈ [1,1+s], Remark <ref> implies that for ^t=D(^t/2)=H^t(Ω)∩ H^1_0(Ω).  
	Thus, for v∈^t ,  
	
    π_h v _h,t≤(π_h-P_h)v_h,t +P_h v_h,t≤ 
    	C (h^1-t(π_h-P_h)v_H^1(Ω) +v_^t(Ω))

	where we also used the inverse inequality for the last inequality
	above.   The inequality (<ref>) for β = t/2∈
	[1/2,(1+s)/2]  follows from  the above inequality,
	the triangle inequality and the well know error estimates
	
		(I-π_h)v_H^1(Ω)+(I-P_h)v_H^1(Ω) ≤ch^t-1
		v_H^t(Ω).
	






§ NUMERICAL EXAMPLES
 


In this section, we present numerical examples for the problem 
coming from (<ref>) with a(x)=1, q(x)=0 and
Ω⊂^d, for d=1,2.




 §.§ One dimensional examples:  Ω=(0,1)

We will consider approximating the solution of  (<ref>)  with the following
choices of f:

	
  * f= exp(-1/x-1/(1-x)+4)  so that f∈^s
	for any s.
	
  * f= x(1-x)  so that f∈^s for s<5/2.
	
  * f=min(x,1-x)  so that f∈^s for s<3/2.
	
  * f=1, so that f∈^s for s<1/2.

We note that f=1 fails to be in ^1/2 since functions in
^1/2 vanish at x=0 and x=1.

The first set of runs demonstrate the time stepping error behavior using the
geometric refined time stepping algorithm (GRM) and the uniform time
stepping scheme (UM) for various α and N.  In this case, we use
a fixed equally spaced mesh with h=1/1000 and 

L=⌈ 2  |log h |/log 2⌉.    
The total number of time
steps for the GRM scheme is thus (L+1)N, with N=1,2,4,8,….  To
make the comparison more meaningful, we report the errors obtained using
the UM and GRM algorithms  as a function of the number of solves.


For the first plot, we use f given by (c) above and α=0.1,0.5,0.9.  
Figure <ref> gives  log-log plots of
the relative L^2 error between u_h=𝒜_h^-απ_h f and the
result obtained using the refined and uniform time stepping schemes
with m=1 (left plot) and m=2 (right plot) as a function of the
number of solves.  Note that the GRM method leads to smaller error
using the same number of solves.   
Plots for f given by (b) and (d) are similar and
are omitted.  




To further demonstrate that the numerical results reflect the
theoretical results proved earlier, we report the approximate order of
convergence going from N=n to N=2n, 
Approximate order:= log(E_n/E_2n)/log(2).
We used n=8 for both the UM runs and the GRM runs. The reason that we chose this n  for this computation is that the errors were getting so small in the GRM algorithm for N≥ 32 that, we suspect, computer round off  was effecting their significance.






We report the approximate order of convergence computed using (<ref>) and  compare it with the theoretical rate (in parenthesis) in 
Table <ref> and Table <ref>. Note that the approximate order of convergence was under the assumption that  inequality (<ref>) holds.  Table <ref> and  Table <ref> give the rates when m=1 and m=2, respectively for varying α and f
given above. In most cases, the computed order is in good agreement with the theoretical rate for both the refinement and uniform time stepping schemes. In contrast, the theoretical rate of the smooth problem (for f given by (a)) would be 2m if (<ref>) held. The results in Table  <ref> suggests that (<ref>) does not hold uniformly for β=4-α. In all of the above examples, the error observed for the refinement scheme as a function of the number of solves was below that of the uniform time stepping scheme.  






Note that the convergence of the GRM schemes is more robust than that of
the UM schemes. The GRM schemes always yield 2m'th order convergence
while the convergence rate of UM schemes are related to the parameter 
α and the the (discrete) regularity of the initial data v as
suggested by the theory.  The advantages of the refinement scheme are
especially evident for the non-smooth initial data problem.



 §.§ A spatial refinement example


The last one dimensional example is for f=1 but uses a sequence of 
refined spatial grids. By  (<ref>), the semi-discrete error 
for an unrefined mesh is O(h^2s) for s<1/4+α.  As the
singular behavior is at 
the endpoints of the interval, it is natural to use refinement there
to try to improve the error behavior.  We consider a mesh
resulting from a geometric refinement near 0 and 1 
similar to the geometric time stepping
refinement at 0. Specifically,  our meshs on [0,1/2] are constructed
by restricting the mesh of Subsection 3.1 to [0,1/2]  as a function of
N,  the number of points per interval. In this
construction,  we choose   L 
so that 2^-L≤
h^-2 where h=1/4N is the mesh size on [1/4,1/2].   The mesh on
[1/2,1] is obtained by reflecting the mesh on [0,1/2] about 1/2. The
number of mesh points in space is O(h^-1log (1/h)).

Table <ref> 
reports errors using the GRM and UM
time stepping schemes applied to the  case when _h comes from a
sequence of refined spatial meshes as discussed above.   For brevity, we only report
results for α=0.5. 
For each spatial mesh, we compute an accurate approximation
u_h,ref to the
semi-discrete solution u_h= _h^-απ_h f
by using a highly refined (in time)
4'th order GRM time-stepping scheme.   We then report the
semi-discrete error norm 
e_semi:=I_h(u) -  u_h,ref where I_h denotes
the finite element interpolation operator on the refined spatial mesh.   The solution u is
computed at the nodes
by using 800000 terms in its Fourier series expansion.
The error e_semi is important as it gives us an idea how small we
need to make the time stepping error so that the overall error is, for
example, less than or equal to 2e_semi.
In Table <ref>,  
nx is  the number of intervals in the
spatially refined  grid,    NS is the number of time steps
used to reduce  the GRM error
E_GRM:=u_GRM-u_h,ref
below e_semi
and E_UM:=u_N-u_h,ref is the UM error for N=10^5 time steps.
It is  clear that the uniform time stepping method is
 inefficient for this problem. Indeed,  in many cases,  
the uniform time stepping fails to reduce the error 
below  e_semi 
even when using 10^5  time steps.






 §.§ Two dimensional examples: Ω=(0,1)^2

In the two dimensional case, we consider f given by: 

	
	
  * f(x,y)= x(1-x)y(1-y)   so that   f∈^s for s<5/2;
	
  * f(x,y)={
	1      0.25≤ x,y≤ 0.75,

	0     . 
	    so that f∈^s for s< 1/2.


For brevity, we only report results for the case of m=2.  For all runs, we use a uniform mesh in space of size 
h=1/100 and L=14 in the GRM case. We report  the relative error of the time stepping solution compared with
u_h=𝒜_h^-απ_h f. Table <ref> gives the  errors as a function of α 
and NS, the number of time steps, for the uniform stepping approximation while those of  Table <ref> are for the 
geometric stepping approximation.  Similar to 1-D case, the reported convergence rates are obtained by (<ref>)
using n:=NS=8× 15 with the theoretical rates in parenthesis.










§ CONCLUSIONS


 We proposed two time-stepping methods based on Padé approximation for solving a special pseudo-parabolic 
 equation  introduced by Vabishchevich for solving equations involving powers of symmetric positive elliptic operators. 
 We consider two schemes that use geometrically refined  and uniform meshes in time.  The scheme that uses geometrically refined mesh has a   
 convergence rate that does not depend on the smoothness of the solution, while the scheme involving
 uniform time-mesh depends crucially on the discrete regularity of the solution. Both, the theoretical estimates and the numerical tests show that the scheme  
 on geometrically refined meshes  is more efficient compared with the uniform time-stepping scheme, especially in the non-smooth data case.



§ ACKNOWLEDGMENTS

B. Duan is supported by China Scholarship Council and the Fundamental Research Funds for the Central Universities of 
Central South University (2016zzts015).  The work of R. Lazarov was supported in part by  NSF-DMS #1620318 grant.
 













abbrv

