Processing math: 85%

Work Packages Models involving memory terms and hybrid PDE+ODE systems

Finite Element approximation of the one-dimensional fractional Laplacian

Authors: - 07 June 2019

We describe here a Finite Element algorithm for the approximation of the one-dimensional fractional Laplacian (d2x)s on the interval (L,L), L>0 and for the numerical resolution of the following fractional Poisson equation

{(d2x)su=f,x(L,L)u=0,x(LL)c=:R(L,L.)

This algorithm has also been emploied in [2,3] for the numerical controllability of fractional parabolic problems (see our previous entries in the DyCon Blog, WP3_P0001 and WP3_P0022).

We recall that the fractional Laplacian is defined, for all s(0,1) and any function u regular enough, as the following singular integral

(d2x)su(x):=csP.V.Ru(x)u(y)|xy|1+2sdy,

with cs an explicit normalization constant given by

cs=s22sΓ(1+2s2)πΓ(1s),

Γ being the usual Euler Gamma function.

The starting point of the FE method is the definition of a bilinear form associated to the fractional Laplace operator (d2x)s and the introduction of the variational formulation corresponding to the elliptic problem (1).

This variational formulation is given as follows: find uHs0(L,L) such that the identity

a(u,v)=LLfvdx

is satisfied for any function vHs0(L,L). In (2), with Hs0(L,L) we indicate the space

Hs0(L,L):={uHs(R):u=0in(L,L)c},

Hs(R) being the usual fractional Sobolev space. Moreover, the bilinear form

a:Hs0(L,L)×Hs0(L,L)R

is defined as

a(u,v)=cs2RR(u(x)u(y))(v(x)v(y))|xy|1+2sdxdy.

Let us describe how, starting from the above bilinear form, we can obtain our FE approximation. For simplicity, in what follows we will take L=1, although the methodology remains the same for any real L>0.

Let us take a uniform partition of the interval (1,1) as follows:

1=x0<x1<<xi<xi+1<<xN+1=1,

with xi+1=xi+h, i=0,N. We call M the mesh composed by the points xi:i=1,,N, while the set of the boundary points is denoted M:=x0,xN+1.

Now, define Ki:=[xi,xi+1] and consider the discrete space

Vh:={vhHs0(1,1)|vh|KiP1},

where P1 is the space of the continuous and piece-wise linear functions. Hence, we approximate (2) with the following discrete problem: find uhVh such that

cs2RR(uh(x)uh(y))(vh(x)vh(y))|xy|1+2sdxdy=11fvhdx,

for all vhVh. If now we indicate with

{ϕi}Ni=1

a basis of Vh, it will be sufficient that the above equality is satisfied for all the functions of the basis, since any element of Vh is a linear combination of them. Therefore the problem takes the following form

cs2RR(uh(x)uh(y))(ϕi(x)ϕi(y))|xy|1+2sdxdy=11fvhdx,i=1,,N.

Clearly, since uhVh, we have uh(x)=Nj=1ujϕj(x), where the coefficients uj are, a priori, unknown. In this way, (4) is reduced to solve the linear system Ahu=F, where the stiffness matrix AhRN×N has components

ai,j=cs2RR(ϕi(x)ϕi(y))(ϕj(x)ϕj(y))|xy|1+2sdxdy,

while the vector FRN is given by F=(F1,,FN) with

Fi=f,ϕiL2(1,1)=11fϕidx,i=1,,N.

Moreover, the basis

{ϕi}Ni=1

that we will employ is the classical one in which each ϕi is the tent function with supp(ϕi)=(xi1,xi+1) and verifying ϕi(xj)=δi,j. In particular, for xxi1,xi,xi+1 the ith function of the basis is explicitly defined as (see Figure 1)

ϕi(x)=1|xxi|h.

Figure 1. Basis function ϕi on its support (xi1,xi+1).

We now start building the stiffness matrix Ah approximating the fractional Laplacian. This will be done in three steps, since the values of the matrix can be computed differentiating among three well defined regions: the upper triangle, corresponding to ji+2, the upper diagonal corresponding to j=i+1 and the diagonal, corresponding to j=i. In each of these regions the intersections among the support of the basis functions are different, thus generating different values of the bilinear form.

We present below an abridged explanation of how to compute the entries of Ah. Complete details may be found in [2].

Step 1: ji+2

In this case we have supp(ϕi)supp(ϕj)= (see also Figure 2). Hence, (5) is reduced to computing only the integral

ai,j=2xj+1xj1xi+1xi1ϕi(x)ϕj(y)|xy|1+2sdxdy.

Figure 2. Basis functions ϕi(x) and ϕj(x) for ji+1. The supports are disjoint.

Taking into account the definition of the basis function (6), the integral (7) becomes

ai,j=2xj+1xj1xi+1xi1(1|xxi|h)(1|yxj|h)|xy|1+2sdxdy.

Let us introduce the following change of variables:

xxih=ˆx,yxjh=ˆy.

Then, rewriting (with some abuse of notations since there is no possibility of confusion) ˆx=x and ˆy=y, we get

ai,j=2h12s1111(1|x|)(1|y|)|xy+ij|1+2sdxdy.

This integral can be computed explicitly, and we obtain

ai,j=h12s4(k+1)32s+4(k1)32s6k32s(k+2)32s(k2)32s2s(12s)(1s)(32s).

Notice that, when s=1/2, both the numerator and the denominator of the expression above are zero. In this case, the value of ai,j can be obtained by taking the limit s1/2 and we get

ai,j=4(k+1)2log(k+1)4(k1)2log(k1)+6k2log(k)+(k+2)2log(k+2)+(k2)2log(k2)

if k2 and

ai,i+2=56ln(2)36ln(3).

Step 2: j=i+1

This is the most cumbersome case, since it is the one with the most interactions between the basis functions (see Figure 3).

Figure 3. Basis functions ϕi(x) and ϕi+1(x). In this case, the intersection of the supports is the interval [xi,xi+1].

Using the symmetry of the integral with respect to the bisector y=x, we have

ai,i+1=RR(ϕi(x)ϕi(y))(ϕi+1(x)ϕi+1(y))|xy|1+2sdxdy=+xi+1+xi+1dxdy+2+xi+1xi+1xidxdy+2+xi+1xidxdy+xi+1xixi+1xidxdy+2xi+1xixidxdy+xixidxdy:=Q1+Q2+Q3+Q4+Q5+Q6.

Also in this case, the terms Qi, i=1,,6, can be computed explicitely, by noticing also the following facts:

  1. Q1=Q6=0 since ϕi=0 on the domain of integration.
  2. Q2=Q5 due to symmetry.

Then, the elements ai,i+1 are given by the sum 2Q2+Q3+Q4 and we have

ai,i+1={h12s332s252s+72s(12s)(1s)(32s),s129ln316ln2,s=12.

Step 3: j=i

As a last step, we fill the diagonal of the matrix Ah, which collects the values corresponding to ϕi(x)=ϕj(x) (see Figure 4).

Figure 4. Basis functions ϕi(x)=ϕj(x).

In this case, we have

ai,i=RR(ϕi(x)ϕi(y))2|xy|1+2sdxdy=+xi+1+xi+1dxdy+2+xi+1xi+1xi1dxdy++xi+1xi1dxdy+xi+1xi1xi+1xi1dxdy+2xi1xi+1xi1dxdy++xi1+xi+1dxdy+xi1xi1dxdy:=R1+R2+R3+R4+R5+R6+R7.

Once again, the terms Ri, i=1,,7 can be computed explicitely, by taking also into account that R1=R3=R6=R7=0 because the corresponding integration domains are all away from the support of the basis functions.

The result of these computations gives the values

ai,i={h12s232s4s(12s)(1s)(32s),s128ln2,s=12.

Step 4: building of Ah

Summarizing, we have the following values for the elements of the stiffness matrix Ah: for s1/2

ai,j=h12s{4(k+1)32s+4(k1)32s6k32s(k+2)32s(k2)32s2s(12s)(1s)(32s),k=ji,k2332s252s+72s(12s)(1s)(32s),j=i+1232s4s(12s)(1s)(32s),j=i.

For s=1/2, instead, we have

ai,j={4(ji+1)2log(ji+1)4(ji1)2log(ji1)+6(ji)2log(ji)+(ji+2)2log(ji+2)+(ji2)2log(ji2),j>i+256ln(2)36ln(3),j=i+2.9ln316ln2,j=i+18ln2,j=i.

Numerical results

We present the numerical simulations corresponding to the algorithm previously described.

First of all, we test numerically the accuracy of our method for the resolution of the elliptic equation (1) by applying it to the following problem

{(d2x)su=1,x(1,1)u0,x(1,1)c.

In this particular case, the unique solution u can be computed exactly and it is given by

u(x)=22sπΓ(1+2s2)Γ(1+s)(1x2)sχ(1,1),

where χ(1,1) indicates the characteristic function of the interval (1,1).

The solution of (9) in the case s=0.1 and N=50 is computed with the following script which emploies the function FEFractionalLaplacian.m ginving the FE discretization of (d2x)s.

s = 0.1;
N = 50;
L = 1;

x = linspace(-L,L,N+2);
x = x(2:end-1);
h = x(2)-x(1);

f = @(x) 1+0*x;
Phi = @(x) 1-abs(x);

F = zeros(N,1);

for i=1:N
    xx = linspace(x(i)-h,x(i)+h,N+1);
    xx = 0.5*(xx(2:end)+xx(1:end-1));
    B1 = f(xx).*(Phi((xx-x(i))/h));
    F(i) = ((2*h)/N)*sum(B1); 
end

A = FEFractionalLaplacian(s,L,N);
sol = A\F;

In Figure 5, we show a comparison for different values of s between the exact solution (10) and the computed numerical approximation.

Figure 5. Real and numerical solution.

One can notice that when s=0.1 (and also for other small values of s), the computed solution is to a certain extent different from the exact solution. Notwithstanding, it is well-known that the notion of trace is not defined for the spaces Hs(1,1) with s1/2. Hence, it is somehow natural that we cannot expect a point-wise convergence in this case.

Furthermore, the accuracy of our algorithm is validated by a simple error analysis.

The computation of the error in the space Hs0(1,1) can be readily done by using the definition of the bilinear form, namely

where have used the orthogonality condition a(v_h,u-u_h)=0 for all v_h \in V_h.

For this particular test, since f\equiv 1 in (-1,1), the problem is therefore reduced to

where the right-hand side can be easily computed, since we have the closed formula

and the term corresponding to \int_{-1}^1u_h can be carried out numerically. Moreover, we have the following convergence result.

Theorem [2] For the solution u of \eqref{WFD} and its FE approximation u_h given by \eqref{WFD}, if h is sufficiently small, the following estimates hold

where \mathcal{C} is a positive constant not depending on h.

In Figure 6, we present the computational errors evaluated for different values of s and h, which can be obtained with the following function.

function [h,e] = error_fl(s,N)

x = linspace(-1,1,N+2);
x = x(2:end-1);
h = x(2)-x(1);

f = @(x) 1+0*x;
Phi = @(x) 1-abs(x);

F = zeros(N,1);

for i = 1:N
    xx = linspace(x(i)-h,x(i)+h,N+1);
    xx = 0.5*(xx(2:end)+xx(1:end-1));
    B1 = f(xx).*(Phi((xx-x(i))/h));
    F(i) = ((2*h)/N)*sum(B1); 
end

A = FEFractionalLaplacian(s,1,N);
sol = A\F;
val = pi/(2^(2*s)*gamma(s+0.5)*gamma(s+1.5));
valnum = h*sum(sol);
e = sqrt(val-valnum);

Figure 6. Computational error in logarithmic scale in terms for different values of s.

The rates of convergence shown are of order (in h) of 1/2, and this convergence rate is maintained also for small values of s. This confirms the accuracy of our approximation. In particular, the behavior shown in Figure 6 is not in contrast with the known theoretical results.

References

[1] G. Acosta, F. Bersetche and J. P. Borthagaray, A short FE implementation for a 2d homogeneous Dirichlet problem of a Fractional Laplacian. Comput. Math. Appl., Vol. 74, No. 4 (2017), pp. 784-816.

[2] U. Biccari and V. Hernández-Santamarı́a, Controllability of a one-dimensional fractional heat equation: theoretical and numerical aspects. IMA J. Math. Control. Inf, to appear, 2018.

[3] U. Biccari, M. Warma and E. Zuazua, Controllability of the one-dimensional fractional heat equation under positivity constraints. Submitted.