2017-12-18 347 views
1

这是我的代码。它是一个函数,用于评估另一个函数在某个x值下的导数。即使对于分数阶导数(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) 

我的问题是,我想整合h1这是(x-z)**(n-a-1)q1(the derivative of f(z))乘法。如果我让f(z)=z^2和手动输入2*zq1它工作正常,但如果我尝试使用q1它说“不能将表达式转换为浮动”。任何想法为什么?

+0

您还可以在代码中包含import语句吗?我猜'sp'是SymPy,'sc'是SciPy,你从哪里导入'gamma()'函数? –

+0

对不起Amit,这是我第一次使用stackoverflow。我希望澄清事情。 –

回答

0

我认为你正在将SymPy与SciPy混合在一起。 SymPy只做符号计算。 SciPy只进行数值计算。如果要计算SymPy中的某些中间结果,然后使用如果进行数值计算,则必须使用lambdifyread more here)函数,这与SymPy的lambda函数相同,该函数提供数字结果。在下面的代码中,我使用了两次lambdify函数将SymPy派生计算转换为给出SciPy预期的数值结果的lambda函数。

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