Feeds:
Posts

## PYTHON IMPLEMENTATION OF MITTAG-LEFFLER AND ITS DERIVATIVES

The Mittag-Leffler function $E_{a,b}(z) = \sum_{k=0}^\infty \frac{z^k}{\Gamma(ak+b)}$ is extremely basic in fractional calculus because $E_{a,1}(-\lambda t^a)$ occurs as the solution of the fractional differential equation $(d/dt)^a f(t) = \lambda f(t)$ and therefore is the generalization of the exponential function.  There exists a MATLAB implementation based on the paper gorenflo-mittag-leffler-computation but there was no open source python implementation.  Here is one.  The feature of this code that goes beyond the material of the wonderful paper is the following simple relation for computing the $n$-th derivative of the Mittag-Leffler function.

$E^(n)_{a,b}(z) = E^{(n-1)}_{a,b-1}(z) + (b-(n-1)a-1) E^{(n-1)}_{a,b}(z)$

which is a simple generalization of the same for $n=1$ due to Dzherbashyan.

 import numpy as np from scipy.special import gamma from scipy.integrate import quad

def Kf(r,alf,bet,z):
res = (1/(np.pi*alf))*r**((1-bet)/alf)*np.exp(-r**(1/alf))
res *= r*np.sin(np.pi*(1-bet))- z*np.sin(np.pi*(1-bet+alf))
res /= r**2 – 2*r*z*np.cos(np.pi*alf) + z**2
return res

def P(phi,alf,bet,eps,z):
w = eps**(1/alf)*np.sin(phi/alf) + phi*(1+(1-bet)/alf)
res = eps**(1+(1-bet)/alf)/(2*np.pi*alf)
res *= np.exp(eps**(1/alf))*(np.cos(w)+1j*np.sin(w))
res /= eps*np.exp(1j*phi)
return res

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 mlf(z,alf,bet,K=10):
rho = 1e-10
if alf>1:
k0 = int(np.ceil(alf)+1)
e = 0
for k in range(k0-1):
z1 = z**(1/k0)
w = np.exp(1j* 2 * np.pi*k/k0)
e += (1/k0)*mlf( z1*w, alf/k0, bet,K)
return e
if abs(z)<1e-6:
return 1/gamma(bet)
elif abs(z)np.floor(10+5*alf):
print(‘case large z’)
k0 = int(np.floor( -np.log(rho)/np.log(abs(z))))
print(k0)
print(np.angle(z))
if abs(np.angle(z)) < alf*np.pi/4 + 0.5*min(np.pi,alf*np.pi):
print(‘real case’)
e = (1/alf)*z**((1-bet)/alf)*np.exp(z**(1/alf))
for k in range(1,k0):
e += z**(-k)/gamma(bet-alf*k)
return(e)
else:
print(‘this case’)
print(‘k0=’,k0)
e = 0
for k in range(1,k0):
t = z**(-k)/gamma(bet-alf*k)
e += -t
print(‘ML=’,e)
return(e)
else:
eps = 0.1
print(‘integral case: ML=’,I1+I2)
print(‘z=’,z)
print(‘alpha=’,alf)
print(‘beta=’,bet)
return I1+I2
return 0

def mlfn( z, a, b, n):
if n==0:
return mlf(z,a,b)
else:
h = mlfn(z,a,b-1,n-1)
A = mlfn(z,a,b,n-1)
h += (b-float(n-1)*a-1)*float(A)
return (1/z)*h

def fpp2( delta, nu,n,t):
z = nu*t**delta
return z**n*mlfn(-z,delta,1,n)