 # Correction of BDFk for fractional Feynman-Kac equation with Lévy flight

In this work, we present the correction formulas of the k-step BDF convolution quadrature at the starting k-1 steps for the fractional Feynman-Kac equation with Lévy flight. The desired kth-order convergence rate can be achieved with nonsmooth data. Based on the idea of [Jin, Li, and Zhou, SIAM J. Sci. Comput., 39 (2017), A3129–A3152], we provide a detailed convergence analysis for the correction BDFk scheme. The numerical experiments with spectral method are given to illustrate the effectiveness of the presented method. To the best of our knowledge, this is the first proof of the convergence analysis and numerical verified the sapce fractional evolution equation with correction BDFk.

## 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

Functionals of Brownian motion have diverse applications in physics, mathematics, and other fields. The probability density function of Brownian functionals satisfies the Feynman-Kac formula, which is a Schrödinger equation in imaginary time. The functionals of non-Brownian motion, or anomalous diffusion, follow the general fractional Feynman-Kac equation

[1, 2], where the fractional substantial derivative is involved . This paper focuses on providing the correction of -step backward differential formulas (BDF) for backward fractional Feynman-Kac equation with Lévy flight [1, 2, 5, 6]

 (1) CsDγtG(t)+AG(t)=f(t)  with  A:=(−Δ)α/2,

where is a given function and the initial condition with the homogeneous Dirichlet boundary conditions. Here with is the fractional Laplacian, the definition of which is based on the spectral decomposition of the Dirichlet Laplacian [20, 32]; and the Caputo fractional substantial derivative with is defined by [4, 5, 11]

 (2) CsDγtG(t)=1Γ(1−γ)∫t0e−σ(t−s)(t−s)γ(σ+∂∂s)G(s)ds

with a constant . It should be noted that reduces to the Caputo fractional derivative if .

High-order schemes for the time discretization of (1) with (Caputo fractional derivative) have been proposed by various authors. There are two predominant discretization techniques in time direction: L1-type approximation [19, 23, 25, 31] and Lubich-Grünwald-Letnikov approximation [3, 21, 24]. For the first group, under the time regularity assumption on the solution, they developed the L1 schemes [19, 31] for the Caputo fractional derivative and strictly proved the stability and convergence rate with

. It is extended to the quadratic interpolation case

 with a convergence rate . Recently, for a layer or blows up at , a sharp new discrete stability result has been considered in 

. In the second group, using fractional linear multistep method and Fourier transform, error analysis of up to sixth order temporal accuracy for fractional ordinary differential equation has been discussed

 with the starting quadrature weights schemes. A few years later, based on operational calculus with sectorial operator, nonsmooth data error estimates for fractional evolution equations have been studied in [8, 22] and developed in [16, 17] to restore . Under the time regularity assumption, high order finite difference method (BDF2) for the anomalous-diffusion equation has been studied in  by analyzing the properties of the coefficients. Application of Grenander-Szegö theorem, stability and convergence for time-fractional sub-diffusion equation have been provided in [13, 15] with weighted and shifted Grünwald operator.

In recent years, the numerical method for backward fractional Feynman-Kac equation (1) with were developed. For example, the time discretization of Caputo fractional substantial derivative was first provided in  with the starting quadrature weights schemes. Spectral methods for substantial fractional ordinary differential equations was presented in . Under smooth assumption, numerical algorithms (finite difference and finite element) for (1) with are considered in . Moreover, the second-convergence analysis are discussed in [6, 12]. In addition, the problem with nonsmooth solution is also discussed in . Recently, the second-order error estimates are presented in  with nonsmooth initial data. However, it seems that there are no published works for more than three order accurate scheme for model (1) with Lévy flight. In this work, we provide a detailed convergence analysis of the correction BDF () for (1) with nonsmooth data.

We first provide the solutions of fractional Feynman-Kac equation with Lévy flight.

### Solution representation for (1)

Let be a bounded domain with a boundary . If , then (1) reduce to the following time-space Caputo-Riesz fractional diffusion equation 

 ⎧⎪⎨⎪⎩C∂γtG(x,t)+AG(x,t)=f(x,t),(x,t)∈Ω×(0,T]G(x,0)=v(x),x∈ΩG(x,t)=0,(x,t)∈∂Ω×(0,T],

