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

September 1, 2015 by zulfahmed

The fractional differential equation whose solution is simply when is given by the Mittag-Leffler function when where

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 where has symbol which is its Fourier multiplier. With and repesenting Laplace and Fourier transforms,

and the Green’s function is

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, . 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)

### Like this:

Like Loading...

*Related*

## Leave a Reply