Numerical Algorithms for the Forward and Backward Fractional Feynman–Kac Equations

Abstract

The Feynman–Kac equations are a type of partial differential equations describing the distribution of functionals of diffusive motion. The probability density function (PDF) of Brownian functionals satisfies the Feynman–Kac formula, being a Schrödinger equation in imaginary time. The functionals of non-Brownian motion, or anomalous diffusion, follow the fractional Feynman–Kac equation (Carmi et al. in J Stat Phys 141:1071–1092, 2010), where the fractional substantial derivative is involved. Based on recently developed discretized schemes for fractional substantial derivatives (Chen and Deng arXiv:1310.3086), this paper focuses on providing algorithms for numerically solving the forward and backward fractional Feynman–Kac equations; since the fractional substantial derivative is non-local time-space coupled operator, new challenges are introduced compared with the ordinary fractional derivative. Two ways (finite difference and finite element) of discretizing the space derivative are considered. For the backward fractional Feynman–Kac equation, the numerical stability and convergence of the algorithms with first order accuracy are theoretically discussed; and the optimal estimates are obtained. For all the provided schemes, including the first order and high order ones, of both forward and backward Feynman–Kac equations, extensive numerical experiments are performed to show their effectiveness.

Appendix

To prove that (1.4) is equivalent to (1.5), we first introduce some properties of the fractional substantial calculus.

For $0<\alpha <1$, using Lemma 2.3, Lemma 2.7 and Lemma 2.6 of [7], there exists

$$\begin{aligned}&{^s\!}D_t^\alpha [{^s\!}I_t^\alpha P(t)] =P(t); \end{aligned}$$

(6.1)

$$\begin{aligned}&{^s\!}I_t^\alpha [{^s\!}D_t^\alpha P(t)] =P(t)-[{^s\!}D_t^{\alpha -1}P(t)]_{t=0}\frac{ t^{\alpha -1} e^{-\rho U(x) t} }{\varGamma (\alpha )};\end{aligned}$$

(6.2)

$$\begin{aligned}&{^s_c}{D}_t^\alpha P(x,t)= {^s\!}D_t^{\alpha -1}[{^s\!}D^1_t P(t)] ={^s\!}D_t^{\alpha }P(t)-\frac{ t^{-\alpha } e^{-\rho U(x) t} }{\varGamma (1-\alpha )}P(0). \end{aligned}$$

(6.3)

Theorem 5.1

Let $0<\alpha <1$ and $P(x,t) \in C_{x,t}^{2,1}[a,b]\times [0,T]$. Then

$$\begin{aligned} {^s\!}D^1_tP(x,t)= {^s\!}D_t^{1-\alpha }\left[ \kappa _\alpha \frac{\partial ^2}{\partial x^2} P(x,t) \right] , \end{aligned}$$

(6.4)

is equivalent to

$$\begin{aligned} {^s_c}{D}_t^\alpha P(x,t)={^s\!}D_t^\alpha P(x,t)-\frac{t^{-\alpha }e^{-\rho U(x) t}}{\varGamma (1-\alpha )} P(x,0) = \kappa _\alpha \frac{\partial ^2}{\partial x^2} P(x,t). \end{aligned}$$

(6.5)

Proof

Derive (6.5) from (6.4). Performing both sides of (6.4) by ${^s\!}D_t^{\alpha -1}$ leads to

$$\begin{aligned} {^s\!}D_t^{\alpha -1} [{^s\!}D_tP(x,t)] ={^s\!}D_t^{\alpha -1}\left\{ {^s\!}D_t^{1-\alpha }\left[ \kappa _\alpha \frac{\partial ^2}{\partial x^2} P(x,t) \right] \right\} . \end{aligned}$$

According to above equation and (6.3), (6.2), we get

$$\begin{aligned} {^s\!}D_t^{\alpha }P(x,t)-\frac{ t^{-\alpha } e^{-\rho U(x) t}P(x,0) }{\varGamma (1-\alpha )} =\kappa _\alpha \frac{\partial ^2}{\partial x^2} P(x,t) -{^s\!}D_t^{-\alpha }\left[ \kappa _\alpha \frac{\partial ^2}{\partial x^2} P(x,t)\right] _{t=0}\frac{ t^{-\alpha } e^{-\rho U(x) t} }{\varGamma (1-\alpha )}. \end{aligned}$$

If a function $P(t)$ is continuously differentiable in the closed interval $[0,t]$, then

$$\begin{aligned}&{^s\!}D_t^{-\alpha }P(t)|_{t=0}\\&\quad =\lim _{ t\rightarrow 0^{+}}\frac{1}{\varGamma (\alpha )}\int \limits _{0}^t{\left( t-\tau \right) ^{\alpha -1}}e^{-\rho U(x)(t-\tau )}{P(\tau )}d\tau \\&\quad =\lim _{ t\rightarrow 0^{+}}\left[ \frac{t^\alpha e^{-\rho U(x) t}}{\varGamma (\alpha +1)} P(0) +\frac{1}{\varGamma (\alpha +1)}\int \limits _{0}^t{\left( t-\tau \right) ^{\alpha }}e^{-\rho U(x)(t-\tau )}{[\rho U(x) P(\tau )+P'(\tau )]}d\tau \right] \\&\quad =0. \end{aligned}$$

Since $P(x,t) \in C_{x,t}^{2,1}[a,b]\times [0,T]$, we have

$$\begin{aligned} {^s\!}D_t^{-\alpha }\left[ \kappa _\alpha \frac{\partial ^2}{\partial x^2} P(x,t)\right] _{t=0}=0. \end{aligned}$$

It implies that (6.5) holds.

Derive (6.4) from (6.5). Performing both sides of (6.5) by ${^s\!}D_t^{1-\alpha }$ results in

$$\begin{aligned} {^s\!}D_t^{1-\alpha }\left[ {^s\!}D_t^\alpha P(x,t)-\frac{t^{-\alpha }e^{-\rho U(x) t}}{\varGamma (1-\alpha )} P(x,0)\right] ={^s\!}D_t^{1-\alpha }\left[ \kappa _\alpha \frac{\partial ^2}{\partial x^2} P(x,t)\right] . \end{aligned}$$

Using (6.3) and (6.1), there exists

$$\begin{aligned} {^s\!}D_t^{1-\alpha }\left[ {^s\!}D_t^\alpha P(x,t)-\frac{t^{-\alpha }e^{-\rho U(x) t}}{\varGamma (1-\alpha )} P(x,0)\right] ={^s\!}D_t^{1-\alpha }\left[ {^s\!}D_t^{\alpha -1}[{^s\!}D^1_t P(x,t)]\right] = {^s\!}D^1_t P(x,t). \end{aligned}$$

That is (6.4) holds.