where is a given function and denotes the Laplacian . Based on the idea of [16, 27, 32] with the eigenpairs of the operator , it is easy to get

 (3) G(x,t)=E(t)v+∫t0¯¯¯¯E(t−s)f(s)ds.

Here the operators and are, respectively, given by

 E(t)=∞∑j=1Eγ,1(−λα/2jtγ)(v,φj)φj(x)

and

 ¯¯¯¯E(t)χ=∞∑j=1tγ−1Eγ,γ(−λα/2jtγ)(v,φj)(χ,φj)φj(x)

with the Mittag-Leffler function [26, p. 17], i.e.,

 Eγ,α/2(z)=∞∑k=0zkΓ(kγ+α/2),z∈C.

If , using (3) and the following property [4, 5]

 CsDγte−σtG(t)=e−σt(C∂γtG(t)),

infer that

 (4) G(x,t)=e−σtE(t)v+e−σt∫t0¯¯¯¯E(t−s)eσsf(s)ds.

From the above solution representation of (1), we known that the smoothness of all the data of (1) do not imply the smoothness of the solution . For example, if and with , the following estimate holds [27, Theorem 2.1]

 ∥C∂γtG(t)∥L2(Ω)≤ct−γ∥G0∥L2(Ω),

which reduces to a classical case if [32, Lemma 3.2 ]. This shows that has an initial layer at (i.e., unbounded near ) . Hence, the high-order convergence rates may not hold for nonsmooth data. Thus, the corrected algorithms are necessary in order to restore the desired convergence rate, even for smooth initial data. In this paper, based on the idea of , we present the correction of the -step BDF convolution quadrature (CQ) at the starting steps for the backward fractional Feynman-Kac equation with Lévy flight (1). The desired th-order convergence rate can be achieved with nonsmooth data. To the best of our knowledge, this is the first proof of the convergence analysis and numerical verified the sapce fractional evolution equation with correction BDF.

The paper is organized as follows. In the next Section, we provide the correction of the -step BDF convolution quadrature at the starting steps for (1). In Section 3, based on operational calculus, the detailed convergence analysis of the correction BDF are provided. To show the effectiveness of the presented schemes, the results of numerical experiments are reported in Section 4.

## 2 Correction of BDFk

Let with the uniform time steplength, and let denote the approximation of and . The convolution quadrature generated by BDF, approximates the Riemann-Liouville fractional substantial derivative by 

 (5) ¯¯¯¯¯Dγτφn:=1τγn∑j=0qjφn−j

with . Here the weights and are the coefficients in the series expansion

 (6) δγτ(ξ)=1τγ∞∑j=0bjξjwithδτ(ξ):=1τk∑j=11j(1−ξ)jand  δ(ξ):=δ1(ξ).

Then the standard BDF for (1) is as following

 (7) ¯¯¯¯¯Dγτ(Gn−e−σtnG(0))+AGn=f(tn).

To obtain the -order accuracy with nosmooth data, we correct the standard BDF (7) at the starting steps by

 (8) ¯¯¯¯¯Dγτ(Gn−e−σtnG(0))+AGn=−a(k)ne−σnτAG0+b(k)nf(0)+k−2∑l=1d(k)l,nτl∂ltf(0)+f(tn)1≤n≤k−1;¯¯¯¯¯Dγτ(Gn−e−σtnG(0))+AGn=f(tn)k≤n≤N.

Here the correction coefficients and are given in Table 1 and is given in Table 2. We noted that the correction coefficients share similarities with the fractional Caputo equations .

### 2.1 Solution representation with CQ for (1)

First we split right-hand side into

 (9) f(t)=f(0)+k−2∑l=1tll!∂ltf(0)+Rk.

Here

 (10) Rk=f(t)−f(0)−k−2∑l=1tll!∂ltf(0)=tk−1(k−1)!∂k−1tf(0)+tk−1(k−1)!∗∂ktf,

and the symbol denotes Laplace convolution. Let with . Then we can rewrite (1) as

 (11) sDγtW(t)+AW(t)=−Ae−σtG(0)+f(0)+k−2∑l=1tll!∂ltf(0)+Rk

with in .

