Search code examples
pythonsympysymbolic-mathderivative

Sympy derivation


I'm working with EPR spectrums and using Dyson function to fit experimental data.

def dyson(x, B_0, delta_B, alpha):
  return (((delta_B + alpha*(x - B_0))/(4*(x - B_0)**2 + delta_B**2)) 
  + ((delta_B - alpha*(x + B_0))/(4*(x + B_0)**2 + delta_B**2)))

I get its derivative using Sympy:

x = Symbol('x')
delta_B = Symbol('delta_B')
B_0 = Symbol('B_0')
alpha = Symbol('alpha')
y = (((delta_B + alpha*(x - B_0))/(4*(x - B_0)**2 + delta_B**2)) 
  + ((delta_B - alpha*(x + B_0))/(4*(x + B_0)**2 + delta_B**2)))
der_dyson = lambdify((x, delta_B, B_0, alpha), y.diff(x))

However, result doesn't match calculations using finite differences method:

def der_dyson_(x, B_0, delta_B, alpha, h=1):
  return (dyson(x + h, B_0, delta_B, alpha) - dyson(x - h, B_0, delta_B, alpha))/(2*h)
fig_test, ax_test = plt.subplots()
x_test = np.linspace(-600, 600, 10000)
ax_test.plot(x_test, der_dyson_(x_test, 200, 250, 0.5, h=0.12), label='finite diff')
ax_test.plot(x_test, der_dyson(x_test, 200, 250, 0.5), label='Sympy')
ax_test.legend()

plots of derivatives

I checked out expression for derivative that Sympy gives and it seems to be correct so I can't figure out what is the problem. In the article where Dyson function defined, author adds graph of Dyson derivative which is of the form corresponding to my finite differences computation.

derivative in original article


Solution

  • Here the order of parameters is incorrect:

    der_dyson = lambdify((x, delta_B, B_0, alpha), y.diff(x))
    

    It should correspond to the order in dyson:

    def dyson(x, B_0, delta_B, alpha):
        ...
    
    def get_dyson_deriv():
        x = Symbol('x')
        delta_B = Symbol('delta_B')
        B_0 = Symbol('B_0')
        alpha = Symbol('alpha')
        y = dyson(x, B_0, delta_B, alpha)
        # NOTE order of parameters in the tuple here:
        return lambdify((x, B_0, delta_B, alpha), y.diff(x))
    
    der_dyson = get_dyson_deriv()
    

    With this change, the two graphs become nearly identical, even with h=0.12.

    Also note that you don't need to retype the expression of dyson when computing the derivative. Simply pass the Symbols directly to the dyson function.