 # Error estimates for backward fractional Feynman-Kac equation with non-smooth initial data

In this paper, we are concerned with the numerical solution for the backward fractional Feynman-Kac equation with non-smooth initial data. Here we first provide the regularity estimate of the solution. And then we use the backward Euler and second-order backward difference convolution quadratures to approximate the Riemann-Liouville fractional substantial derivative and get the first- and second-order convergence in time. The finite element method is used to discretize the Laplace operator with the optimal convergence rates. Compared with the previous works for the backward fractional Feynman-Kac equation, the main advantage of the current discretization is that we don't need the assumption on the regularity of the solution in temporal and spatial directions. Moreover, the error estimates of the time semi-discrete schemes and the fully discrete schemes are also provided. Finally, we perform the numerical experiments to verify the effectiveness of the presented algorithms.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The Feynman-Kac equation describes the distribution of the functionals of the trajectories of the particles, where the functional is defined as with being a trajectory of a particle and a prescribed function depending on specific applications Kac1949

. There are two kinds of Feynman-Kac equations: one is for the forward Feynman-Kac equation, governing the joint probability density of the functional and position; and another one is for the backward equation, just focusing on the distribution of the functionals. If the particles are with power-law waiting time and/or jump length distribution(s), the governing equations for the distribution of the functionals are so-called fractional Feynman-Kac equations

Agmon1984 ; Carmi2010 ; Wang2018 , since the fractional substantial derivative is involved in the equations. More generalizations of the Feynman-Kac equations include the models governing the distribution of the functionals of the particles undergoing the reaction and diffusion processes and of the particles with multiple internal states Hou2018 ; Xu2018 ; Xu2018-2 .

Here we solve the following backward fractional Feynman-Kac equation, presented in Carmi2010 , describing the functional distribution of the particles with power-law waiting time, i.e.,

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩∂G(x0,ρ,t)∂t=0D1−α,x0tΔG(x0,ρ,t)−ρU(x0)G(x0,ρ,t),(x0,t)∈Ω×(0,T],G(x0,ρ,0)=G0(x0),x0∈Ω,G(x0,ρ,t)=0,(x0,t)∈∂Ω×(0,T], (1.1)

where and

is the joint probability density function of finding the particle on

at time with the initial position of the particle at ; is the Fourier pair of ; ; stands for Laplace operator; is a bounded domain and is assumed to be bounded in in this paper; is a fixed final time; denotes the Riemann-Liouville fractional substantial derivative, whose definition Li2015 is

 0D1−α,x0tG(x0,ρ,t)= 1Γ(α)[∂∂t+ρU(x0)]∫t0(t−ξ)α−1e−(t−ξ)ρU(x0)G(x0,ρ,ξ)dξ (1.2) = e−tρU(x0) 0D1−αt(etρU(x0)G(x0,ρ,t)),

where denotes the Riemann-Liouville fractional derivative with the definition Podlubny1999

 0DαtG(x0,ρ,t)=1Γ(1−α)∂∂t∫t0(t−ξ)−αG(x0,ρ,ξ)dξ,α∈(0,1).

So far there have been many works for fractional partial differential equations, including the finite difference method, finite element method, spectral method, and so on

Acosta2018 ; Bazhlekova2015 ; chen2009 ; Cheng2015 ; Deng2008 ; Deng2013 ; Ervin2006 ; li2012 ; sun2006 , but there are relatively less researches on solving fractional Feynman-Kac equation numerically Chen2018 ; Deng2015 ; Deng:17 ; Deng2017 ; Nie2019 . The main reasons are that fractional substantial derivative is a time-space coupled non-local operator and the equation covers the complex parameters which bring about many challenges on regularity and numerical analyses. To our best knowledge, numerical approximation on fractional substantial derivative is given in Chen2015 ; Ref. Deng2015 numerically solves the forward and backward fractional Feynman-Kac equations with the assumptions that the solution is regular, is a positive constant, and ( means the real part of ); Ref. Deng2017 presents the error estimate for the backward fractional Feynman-Kac equation with and ; Ref. Deng:17 provides an efficient time-stepping method to solve the forward fractional Feynman-Kac equation and makes error analysis in the measure norm. In this paper, we use the finite element method in space and convolution quadrature introduced in Lubich1988 ; Lubich1988-2 in time to solve the backward fractional Feynman-Kac equation (1.1). The main contributions are as follows.