Applying Laplace transform in  to both sides of the above equation, we have

 (σ+z)γˆW(z)+AˆW(z)=−A(σ+z)−1G(0)+z−1f(0)+k−2∑l=11zl+1∂ltf(0)+ˆRk(z),

where and denote the Laplace transform, i.e, . Then

 ˆW(z)=((σ+z)γ+A)−1[−A(σ+z)−1G(0)+z−1f(0)+k−2∑l=11zl+1∂ltf(0)+ˆRk(z)].

By the inverse Laplace transform, we can obtain the solution as following

 (12) W(t)=12πi∫Γθ,κeztK(σ+z)(−AG(0)+σ+zzf(0))dz+12πi∫Γθ,κezt(σ+z)K(σ+z)(k−2∑l=11zl+1∂ltf(0)+ˆRk(z))dz

with

 (13) K(σ+z)=((σ+z)γ+A)−1(σ+z)−1.

Here the is defined by 

 (14) Γθ,κ={z∈C:|z|=κ,|argz|≤θ}∪{z∈C:z=re±iθ,κ≤r<∞}.

It should be noted that we have the following resolvent bound 

 (15) ||(zγ+A)−1||≤c|z|−γand||K(z)||≤c|z|−1−γ

with a positive constant .

### 2.2 Discrete solution representation with CQ for (8)

In this subsection, we provide the discrete solution of (8). Let and . Let discrete solution with . Then

 (16) Wn=12πi∫Γτθ,κetnzK(δτ(e−(σ+z)τ))[−μ1(e−(σ+z)τ)AG0+δτ(e−(σ+z)τ)δτ(e−zτ)μ2(e−zτ)f(0)]dz+12πi∫Γτθ,κetnzK(δτ(e−(σ+z)τ))δτ(e−(σ+z)τ)×k−2∑l=1(γl(e−zτ)l!+k−1∑j=1d(k)l,ne−zjτ)τl+1∂ltf(0)dz+12πi∫Γτθ,κetnzK(δτ(e−(σ+z)τ))δτ(e−(σ+z)τ)τ˜Rk(e−zτ)dz.

Here the contour is

 Γτθ,κ={z∈C:|z|=κ,|argz|≤θ}∪{z∈C:z=re±iθ,κ≤r<πτsin(θ)},

and , and

 (17) μ1(ξ)=δ(ξ)(ξ1−ξ+k−1∑j=1a(k)jξj),μ2(ξ)=δ(ξ)(ξ1−ξ+k−1∑j=1b(k)jξj)  and  γl(ξ)=∞∑n=1nlξn.

Let . From (8), it holds

 (18) ¯¯¯¯¯DγτWn+AWn=−(1+a(k)n)e−σnτAG0+(1+a(k)n)f(0)+k−2∑l=1(tlnl!+d(k)l,nτl)∂ltf(0)+Rk(tn)1≤n≤k−1;¯¯¯¯¯DγτWn+AWn=−e−σnτAG0+f(0)+k−2∑l=1tlnl!∂ltf(0)+Rk(tn)k≤n≤N.

Multiplying (18) by and summing over , we have

 (19) ∞∑n=1ξn(¯¯¯¯¯DγτWn+AWn)=−(∞∑n=1ξne−σnτ+k−1∑j=1ξja(k)je−σjτ)AG0+(∞∑n=1ξn+k−1∑j=1ξjb(k)j)f(0)+k−2∑l=1(∞∑n=1ξntlnl!+k−1∑j=1ξjd(k)l,nτl)∂ltf(0)+˜Rk(ξ)=−(ξe−στ1−ξe−στ+k−1∑j=1ξja(k)je−σjτ)AG0+(ξ1−ξ+k−1∑j=1ξjb(k)j)f(0)+k−2∑l=1(γl(ξ)l!+k−1∑j=1ξjd(k)l,n)τl∂ltf(0)+˜Rk(ξ),

where and . Here we use the property . Since

 (20) ∞∑n=1ξn¯¯¯¯¯DγτWn=∞∑n=1ξn1τγn∑j=1qn−jWj=∞∑j=1∞∑n=jξn1τγqn−jWj=∞∑j=1∞∑n=0ξn+j1τγqnWj=1τγ∞∑n=0ξnqn∞∑j=1ξjWj=δγτ(e−στξ)˜W(ξ)

