Finite Difference Methods for the generator of 1D asymmetric alpha-stable L\'{e}vy motions

TL;DR

This paper develops finite difference methods to approximate the generator of 1D asymmetric alpha-stable Levy motions, facilitating numerical analysis of related stochastic processes.

Contribution

It introduces spectral-based finite difference schemes with specific convolution weights for the generator of asymmetric alpha-stable Levy motions, analyzing their accuracy and practical applicability.

Findings

Methods effectively approximate the generator in spectral space.

Different schemes have distinct accuracy and computational advantages.

Guidance provided for selecting suitable schemes for practical problems.

Abstract

Several finite difference methods are proposed for the infinitesimal generator of 1D asymmetric $\alpha$-stable L\'{e}vy motions, based on the fact that the operator becomes a multiplier in the spectral space. These methods take the general form of a discrete convolution, and the coefficients (or the weights) in the convolution are chosen to approximate the exact multiplier after appropriate transform. The accuracy and the associated advantages/disadvantages are also discussed, providing some guidance on the choice of the right scheme for practical problems, like in the calculation of mean exit time for random processes governed by general asymmetric $\alpha$-stable motions.

Equations

$$\psi(\xi)=i\mu\xi-\frac{1}{2}q^{2}\xi^{2}+\int_{\mathbb{R}\setminus\{0\}}\big{(}e^{iy\xi}-1-iy\xi\mathbf{1}_{|y|\leq 1}\big{)}\nu({\rm d}y),$$

$${\int }_{\mathbb{R}\setminus \{0\}}\min (1,y^2)\nu (\mathrm{d}y)<\infty .$$

$$\mathcal{A}f(x):=\lim_{t\to 0^{+}}\frac{\mathbb{E}\big{(}f(X_{t})\mid X_{0}=x\big{)}-f(x)}{t}.$$

$$\displaystyle\mathcal{A}f(x)=\mu f^{\prime}(x)+\frac{q^{2}}{2}f^{\prime\prime}(x)+\int_{\mathbb{R}\setminus\{0\}}\big{[}f(x+y)-f(x)-y\mathbf{1}_{|y|\leq 1}f^{\prime}(x)\big{]}\ \nu({\rm d}y).$$

$$\widehat{\mathcal{A}f}(\xi )=\psi (-\xi )\hat{f}(\xi ).$$

$$\mathcal{A}u=-1\text{mbox}on\Omega ,$$

$$\nu (\mathrm{d}y)=\frac{C_{+}{\chi }_{(0,\infty )}(y)+C_{-}{\chi }_{(-\infty ,0)}(y)}{|y{|}^{1+\alpha }}(\mathrm{d}y),$$

$$\psi (t)=i(\mu +K)\xi -\frac{1}{2}{q}^2{\xi }^2+\Gamma (-\alpha )C_{+}(-i\xi {)}^{\alpha }+\Gamma (-\alpha )C_{-}(i\xi {)}^{\alpha }$$