• We first provide Sobolev regularity for the solution of Eq. (1.1), i.e., Theorem 2.1 gives that the solution when is bounded in and . Compared with the previous works Chen2018 ; Deng2015 ; Deng2017 , we construct numerical scheme without any assumption on the regularity of solution in temporal and spatial directions.

• Then we modify the approximation of the Riemann-Liouville fractional derivative got by convolution quadrature to approximate the Riemann-Liouville fractional substantial derivative, which skillfully overcome the trouble brought by the non-commutativity of the Riemann-Liouville fractional derivative and in error estimate for fully discrete scheme, i.e., in Eq. (3.1).

• Next, a suitable modify based on the Laplace transform representation of solution is presented to guarantee the accuracy of second-order backward difference scheme (3.9) (see Sec. 3).

• Besides, motivated by the error estimate in space in Bazhlekova2015 ; Jin2016 , a general idea is to get the estimate of the difference between and (for their detailed definitions, see Sec. 3 and Sec. 4). Generally, the sufficient regularity on is required to ensure the accuracy of the approximation. Here, we use for the fully discrete scheme (4.1) in Sec. 4 instead of , which weakens the requirement of regularity on to keep the accuracy of the numerical scheme.

• Finally, we provide a complete error analysis for the proposed numerical scheme and obtain the optimal convergence rates in - and -norm.

The rest of the paper is organized as follows. We first provide some preliminaries and a regularity estimate for the solution of Eq. (1.1) in Sec. 2. Section 3 presents the approximation of the Riemann-Liouville fractional substantial derivative by backward Euler and second-order backward difference convolution quadratures and gives the error estimates of the time semi-discrete schemes. In Sec. 4, we use the finite element method to discretize the Laplace operator and provide the error estimate for the fully discrete scheme with the non-smooth initial data. In Sec. 5, we verify the effectiveness of the algorithm by numerical experiments. We conclude the paper with some discussions in the last section.

## 2 Preliminaries

First, we introduce with a zero Dirichlet boundary condition. For any , denote the space with the norm Thomee2006

 ∥v∥2˙Hr(Ω)=∞∑j=1λrj(v,φj)2,

where

are the eigenvalues ordered non-decreasingly and the corresponding eigenfunctions (normalized in the

norm) of operator . Thus , , and . For and , we define sectors and in the complex plane as

 Σθ={z∈C∖{0},|argz|≤θ}, Σθ,κ={z∈C:|z|≥κ,|argz|≤θ},

and the contour is defined by

 Γθ,κ={z∈C:|z|=κ,|argz|≤θ}∪{z∈C:z=re±iθ:r≥κ},

oriented with an increasing imaginary part, where denotes the imaginary unit and . Then we denote as the operator norm from to and define and as and respectively in the following. Throughout this paper, denotes a generic positive constant, whose value may differ at each occurrence; and let arbitrary small.

Similar to the skill used in Chen2018 ; Deng2015 ; Deng2017 , Eq. (1.1) can also be converted into

 ⎧⎪⎨⎪⎩C0Dα,x0tG(x0,ρ,t)=ΔG(x0,ρ,t),(x0,t)∈Ω×(0,T],G(x0,ρ,0)=G0(x0),x0∈Ω,G(x0,ρ,t)=0,(x0,t)∈∂Ω×(0,T], (2.1)