with and . We obtain

 (21) ˜W(ξ)=K(δτ(e−στξ))δτ(e−στξ)[−(e−στξ1−e−στξ+k−1∑j=1a(k)je−σjτξj)AG0+(ξ1−ξ+k−1∑j=1b(k)jξj)f(0)+k−2∑l=1(γl(ξ)l!+k−1∑j=1d(k)l,nξj)τl∂ltf(0)+˜Rk(ξ)]=K(δτ(e−στξ))[−τ−1μ1(e−στξ)AG0+τ−1δτ(e−στξ)δτ(ξ)μ2(ξ)f(0)+δτ(e−στξ)k−2∑l=1(γl(ξ)l!+k−1∑j=1d(k)l,nξj)τl∂ltf(0)+δτ(e−στξ)˜Rk(ξ)],

where is given by (13) and

 (22) μ1(e−στξ)=δ(e−στξ)(e−στξ1−e−στξ+k−1∑j=1a(k)je−σjτξj),μ2(ξ)=δ(ξ)(ξ1−ξ+k−1∑j=1b(k)jξj).

According to Cauchy’s integral formula, and the change of variables , and Cauchy’s theorem of complex analysis, we have

 (23) Wn=12πi∫|ξ|=ϱκξ−n−1˜W(ξ)dξ=12πi∫Γτetnzτ˜W(e−zτ)dz=12πi∫Γτθ,κetnzτ˜W(e−zτ)dz

with the Bromwich contours 

 Γτ={z=κ+1+iy:y∈Rand|y|≤π/τ}.

The proof is completed.

## 3 Convergence analysis

In this section, based on the idea of , we provided the detailed convergence analysis of the correction BDF for (1) with Lévy flight.

### 3.1 A few technical Lemmas

We first introduce a few technical lemmas. Let is given by (6) for . Then

 δ(e−y)=y−1k+1yk+1+O(yk+2).

According to (6) and the Taylor series expansion , we have

The proof is completed.

Let be given by (6) for . Then there exist the positive constants and such that

 c1|z|≤|δτ(e−zτ)|≤c2|z|and|δγτ(e−zτ)−zγ|≤cτk|z|k+γ∀z∈Γτθ,κ.

From Lemma 3.1, we have . On the other hand, according to (39) of  and Lemma 3.1, we have

 |δγτ(e−zτ)−zγ|=γ∣∣ ∣∣∫δτ(e−zτ)zζγ−1dζ∣∣ ∣∣≤maxζ|ζ|γ−1∣∣δτ(e−zτ)−z∣∣≤cτk|z|k+γ.

The proof is completed.

Let and be given by (6) and (13), respectively. Then there exist a positive constants such that

 ∥∥K(δτ(e−(σ+z)τ))−K(σ+z)∥∥≤cτk|σ+z|k−1−γ≤cτk|z|k−1−γ,

where in (14) is large enough compared with . Using the triangle inequality and (15), we have

 ∥∥K(δτ(e−(σ+z)τ))−K(σ+z)∥∥=∥∥(δτ(e−(σ+z)τ))−1(δγτ(e−(σ+z)τ)+A)−1−(σ+z)−1((σ+z)γ+A)−1∥∥≤∣∣((δτ(e−(σ+z)τ))−1−(σ+z)−1)∣∣∥∥(δγτ(e−(σ+z)τ)+A)−1∥∥+∣∣(σ+z)−1∣∣∥∥∥(δγτ(e−(σ+z)τ)+A)−1−((σ+z)γ+A)−1∥∥∥≤cτk|σ+z|k−1−γ.

Here we use 

 ∥∥(δγτ(e−(σ+z)τ)+A)−1∥∥≤c|σ+z|−γ

and

 (δγτ(e−(σ+z)τ)+A)−1−((σ+z)γ+A)−1=((σ+z)γ−δγτ(e−(σ+z)τ))(δγτ(e−(σ+z)τ)+A)−1((σ+z)γ+A)−1.

The proof is completed.

Let be given by (6) with and be given by (17). Let , and be given in Table 1 and Table 2. Then

 ∣∣δτ(e−(σ+z)τ)−(σ+z)∣∣≤cτk|σ