Here is my code. It is a function that evaluates the derivative of another function at some x-value. I want it to return a valid output even for fractional order derivatives(a).
from scipy.special import gamma
import scipy.integrate as integrate
import sympy as sp
import scipy as sc
import math
def f(z):
return z**2
def fracdiff(f,x,a):
if a==0:
return f(x)
else:
if math.ceil(a)-a==0:
q=sp.diff(f(z),z,a)
h=q.subs(z,x)
return h
else:
n=math.ceil(a)
g1=(1/(gamma(n-a)))
q1=sp.diff(f(z),z,n)
print(q1) # showing that q1 equals 2*z
h1= lambda z:(x-z)**(n-a-1)*2*z # for z^2 the derivative is 2*z
ans=sc.integrate.quad(h1,0,x)
r=ans[0]*g1
return r
ss=fracdiff(f,1,0.5)
My problem is that I want to integrate h1 which is the multiplication of (x-z)**(n-a-1) and q1(the derivative of f(z)). If I let f(z)=z^2 and manually input 2*z for q1 it works fine but if I try using q1 it says "can't convert expression to float". Any ideas why?
I think you are mixing up SymPy with SciPy. SymPy does only symbolic calculations. SciPy does only numerical calculations. If you want to calculate some intermediate results in SymPy and then use if for numerical calculations you will have to use lambdify (read more here) function which is like a lambda function for SymPy that gives numeric results. In the following code, I used lambdify function twice to convert your SymPy derivative calculations into lambda functions that give numeric results which SciPy expects.
import sympy as sp
import scipy as sc
import math
from scipy.special import gamma
def f(z):
return z**2
def fracdiff(f,x,a):
if a==0:
return f(x)
else:
if isinstance(a,int):
z = sp.Symbol('z')
q = sp.diff( f(z), z, a )
qf = sp.lambdify( z, q ) #Lambdify used here 1
h = qf(x)
return h
else:
n = math.ceil(a)
g1 = (1/(gamma(n-a)))
z = sp.Symbol('z')
q1 = sp.diff( f(z), z, n )
q1f = sp.lambdify( z, q1 ) #Lambdify used here 2
h1 = lambda p: q1f(p)*(x-p)**(n-a-1)
ans = sc.integrate.quad(h1,0,x)
r = ans[0]*g1
return r
ss=fracdiff(f,1,0.5) # returns 1.5045055561272689