where denotes Caputo fractional substantial derivative defined by Li2015

 C0Dα,x0tG(x0,ρ,t)=e−tρU(x0) C0Dαt(etρU(x0)G(x0,ρ,t)),α∈(0,1),

and means the Caputo fractional derivative with its definition Podlubny1999

 C0DαtG(x0,ρ,t)=1Γ(1−α)∫t0(t−ξ)−α∂∂ξG(x0,ρ,ξ)dξ,α∈(0,1).

Then we recall the Laplace transform for the fractional substantial derivative.

###### Lemma 2.1 (Li2015 )

The Laplace transform of the Riemann-Liouville fractional substantial derivative with is given by

 ˜ 0Dα,x0tG(z)=(β(z,x0))α~G(z),

and the Laplace transform of the Caputo fractional substantial derivative with is given by

 ˜ C0Dα,x0tG(z)=(β(z,x0))α~G(z)−(β(z,x0))α−1G(0),

where and ‘’ stands for taking the Laplace transform. And in the following we denote as .

According to Lemma 2.1, the solution of Eq. can be written as

 ~G(z)=((β(z))α+A)−1(β(z))α−1G0. (2.2)
###### Remark 2.1

By the definition of , it is easy to see that

 β(z)A≠Aβ(z),A((β(z))α+A)−1≠((β(z))α+A)−1A, ((β(z))α+A)−1(β(z))α−1≠(β(z))α−1((β(z))α+A)−1.

Before we provide the regularity estimate for the solution of Eq. (2.1), the following lemma about is also needed.

###### Lemma 2.2 (Deng:17 )

Let be bounded in . By choosing sufficiently close to and sufficiently large (depending on the value ), we have the following results:

1. For all and , we have , and

 C1|z|≤|β(z)|≤C2|z|, (2.3)

where and denote two positive constants. So and are both analytic function of .

2. The operator is well-defined, bounded, and analytic with respect to , satisfying

 ∥A((β(z))α+A)−1∥≤C     for all  z∈Σθ,κ, (2.4)
 ∥((β(z))α+A)−1∥≤C|z|−α   for all  z∈Σθ,κ. (2.5)
###### Theorem 2.1

Assume is bounded in . If and is the solution of Eq. (2.1), then we have the estimate

 ∥G(t)∥˙Hσ(Ω)≤Ct−σα/2∥G0∥L2(Ω),σ∈[0,2].
###### Proof

Taking inverse Laplace transform and norm on both sides of (2.2), according to Lemma 2.2, we have

 ∥G(t)∥L2(Ω)≤ C∥∥ ∥∥∫Γθ,κezt((β(z))α+A)−1(β(z))α−1G0dz∥∥ ∥∥L2(Ω) ≤ C∫Γθ,κ|ezt|∥∥((β(z))α+A)−1(β(z))α−1∥∥∥G0∥L2(Ω)|dz| ≤ C∫Γθ,κ|ezt||z|−1∥G0∥L2(Ω)|dz| ≤ C(∫∞κercos(θ)tr−1dr+∫θ−θeκcos(φ)tdφ)∥G0∥L2(Ω) ≤ C∥G0∥L2(Ω).

Applying on both sides of (2.2), taking inverse Laplace transform, and acting norm on both sides, from Lemma 2.2, there is

 ∥AG(t)∥L2(Ω)≤ C∥∥ ∥∥∫Γθ,κeztA((β(z))α+A)−1(β(z))α−1G0dz∥∥ ∥∥L2(Ω) ≤ C∫Γθ,κ|ezt|∥∥A((β(z))α+A)−1(β(z))α−1∥∥∥G0∥L2(Ω)|dz| ≤ C∫Γθ,κ|ezt||z|α−1∥G0∥L2(Ω)|dz| ≤ C(∫∞κercos(θ)trα−1dr+∫θ−θeκcos(φ)tκαdφ)∥G0∥L2(Ω) ≤ Ct−α∥G0∥L2(Ω).

