pythonscipysympyderivativecalculus

Use of sympy and scipy results in bad behaviour, where is my error?

I am writing a code in which I calculate derivatives and integrals. For that I use sympy and scipy.integrate, respectively. But the use results in a strange erroneous, behaviour. Below is a minimal code that reproduces the behaviour:

from scipy.integrate import quadrature
import sympy
import math

def func(z):
    return math.e**(2*z**2)+3*z

z = symbols('z')
f = func(z)
dlogf_dz = sympy.Lambda(z, sympy.log(f).diff(z))
print(dlogf_dz)
print(dlogf_dz(10))

def integral(z):
    # WHY DO I HAVE TO USE [0] BELOW ?!?!
    d_ARF=dlogf_dz(z[0])
    return d_ARF

result = quadrature(integral, 0, 3)
print(result)

>>> Lambda(z, (4.0*2.71828182845905**(2*z**2)*z + 3)/(2.71828182845905**(2*z**2) + 3*z))
>>> 40.0000000000000
>>> (8.97457203290041, 0.00103422711649337)

The first 2 print statements deliver mathematically correct results, but the integration outcome result is wrong - instead of ~8.975 it should be exactly 18. I know this, because I double-checked with WolframAlpha, since dlogf_dz is mathematically quite simple, you can check for yourself here.

Strangely enough, if I just hard-code the math-expression from my first print statement into integral(z), I get the right answer:

def integral(z):
    d_ARF=(4.0*2.71828182845905**(2*z**2)*z + 3)/(2.71828182845905**(2*z**2) + 3*z)
    return d_ARF

result = quadrature(integral, 0, 3)
print(result)

>>> (18.000000063540558, 1.9408245677254854e-07)

I think the problem is that I don't provide dlogf_dz to integral(z) in a correct way, I have to somehow else define dlogf_dz as a function. Note, that I defined d_ARF=dlogf_dz(z[0]) in integral(z), else the function was giving errors.

What am I doing wrong? How can I fix the code, more precisely make dlogf_dz compatible with sympy integrations? Tnx

Solution

  • In a sympy environment where z is a symbol, I can create a numpy compatible function from a sympy expression using lambdify (not Lambda).

    In [22]: expr = E**(2*z**2)+3*z

In [23]: fn = lambdify(z, log(expr).diff(z))

In [24]: quad(fn,0,3)

Out[24]:

    enter image description here

    lambdify produced this function:

    In [25]: help(fn)
Help on function _lambdifygenerated:

_lambdifygenerated(z)
    Created with lambdify. Signature:
    
    func(z)
    
    Expression:
    
    (4*z*exp(2*z**2) + 3)/(3*z + exp(2*z**2))

    lambdify isn't fool proof but is the best tool for using sympy and scipy together.

    https://docs.sympy.org/latest/modules/utilities/lambdify.html