$$K=\{\begin{array}{ll}(C_{-}-C_{+})/(1-\alpha ),& \alpha \ne 1,\\ (1-\gamma )(C_{+}-C_{-}),& \alpha =1,\end{array}$$

$$\psi(\xi)=i(\mu+K)\xi-\frac{1}{2}q^{2}\xi^{2}-\sigma^{\alpha}|\xi|^{\alpha}\big{(}1-i\beta\mbox{sgn}(\xi)\tan{({\pi\alpha}/2)}\big{)},$$

$$C_{+}=C_{\alpha }\frac{1+\beta }{2},C_{-}=C_{\alpha }\frac{1-\beta }{2},-1\le \beta \le 1,$$

$$C_{\alpha }=-\frac{{\sigma }^{\alpha }}{\Gamma (-\alpha )\cos (\pi \alpha /2)}.$$

$$\widehat{\mathcal{A}}=\psi(-\xi)=\begin{cases}-|\xi|^{\alpha}\big{(}1+i\beta\operatorname{sgn}(\xi)\tan({\alpha\pi}/{2})\big{)},\qquad&\alpha\neq 1,\cr-|\xi|\big{(}1+i2\beta\operatorname{sgn}(\xi)\ln|\xi|/\pi\big{)},&\alpha=1,\end{cases}$$

$$\hat{v}(\xi )=\mathcal{F}[v](\xi )=h\sum _{j=-\infty }^{\infty }{e}^{-i\xi x_j}{v}_j,\xi \in I_h:=\left[-\frac{\pi }{h},\frac{\pi }{h}\right].$$

$$v_j={\mathcal{F}}^{-1}[\hat{v}](x_j)=\frac{1}{2\pi }{\int }_{-\pi /h}^{\pi /h}{e}^{i\xi x_j}\hat{v}(\xi )\mathrm{d}\xi .$$

$$\mathcal{\ell}^{2}(\mathbb{Z}_{h})=\Big{\{}v:\mathbb{Z}_{h}\to\mathbb{C}\mid h\sum_{j=-\infty}^{\infty}|v_{j}|^{2}<\infty\Big{\}},$$

$$L^{2}(I_{h})=\Big{\{}\hat{v}:I_{h}\to\mathbb{C}\mid\int_{-\pi/h}^{\pi/h}|\hat{v}(\xi)|^{2}{\rm d}\xi<\infty\Big{\}}.$$

$$⟨u,v{⟩}_{{\ell }^2({\mathbb{Z}}_h)}=\frac{1}{2\pi }⟨\hat{u},\hat{v}{⟩}_{L^2(I_h)},$$

$$⟨u,v{⟩}_{{\ell }^2({\mathbb{Z}}_h)}=h\sum _{j=-\infty }^{\infty }\overline{u_j}{v}_j,⟨\hat{u},\hat{v}{⟩}_{L^2(I_h)}={\int }_{-\pi /h}^{\pi /h}\overline{\hat{u}(\xi )}\hat{v}(\xi )\mathrm{d}\xi .$$

$$(\mathcal{D}v{)}_j=h\sum _{k=-\infty }^{\infty }{v}_{j-k}{w}_k$$

$$\widehat{\mathcal{D}v}(\xi )=\hat{w}(\xi )\hat{v}(\xi ).$$

$${\mathcal{A}}_h{u}_j=\sum _{k=-\infty }^{\infty }{w}_{j-k}{u}_k,$$

$${\psi }_h(\xi )=-{\tilde{M}}_{eh}(\xi )+i\beta {\tilde{M}}_{oh}(\xi )\mathrm{sgn}(\xi )\tan (\alpha \pi /2),$$

$$\hat{w}(\xi )=h{\psi }_h(-\xi )=-h{\tilde{M}}_{eh}(\xi )-ih\beta {\tilde{M}}_{oh}(\xi )\mathrm{sgn}(\xi )\tan (\alpha \pi /2)$$

$$w_j$$

$$=-\frac{h^{-\alpha }}{2\pi }{\int }_{-\pi }^{\pi }{e}^{ij\xi }\left[h^{\alpha }{\tilde{M}}_{eh}\left(\frac{\xi }{h}\right)+i\beta h^{\alpha }{\tilde{M}}_{oh}\left(\frac{\xi }{h}\right)\mathrm{sgn}(\xi )\tan \left(\frac{\alpha \pi }{2}\right)\right]\mathrm{d}\xi .$$

$$w_j$$

$$=-\frac{h^{-\alpha }}{\pi }{\int }_0^{\pi }{M}_e(\xi )\cos (j\xi )\mathrm{d}\xi +\frac{h^{-\alpha }}{\pi }\beta \tan \left(\frac{\alpha \pi }{2}\right){\int }_0^{\pi }{M}_o(\xi )\sin (j\xi )\mathrm{d}\xi .$$

$$\sum_{k=-\infty}^{\infty}w_{k}=\left.\left(\sum_{k=-\infty}^{\infty}e^{-i\xi x_{k}}w_{k}\right)\right|_{\xi=0}=\frac{\hat{w}(0)}{h}=\left.-\Big{(}\tilde{M}_{eh}(\xi)+i\beta\tilde{M}_{oh}(\xi)\operatorname{sgn}(\xi)\tan({\alpha\pi}/{2})\Big{)}\right|_{\xi=0}=0.$$

$${\mathcal{A}}_h{u}_j=\sum _{k=-\infty }^{\infty }{w}_{j-k}(u_k-u_j).$$

$$\hat{w}(\xi )=-h^{1-\alpha }{M}_e(h\xi )-i{h}^{1-\alpha }\beta M_o(h\xi )\mathrm{sgn}(\xi )\tan \left(\frac{\alpha \pi }{2}\right),$$

Full text

*AMS Mathematics Subject Classification 2000. * 35R09, 60G51, 65N06

Keywords: Finite difference method, $\alpha$-stable Lévy processes, psudo-differential operator

1 Introduction

While most random dynamical systems studied during the past century are driven by Gaussian noise, nowadays many complex systems are modelled with the presence of non-Gaussian noise generated from even more exotic Lévy processes [23].

A Lévy process is a stochastic process with stationary and independent increments [1, 24]. As a generalisation of the well-known Wiener process (also called Brownian motion), a (scalar) Lévy process $\{X_{t}\}_{t\geq 0}$ is characterised by the Lévy-Khintchine theorem, which states that the characteristic function of $X_{t}$ can be written as \mathbb{E}\big{(}\exp(i\xi X_{t})|X_{0}=0\big{)}=\exp\big{(}t\psi(\xi)\big{)}. Here, $\psi$ takes the general form

[TABLE]

where $\mu\in\mathbb{R},q\geq 0$, $\mathbf{1}_{S}$ is the indicator function on the set $S$ and $\nu$ is the so called Lévy measure that satisfies the condition

[TABLE]

Intimately related to a Lévy process is its infinitesimal generator $\mathcal{A}$, defined for a smooth function $f$ as

[TABLE]

If $X_{t}$ is a Lévy process starting from $x$ at time $t=0$, then \mathbb{E}\big{(}\exp(i\xi(X_{t}-x))|X_{0}=x\big{)}=\exp(t\psi(\xi)), and the generator $\mathcal{A}f$ can also be represented as [5]

[TABLE]

If the Fourier transform of $f$ is defined to be $\hat{f}(\xi)=\mathcal{F}[f](\xi)=\int_{\mathbb{R}}f(x)\exp(-i\xi x){\rm d}\xi$, then $\mathcal{A}$ becomes a multiplier in the spectral space111Notice the minus sign in $\psi(-\xi)$, because of the factor $\exp(-i\xi x)$ used in the definition of the Fourier transform and $\exp(i\xi x)$ is used in the inverse transform $f(x)=\frac{1}{2\pi}\int_{\mathbb{R}}\hat{f}(\xi)\exp(i\xi x)\mathrm{d}x$., i.e.,

[TABLE]

This operator $\mathcal{A}$ plays a similar role as the adjoint of the classical Fokker-Plank operator associated with systems driven by Brownian motions, and appears in the governing non-local differential equations for quantities that characterise the underlying stochastic processes. For instance, the mean exit time $u(x)$ for the expect time of a particle starting at $x$ and leaving the domain $\Omega$ satisfies the equation

[TABLE]

with the boundary condition $u=0$ on $\Omega^{c}$. The adjoint operator $\mathcal{A}^{*}$ also arises in other contexts, like the evolution of the probability density distribution and escape probability [5].

In the modelling of practical problems driven by non-Gaussian noise, the final Lévy process usually appears as the limit of many small Lévy processes [24], giving rise to $\alpha$-stable Lévy processes that enjoy better scaling properties. In these new processes, the Lévy measure is simply

[TABLE]

for some non-negative constants $C_{+}$ and $C_{-}$, and the stability parameter $\alpha\in(0,2)$. Using special integrals (22) in Appendix B, $\psi(t)$ defined in (1) becomes

[TABLE]

where

[TABLE]

and $\gamma(\approx 0.57721)$ is the Euler-Mascheroni constant. By isolating the real and imaginary parts, the function $\psi(t)$ in (4) can be written in the following form commonly seen in the literature

[TABLE]

with the scaling parameter $\sigma>0$, the skewness parameter $\beta\in[-1,1]$. The constants in (4) and (5) are related by

[TABLE]

with

[TABLE]

This connection between (3) and (4) turns out to motivate a particular type of finite difference schemes originated from fractional derivatives [9, 10].

Parallel to the progress made in the theoretical study of Lévy processes [1], related fields have been heavily investigated during the past two decades, including random simulations of these more exotic processes [16] and numerical approximations of the underlying non-local differential equations for quantities like the mean exit time and escape probability [5]. Numerous schemes have also been proposed to discretize these more complicated integro-differential equations, like those originated from random walk models [8, 9, 10] with fractional derivatives [3], those based on quadrature of the singular integral representation [4, 7, 13, 26] and harmonic extension [19].

In this paper, several finite difference schemes for the infinitesimal generator $\mathcal{A}$ will be constructed, by exploring their representations in the appropriate spectral space. For simplicity, only the normalised operator corresponding to

[TABLE]

will be considered; general cases can be adapted easily by adding discretization of classical first or second order derivatives for convection (if $\mu\neq 0$) or diffusion (if $q\neq 0$), and by applying appropriate scaling (if $\sigma\neq 1$). To motivate the general form of the schemes, basic theoretical tools related to semi-discrete Fourier transform will be reviewed in Section 2. In Section 3, specific schemes are constructed based on the approximate multipliers in the spectral space, followed by several numerical experiments in Section 4. Finally, special functions and integrals, as well as lengthy derivations of the weights in certain schemes are collected in the appendix.

2 The discrete operator and semi-discrete Fourier transform

For local operators like classical integer order derivatives, Taylor expansion is usually employed to derive difference schemes and to assess their order of accuracy. However, for nonlocal operators like fractional order derivatives [20] and infinitesimal generators like $\mathcal{A}$ above, local expansions around the point of interest soon become inadequate. The right theoretical framework has to be built upon new tools, like the semi-discrete Fourier transform for functions defined on a uniform infinite lattice.

Let $v$ be a discrete function (usually sampled from a continuous one) defined on the grid $\mathbb{Z}_{h}=\{x_{j}=hj\mid j\in\mathbb{Z}\}$ with uniform spacing $h$. Then its semi-discrete Fourier transform is defined as

[TABLE]

Here the domain is restricted to be the finite interval $I_{h}$, as $\hat{v}$ is periodic with period $2\pi/h$. The inverse transform can also be readily worked out as

[TABLE]

When the grid size $h$ goes to zero, both (8) and (9) converge to the continuous Fourier transform and inverse Fourier transform, respectively.

Though not as popular as the closely related Fourier transform and Fourier series, the transform (8) has been used for a long time in numerical analysis like in spectral methods [25] and more recently in the discretization of the generator of symmetric $\alpha$-stable Lévy processes (or better known as Fractional Laplacian) [14]. By the well-known Plancherel type theorem, if $v$ is in

[TABLE]

then its semi-discrete transform belongs to

[TABLE]

Similarly, the following Parseval’s identity holds

[TABLE]

with the usual inner products:

[TABLE]

Finally the following convolution theorem is expected [14].

Theorem 2.1** (Convolution Theorem).**

Let $v$ and $w$ be functions in $\ell^{2}(\mathbb{Z}_{h})$ together with their semi-discrete transforms $\hat{v}$ and $\hat{w}$ in $L^{2}(I_{h})$. Then the operator $\mathcal{D}:\ell^{2}(\mathbb{Z}_{h})\to\ell^{2}(\mathbb{Z}_{h})$ defined by

[TABLE]

becomes a multiplier in the spectral space $L^{2}(I_{h})$, that is,

[TABLE]

Once the foundations have been firmly laid, finite difference schemes for the infinitesimal generator $\mathcal{A}$ can be constructed based on their representations in the spectral space: since $\mathcal{A}$ corresponds to the multiplier $\psi(-\xi)$, any reasonable scheme is expected to be a multiplier $\psi_{h}(-\xi)$ in the spectral space that approximates $\psi(-\xi)$. Therefore, if the discrete operator $\mathcal{A}_{h}u$ is represented as $\widehat{\mathcal{A}_{h}u}(\xi)=\psi_{h}(-\xi)\hat{u}(\xi)$ in the spectral space, then in the physical space the scheme reads

[TABLE]

for some weights $\{w_{j}\}_{j=-\infty}^{\infty}$ with $\hat{w}(\xi)=h\psi_{h}(-\xi)$. The equivalence between these two representations can be established by the Convolution Theorem 2.1, and schemes can be constructed directly from appropriate multiplier $\psi_{h}(-\xi)$, or alternatively from weights with suitable $\psi_{h}(-\xi)$ that approximates $\psi(-\xi)$. Below we focus on the case $\alpha\neq 1$, and the case $\alpha=1$ will be treated separately.

Before discussing explicit weights detailed in the next section, we make several observations and assumptions to simplify the presentation. Motivated from \psi(\xi)=-|\xi|^{\alpha}\big{(}1-i\beta\operatorname{sgn}(\xi)\tan({\alpha\pi}/{2})\big{)}, we can assume that the most general form of $\psi_{h}$ reads

[TABLE]

where both $\tilde{M}_{eh}(\xi)$ and $\tilde{M}_{oh}(\xi)$ behave like $|\xi|^{\alpha}$ near the origin. Because of the consistency in the accuracy and the appearance of related special integrals in the expressions of the weights, it happens that in many cases $\tilde{M}_{eh}(\xi)=\tilde{M}_{oh}(\xi)$, but there is no need to impose further restrictions like this. Following the above assumption for $\psi_{h}$,

[TABLE]

and from the inverse transform (9)

[TABLE]

Next, because $\mathcal{A}$ is a spatial derivative of order $\alpha$, the only dependence of the weights $\{w_{j}\}$ on the grid size $h$ is expected to be the factor $h^{-\alpha}$. This observation implies that the rescaled symbols or multipliers $M_{e}(\xi):=h^{\alpha}\tilde{M}_{eh}(\xi/h)$ and $M_{o}(\xi):=h^{\alpha}\tilde{M}_{oh}(\xi/h)$ can be assumed to be independent on $h$. Under this final assumption, the above integral for $w_{j}$ can be further simplified by isolating the real and imaginary parts of the integrand. That is,

[TABLE]

Therefore, the weight $w_{j}$ consists of two Fourier integrals, the symmetric part $\int_{0}^{\pi}M_{e}(\xi)\cos(j\xi){\rm d}\xi$ and the anti-symmetric part $\int_{0}^{\pi}M_{o}(\xi)\sin(j\xi){\rm d}\xi$; the symmetric part corresponds to the well-known Fractional Laplacian studied thoroughly in [14] under the same framework.

Remark*.*

From the conditions $M_{e}(\xi)\sim|\xi|^{\alpha}$ and $M_{o}(\xi)\sim|\xi|^{\alpha}$, we have $\tilde{M}_{eh}(0)=\tilde{M}_{oh}(0)=0$ and

[TABLE]

As a result, the general scheme (10) can also be written as

[TABLE]

This is the form adopted in [14], for its close resemblance with the singular integral representation of the Fractional Laplacian operator.

3 Finite difference discretization of the generator $\mathcal{A}$

In this section, several finite difference schemes will be constructed with appropriate rescaled multipliers $M_{e}(\xi)$ and $M_{o}(\xi)$, such that the semi-discrete transform of the weights $\{w_{j}\}$ in the general scheme (10) becomes

[TABLE]

where the substitutions $\tilde{M}_{eh}(\xi)=h^{-\alpha}M_{e}(h\xi)$ and $\tilde{M}_{oh}(\xi)=h^{-\alpha}M_{o}(h\xi)$ are used in (11). Once $M_{e}(\xi)$ and $M_{o}(\xi)$ are chosen, the weights $\{w_{j}\}$ are then expressed explicitly as in (12), and vice versa.

Although no restrictions placed on $M_{e}(\xi)$ and $M_{o}(\xi)$ other than the behaviour $|\xi|^{\alpha}$ near the origin, there is one main constraint in practice: the Fourier integrals in (12) can not be evaluated in closed form expressions in general, while numerical quadrature could introduce large error, especially when the index $j$ of $w_{j}$ is large. In this section, we document several choices of $M_{e}(\xi)$ and $M_{o}(\xi)$, mainly motivated from schemes for the symmetric fractional Laplacian in [14]. The key features of these schemes are explicit weights expressed with special functions available in standard packages like MATLAB, making any numerical quadrature of highly oscillatory integrals dispensable.

3.1 Spectral weights

The most natural choice is $M_{e}(\xi)=M_{o}(\xi)=|\xi|^{\alpha}$, such that $\psi_{h}(\xi)$ coincides with $\psi(\xi)$ on the interval $[-\pi/h,\pi/h]$. Consequently, the resulting weights are denoted by $\{w_{j}^{SP}\}$, called spectral weights. The corresponding scheme has been studied for the symmetric fractional Laplacian ($\beta=0$): the expressions $\alpha=2$ (the classical second order derivative) appeared in the context of spectral methods [25]; general cases for $\alpha\in(0,2)$ were extensively discussed in [14] under a similar framework.

From (12), the weights $\{w_{j}^{SP}\}$ depend essentially on $\int_{0}^{\pi}\xi^{\alpha}\cos(j\xi){\rm d}\xi$ and $\int_{0}^{\pi}\xi^{\alpha}\sin(j\xi){\rm d}\xi$, which can be simplified using series expansions of the trigonometric functions. For instance,

[TABLE]

The last series can be represented using the generalised hypergeometric function ${}_{1}F_{2}$ (not the more popular Gauss hypergeometric function ${}_{2}F_{1}$, see the Appendix A for more details), leading to

[TABLE]

where $(a)_{n}=a(a+1)\cdots(a+n-1)=\Gamma(a+n)/\Gamma(a)$ is the Pochhammer symbol. Similarly,

[TABLE]

Therefore, from these explicit expressions of the two integrals,

[TABLE]

For the special case $\alpha=1$ is even simpler: $w_{0}^{SP}=-\pi/2h$ and for $j\neq 0$,

[TABLE]

where $\text{Si}(z)=\int_{0}^{z}\frac{\sin{t}}{t}{\rm d}t$ is the Sine integral.

3.2 Grünwald-Letnikov weights and fractional derivatives

Although the above scheme with the exact multiplier $M_{e}(\xi)=M_{o}(\xi)=|\xi|^{\alpha}$ is natural, the most popular one in the literature is based on classical Grünwald-Letnikov finite differences for Riemann-Liouville fractional derivatives defined on bounded intervals [3, 20]. For functions defined on the whole real line, the original Grünwald-Letnikov differences can easily be generalised by allowing integration limits in the definition going to infinity, leading to the so called Weyl fractional derivatives. Here we follow the historical development of this method and quote the weights first, before showing the multipliers. For $\alpha\in(0,1)$, the Weyl fractional derivatives $\mathcal{D}_{+}^{\alpha}$ and $\mathcal{D}_{-}^{\alpha}$ for a smooth function $f$ on the real line are defined as [22]

[TABLE]

Numerous schemes have been proposed to approximate fractional derivatives like these, including the classical Grünwald-Letnikov finite differences (see [3, 20] for more details):

[TABLE]

where $f_{j}=f(x_{j})$ and $\binom{\alpha}{k}=\frac{\Gamma(k-\alpha)}{k!\Gamma(-\alpha)}=(-1)^{k}\frac{\alpha(\alpha-1)\cdots(\alpha-k+1)}{k!}$ is the generalised binomial coefficients. The infinitesimal generator222The term $-iy\xi\mathbf{1}_{|y|\leq 1}$ is omitted in the expression of the integral, contributing to the constant $iK\xi$ in (4). $\mathcal{A}$ is then related to these fractional derivatives by

[TABLE]

where the relations (6) and (7) are used in the last step. As a result, we obtain a scheme for $\mathcal{A}$ with the fractional derivatives $\mathcal{D}_{\pm}^{\alpha}$ approximated by (14), and the weights $\{w_{j}^{GL}\}$, called Grünwald-Letnikov weights, are collected as

[TABLE]

It is easy to check that $w_{j}^{GL}$ is positive for all $\alpha\in(0,1)$ and $j\neq 0$. The associated rescaled multipliers $M_{e}(\xi)$ and $M_{o}(\xi)$ can then be obtained by using the definition of the semi-discrete transform (8) and by recognizing that the (positive) weights are proportional to the binomial coefficients of $(1-z)^{\alpha}$, i.e.,

[TABLE]

By choosing the principal branch of the power function $(1-e^{\pm i\xi})^{\alpha}$, we can further verify that $M_{e}(\xi)\sim|\xi|^{\alpha}$ and $M_{o}(\xi)\sim|\xi|^{\alpha}$ as $|\xi|\to 0$, and confirm that the scheme (10) with the weights $\{w_{k}^{GL}\}$ provides a reasonable discretization of $\mathcal{A}$.

The case when $\alpha\in(1,2)$ can be worked out similarly, starting from the fractional derivatives

[TABLE]

and their finite difference approximations

[TABLE]

where the indices are shifted by one to preserve the non-negativity of the coefficient of $f_{j\pm k}$ with $k\neq 0$. This modification is important in time-dependent problems like random walk models [9]. Hence the weights collected from (15) become

[TABLE]

Similarly, the rescaled multipliers can also be obtained as

[TABLE]

One reason for the popularity of this scheme is the non-negativity of the weights for all $\beta\in[-1,1]$, important to non-negativity of the probability densities and other stability properties when used for evolutionary equations. However, the weights becomes singular as $\alpha$ goes to $1$ (in either direction), which might be suggested by above two distinct forms of $\mathcal{D}_{\pm}^{\alpha}$; the case $\alpha=1$ has to be treated differently, for instance in [11].

3.3 Regularized Spectral weights

Finally, we consider the scheme with M_{e}(\xi)=M_{o}(\xi)=\big{(}2-2\cos(\xi)\big{)}^{\alpha/2}, which is motivated from the rescaled multiplier $2-2\cos(\xi)$ for the standard three-point central difference method for the second order derivative and can be thought as a regularised function of the exact multiplier $|\xi|^{\alpha}$ near the origin. The scheme for the symmetric case ($\beta=0$) was already studied in [27], together with different boundary conditions associated with a bounded domain. The weights for general $\beta\in(-1,1)$ are

[TABLE]

where detailed derivation is given in Appendix C.

For $\alpha=1$, we have (by taking the limit as $\alpha$ goes to $1$)

[TABLE]

but the other integral $\int_{0}^{\pi}\sqrt{2-2\cos\xi}\ln\xi\sin(k\xi){\rm d}\xi$ does not seem to possess a closed form expression and hence has to be calculated using numerical quadrature [15].

Here we only focus on schemes whose weights can be expressed explicitly as above, while in principle more general ones can be constructed from given pairs rescaled multipliers $M_{e}(\xi)$ and $M_{o}(\xi)$ that behaves like $|\xi|^{\alpha}$ near the origin. When the schemes are used in practical problems, other issues like accuracy and non-negativity of the weights start to play important roles, as commented in the numerical experiments in the next section.

4 Numerical experiments

In this section, several numerical experiments are performed to compare the proposed schemes, for their accuracy and other related issued in the discretization of practical problems.

4.1 Accuracy of the schemes

Because of the non-locality of the operator, the accuracy of different schemes is better to be assessed through the spectral space. Let $\{u_{j}\}$ be sampled from a function $u$ defined on $\mathbb{R}$. If the Fourier transform of $u$ is defined as $\hat{u}(\xi)=\int_{-\infty}^{\infty}u(x)e^{-ix\xi}\mathrm{d}x$, then

[TABLE]

Hence the error between $\mathcal{A}u(x_{j})$ and its difference approximation at $x_{j}$, provided that the infinite sum $\sum_{k}w_{j-k}u_{k}$ can be performed accurately, is given by

[TABLE]

If $\hat{u}(\xi)$ decays to zero fast enough as $|\xi|$ goes to infinity, the error is dominated by the second integral. In other words, the accuracy depends on how well $\psi(-\xi)$ is approximated by $\psi_{h}(-\xi)$, or equivalently how well $|\xi|^{\alpha}$ is approximated by $M_{e}(\xi)$ and $M_{o}(\xi)$.

For the scheme with spectral weights $\{w_{j}^{SP}\}$, $M_{e}(\xi)=M_{o}(\xi)=|\xi|^{\alpha}$ and the accuracy is expected to be spectral—the error usually decays like $O(e^{-C/h})$ for some positive constant $C$.

For the scheme with Grünwald-Letnikov weights $\{w_{j}^{GL}\}$, the rescaled multipliers can be rewritten to better illustrate their behaviours near the the origin: for $\alpha\in(0,1)$,

[TABLE]

and for $\alpha\in(1,2)$,

[TABLE]

After some complicated algebra, all these multipliers can be expanded as $|\xi|^{\alpha}(1+c_{1}\xi+c_{2}\xi^{2}+\cdots)$ for some non-zero constant $c_{1}$ (different for different multipliers). Therefore, the second integral in (18) becomes

[TABLE]

Upon substitution of the expansion for $M_{e}(\xi)$ and $M_{o}(\xi)$, the above error is bounded by $O(h)$, provided that $\int_{-\pi/h}^{\pi/h}|\xi|^{\alpha+1}|\hat{u}(\xi)|{\rm d}\xi$ is finite.

Finally for the scheme with regularised spectral weights $\{w_{j}^{RS}\}$,

[TABLE]

A similar procedure shows the leading $O(h^{2})$ error, provided that $\int_{-\pi/h}^{\pi/h}|\xi|^{\alpha+2}|\hat{u}(\xi)|{\rm d}\xi$ is finite.

Different approximations to $\mathcal{A}u$ (with $\alpha=0.5,\beta=0.5$) for the Gaussian $u(x)=\frac{1}{\sqrt{2\pi}}e^{-x^{2}/2}$ on a grid size $h=0.4$ are shown in Figure 1(a), where the ’exact’ approximation is computed using high order numerical quadrature for the inverse Fourier transform. While small error could still be observed when the Grünwald-Letnikov weights are used, the approximations are almost indistinguishable from the exact value $\mathcal{A}u$ for the spectral weights and regularised spectral weights, even on a coarse grid with $h=0.4$. The above analysis about the order of accuracy is further verified in Figure 1(b), when the $L^{\infty}$ error (computed on the interval $[-4,4]$) decreases with the expected order as the grid size $h$ is refined. The spectral convergence for $w^{SP}$ is omitted in the figure, because the error is already less than $10^{-6}$ on a grid size $h=0.4$.

Notice that the above convergence rates are only expected under certain restrictions: the function $u(x)$ should be smooth enough such that the integral $\int_{|\xi|>\pi/h}\psi(-\xi)\hat{u}(\xi)e^{ix\xi}{\rm d}\xi$ can be safely ignored; care must be taken to avoid introducing any error when the infinite sum in the convolution of the scheme (10) is truncated.

4.2 Application to mean exit time

In practice, we are more interested in solutions to non-local differential equations with the generator $\mathcal{A}$, rather than approximations of the operator alone. For example, the mean exit time appears in many systems driven by various random noises [2, 6, 17]. Let the first exit time starting at $x$ from a bounded domain $\Omega$ is defined as $\tau{(\omega,x)}:=\inf\{t\geq 0,X_{0}=x,X_{t}(\omega,x)\notin\Omega\}$, and the mean first exit time (in short, mean exit time) is $u(x)=\mathbb{E}[\tau(\omega,x)].$ If $X_{t}$ is an $\alpha$-stable Lévy motion, such that \mathbb{E}\big{(}\exp(i\xi(X_{t}-X_{0}))\big{)}=\exp(t\psi(\xi)) with

[TABLE]

then the mean exit time satisfies the following nonlocal partial differential equation [5]

[TABLE]

subject to the Dirichlet-type exterior condition $u(x)\equiv 0$ for $x\in\Omega^{c}$.

Because of the zero boundary condition outside the domain, $u_{j}\equiv 0$ for $|j|\geq 1/h$ and the unknowns $\{u_{j}\}$ are governed by a linear system of equations

[TABLE]

This system is Toeplitz and can usually be solved very efficiently [18].

The numerical solutions for $\alpha=0.5$ and $\alpha=1.5$ ($\beta=0.5$ in both cases) are show in Figure 2 with grid size $h=0.1$. As a comparison, mean exit time estimated from Monte Carlo simulation (see [16] for related random number generation) is also plotted, with 10000 sample paths for each process starting at $x$ and time step $0.0001$ for advancing the evolution. While all numerical solutions agree with the Monte Carlo simulation to some degree, new problems arise as a result of the non-smoothness of the solution $u$ near the boundaries $x=\pm 1$. For instance, the mean exit time $u(x)$ for the symmetric operator behaves like $\mbox{dist}(x,\partial\Omega)^{\alpha/2}$ near the boundary [21], which is not even Lipschitz continuous for any $\alpha\in(0,2)$. Because of this loss of regularity, the order of convergence derived in the previous subsection for each scheme no longer holds. Moreover, because of the appearance of negative weights, the numerical solutions become oscillatory, which is more pronounced when $\alpha$ is less than one, or when $|\beta|$ is close to one.

The situation is even worse for evolutionary problems when schemes with negative weights (except $w_{0}$, which is always negative yet does not affect the stability) are used—probability densities 333Strictly speaking, it is the formal adjoint operator $\mathcal{A}^{*}$ should be used for the evolution of probability densities. could become negative during the evolution. In cases the underlying process is genuinely asymmetric ($\beta\neq 0$), it is common to have negative weights, because the magnitude of the integral $\int_{0}^{\pi}M_{o}(\xi)\sin(j\xi){\rm d}\xi$ usually decays to zero slower than that of $\int_{0}^{\pi}M_{e}(\xi)\cos(j\xi){\rm d}\xi$. In this regard, only the scheme with Grünwald-Letnikov weights has the desired stability property for all $\beta\in[-1,1]$.

5 Conclusion

In this paper, we proposed several difference schemes of the general form (10) for the infinitesimal generator of the $\alpha$-stable Lévy process, mainly based on the corresponding representations under the semi-discrete Fourier transform, such that the exact multiplier is properly approximated. The scheme with spectral weights $\{w_{j}^{SP}\}$ is spectrally accurate for smooth functions, but the accuracy becomes degenerate for non-smooth functions, and the solutions show spurious oscillations in the discretization of practical problems. The scheme with classical Grünwald-Letnikov weights $\{w_{j}^{GL}\}$ exhibits good numerical stability, but it is only first order accurate and a “singularity” appears when $\alpha$ approaches one. A compromise between accuracy and stability is reached in the scheme with regularised spectral weights $\{w_{j}^{RS}\}$: the accuracy is second order, adequate for most problems, while non-physical oscillations are not as pronounced as those when the spectral weights are used. The construction of high order schemes with good stability stability remains a challenging task. In practice, if the function under consideration does not decay to zero fast enough, the accuracy of the scheme could be reduced when the discrete operator $\mathcal{A}_{h}$ is just truncated to finite sums. The next order may still be recovered in many circumstances, by combining far-field asymptotic behaviour of the function and the integral representation in (2), in a similar way as in [13, Section 5]. However, this treatment of boundary condition requires certain information like the decay rate of the function that may demand deep theoretical analysis.

Appendix A The generalised hypergeometric function

The generalised hypergeometric function appears at several places in this paper, and near the origin, it is defined as a series

[TABLE]

where $(a)_{n}=a(a+1)\cdots(a+n-1)$ is the Pochhammer symbol. The most common ones include Kummer’s confluent hypergeometric function ${}_{1}F_{1}$ and Gauss’s hypergeometric function ${}_{2}F_{1}$. From the series representation (20), it is easy to verify that

[TABLE]

Appendix B Special definite integrals

In the simplification of the function $\psi(\xi)$ when the Lévy measure is given by (3), the following special definite integrals are used (see [12, p.436-437]):

[TABLE]

Appendix C Evaluation of the integrals in the weights $w^{RS}_{j}$

Here explicit expressions of the weights $\{w_{j}^{RS}\}$ associated with the regularised multipliers M_{e}(\xi)=M_{o}(\xi)=\big{(}2-2\cos\xi)^{\alpha/2} are derived. In fact, $M_{e}(\xi)=M_{o}(\xi)=2^{\alpha}\sin^{\alpha}\left(\frac{\xi}{2}\right)$, and we can verify the following indefinite integral

[TABLE]

First, using (21) for the derivative of generalised hypergeometric functions, we get

[TABLE]

By the series representation of ${}_{2}F_{1}$, the right hand side of (24) becomes

[TABLE]

which is exactly $(1-e^{i\xi})^{\alpha}$. Therefore, the indefinite integral (23) is established by

[TABLE]

while in the last step the principal branch of the fractional power is taken with $\xi$ assumed to be in the interval $[0,\pi]$.

Extract the real and the imaginary parts of the following integral

[TABLE]

we get

[TABLE]

and

[TABLE]

where are simplified by Gauss’s identity

[TABLE]