Using interpolation properties

 ∥G(t)∥˙Hσ(Ω)≤Ct−σα/2∥G0∥L2(Ω),σ∈[0,2].

## 3 Temporal discretization and error analysis

In this section, we first use the backward Euler and second-order backward difference convolution quadratures introduced in Lubich1988 ; Lubich1988-2 to discretize the Riemann-Liouville fractional substantial derivative and get the first- and second-order schemes in time. Then we provide the complete error analysis.

Let the time step size with , , , and . Firstly, we use the relationship between Caputo fractional derivative and Riemann-Liouville fractional derivative

 C0Dαtu(t)=0Dαt(u(t)−u(0))  with  α∈(0,1)

to reformulate Eq. (2.1) with Riemann-Liouville fractional substantial derivative, i.e.,

 ⎧⎪⎨⎪⎩0Dα,x0tG(x0,ρ,t)+AG(x0,ρ,t)=e−ρU(x0)t0DαtG(x0,ρ,0),(x0,t)∈Ω×(0,T],G(x0,ρ,0)=G0(x0),x0∈Ω,G(x0,ρ,t)=0,(x0,t)∈∂Ω×(0,T]. (3.1)

### 3.1 Backward Euler scheme and error estimate

We use backward Euler convolution quadrature to discretize the time fractional substantial derivative and get the first-order accuracy in time. Introduce as the numerical approximation of solution . Then we can obtain the temporal semi-discrete scheme

 ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩n−1∑i=0dα,1ie−tiρU(x0)Gn−i+AGn=e−tnρU(x0)n−1∑i=0dα,1iG0,x0∈Ω,n=1,…,N,G0=G0,x0∈Ω,Gn(x0)=0,x0∈∂Ω,n=1,…,N, (3.2)

where

 (δτ,1(ζ))α=∞∑i=0dα,1iζi (3.3)

and

 δτ,1(ζ)=1−ζτ.

Multiplying on both sides of the first formula of (3.2) and summing from to lead to

 ∞∑n=1n−1∑i=0dα,1ie−tiρU(x0)Gn−iζn+∞∑n=1AGnζn=∞∑n=1e−tnρU(x0)n−1∑i=0dα,1iG0ζn. (3.4)

Simple calculation implies

 (δτ,1(e−τρU(x0)ζ))α∞∑n=1Gnζn+A∞∑n=1Gnζn=(δτ,1(e−τρU(x0)ζ))α∞∑n=1e−tnρU(x0)G0ζn,

which is followed by . Furthermore, we have

 (δτ,1(e−τρU(x0)ζ))α∞∑n=1Gnζn+A∞∑n=1Gnζn=(δτ,1(e−τρU(x0)ζ))α−1G0e−τρU(x0)ζτ.

Using Cauchy’s integral formula yields

 Gn= 12πi∫|ζ|=ξτζ−n−1((δτ,1(e−τρU(x0)ζ))α+A)−1(δτ,1(e−τρU(x0)ζ))α−1G0e−τρU(x0)ζτdζ (3.5) = 12πi∫Γτeztn((βτ,1(z))α+A)−1(βτ,1(z))α−1e−τ(z+ρU(x0))G0dz,

where , , , and the second equality follows by taking . Deforming the contour to , one has

 Gn=12πi∫Γτθ,κeztn((βτ,1(z))α+A)−1(βτ,1(z))α−1e−τ(z+ρU(x0))G0dz. (3.6)

Next, we provide a lemma about defined in (3.5).

###### Lemma 3.1 (Deng:17 )

Let be bounded in . By choosing sufficiently close to and sufficiently large depending on , there exists a positive constant depending on and such that the following estimates hold when :

1. For , we have , and

 C1|z|≤|βτ,1(z)|≤C2|z|.
2. The operator is well-defined, bounded, and analytic with respect to , satisfying

 ∥A((βτ,1(z))α+A)−1∥≤C     for all  z∈Στθ,κ,
 ∥((βτ,1(z))α+A)−1∥≤C|z|−α     for all  z∈Στθ,κ,

