Feeds:
Posts

## HURST EXPONENTS OF THE SIMPLEST TIME-FRACTIONAL DIFFERENTIAL EQUATION

The fractional differential equation $(d/dt)^\gamma y + y = 0$ whose solution is simply $y(t) = e^{-t}$ when $\gamma=1$ is given by the Mittag-Leffler function $E_{\gamma,1}(-t^\gamma)$ when $\gamma \ne 1$ where

$E_{\alpha,\beta}( z ) = \sum_{j=0}^\infty \frac{z^j}{\Gamma(\alpha j + \beta)}$

Efficient open source code for this exists in MATLAB but it is annoying to translate and debug into python.  The straightforward non-efficient approximation by direct evaluation of the series is quick and allows us to take a look at the Hurst exponent of the solution vector and the series of increments.  Our interest is to look at the Hurst exponents of fractional partial differential equations, i.e. $((d/dt)^\gamma + (d/dx)^\delta)u=0$ whose fundamental solutions are generalizations of the Gaussians for brownian motion.  But to ensure that we have working code, let’s do the simpler exercise first:

import numpy as np
from scipy.special import gamma
import matplotlib.pyplot as plt

def mlf2(z,alf,bet,K=10):
e = 1/gamma(bet)
for k in range(1,K):
e += z**k/gamma(alf*k+bet)
return e

def genFD():
t = np.array(range(1,501))
X = np.zeros((500,1)).flatten()
alf = 0.6
bet = 1.0
for j in range(500):
X[j] = mlf2( -t[j]**alf, alf,bet)
return X

>> x=genFD()
>>> res=rescaled_range(np.diff(x))
>>> res[1]
0.88576770698609264

So the Hurst exponent is 0.88 for the difference of the ordinary fractionals diffusion.  Now let’s do the exercise of taking the fundamental solution of the space-time equation as probability densities, sampling from these, and then calculating the Hurst exponent of the resulting sequence.

Consider how now the fractional diffusion equation $(D^\gamma_t + A)u(x,t)$ where $A$ has symbol $\psi$ which is its Fourier multiplier.  With $L$ and $F$ repesenting Laplace and Fourier transforms,

$s^\gamma LF u(s,k) = s^{\gamma-1} Fu(0^+,k) - \psi(k) LFu(s,k)$

and the Green’s function is

$G(t,x) = F^{-1} ( E_{\gamma,\rho}( -(1+\psi(k))t^\gamma)/E_{\gamma,\rho}(-t^\gamma))$

So we need to use this as probability density from which to sample, and then calculate the Hurst exponent of the resulting sample.  For Riesz-Feller derivative, $\psi(k) = |k|^\alpha e^{i\theta\pi/2}$.  More details can be found in Mainardi et. al. 2007 paper green-function-fracdiff-mainardi.

def green(t,gm,alf):
N = 100
k = np.zeros((N,1)).flatten()
E = np.zeros((N,1)).flatten()
et = mlf(-t**gm)
for l in range(N):
k[l] = 2*np.pi*l/N
psi = abs(k[l])**alf
ek = mlf2( -(1+psi)*t**alf)
E[k] = ek/et

return np.fft.ifft(E)