where . Here, means the imaginary part of and the real part of .

3. For the real number , the following estimate holds

 |(β(z))γ−(βτ,1(z))γ|≤Cτ|z|γ+1,z∈Γτθ,κ.
###### Theorem 3.1

Let and be the solutions of Eqs. (2.1) and (3.2) respectively and assume . Then we obtain

 ∥G(x0,ρ,tn)−Gn∥L2(Ω)≤Ct−1nτ∥G0∥L2(Ω).
###### Proof

Subtracting (3.6) from the inverse Laplace transform of (2.2), we have

 ∥G(x0,ρ,tn)−Gn∥L2(Ω) ≤ C∥∥ ∥∥∫Γθ,κeztn((β(z))α+A)−1(β(z))α−1G0dz −∫Γτθ,κeztn((βτ,1(z))α+A)−1(βτ,1(z))α−1e−τ(z+ρU(x0))G0dz∥∥ ∥∥L2(Ω) ≤ C∥∥ ∥∥∫Γθ,κ∖Γτθ,κeztn((β(z))α+A)−1(β(z))α−1dz∥∥ ∥∥∥G0∥L2(Ω) +C∥∥∫Γτθ,κeztn(((β(z))α+A)−1(β(z))α−1 −((βτ,1(z))α+A)−1(βτ,1(z))α−1e−τ(z+ρU(x0)))dz∥∥∥G0∥L2(Ω) ≤ C(I+II)∥G0∥L2(Ω).

For , using Lemma 2.2, there is

 I≤ ∫Γθ,κ∖Γτθ,κ|eztn||z|−1|dz|≤Cτ∫∞πτsin(θ)etnrcos(θ)dr≤Ct−1nτ.

As for , one can split it into

 II≤ ∥∥ ∥∥∫Γτθ,κeztn(((β(z))α+A)−1(β(z))α−1−((βτ,1(z))α+A)−1(β(z))α−1)dz∥∥ ∥∥ +∥∥ ∥∥∫Γτθ,κeztn(((βτ,1(z))α+A)−1(β(z))α−1−((βτ,1(z))α+A)−1(βτ,1(z))α−1)dz∥∥ ∥∥ +∥∥ ∥∥∫Γτθ,κeztn((βτ,1(z))α+A)−1(βτ,1(z))α−1(1−e−τ(z+ρU(x0)))dz∥∥ ∥∥.

Then by Lemmas 2.2, 3.1 and the fact

 ∥∥((β(z))α+A)−1(β(z))α−1−((βτ,1(z))α+A)−1(β(z))α−1∥∥ = ∥∥((β(z))α+A)−1((βτ,1(z))α−(β(z))α)((βτ,1(z))α+A)−1(β(z))α−1∥∥≤Cτ,

it has

 II≤ Cτ∫Γτθ,κ|eztn||dz|≤Cτt−1n.

Thus

 ∥G(x0,ρ,tn)−Gn∥L2(Ω)≤Ct−1nτ∥G0∥L2(Ω).

### 3.2 Second-order backward difference scheme and error estimate

In this subsection, we use second-order backward difference convolution quadrature to discretize the time fractional substantial derivative and obtain the second-order accuracy in time. Similarly, introduce as the numerical approximation of the solution , and let

 δτ,2(ζ)=(1−ζ)+(1−ζ)2/2τandν(ζ)=(3−ζ2(1−ζ))ζ=ζ(32+∞∑n=1ζn). (3.7)

According to (2.2), we have

 ~G(z)−(β(z))−1G0=−((β(z))α+A)−1A(β(z))−1G0. (3.8)

Using , , and to approximate , , and respectively, we have

 ∞∑n=1(Gn−e−tnρU(x0)G0)e−ztn=−((δτ,2(e−β(z)