Problems with quad when using lambdify

Question

I'm trying to solve these two integrals, I want to use a numerical approach because C_i will eventually become more complicated and I want to use it for all cases. Currently, C_i is just a constant so _quad is not able to solve it. I'm assuming because it is a Heaviside function and it is having trouble finding the a,b. Please correct me if I'm approaching this wrongly.

Equation 33

In [1]: import numpy as np
   ...: import scipy as sp
   ...: import sympy as smp
   ...: from sympy import DiracDelta
   ...: from sympy import Heaviside

In [2]: C_i = smp.Function('C_i')

In [3]: t, t0, x, v = smp.symbols('t, t0, x, v', positive=True)

In [4]: tot_l = 10

In [5]: C_fm = (1/tot_l)*v*smp.Integral(C_i(t0), (t0, (-x/v)+t, t))

In [6]: C_fm.doit()
Out[6]: 
0.1*v*Integral(C_i(t0), (t0, t - x/v, t))

In [7]: C_fm.doit().simplify()
Out[7]: 
0.1*v*Integral(C_i(t0), (t0, t - x/v, t))

In [8]: C_fms = C_fm.doit().simplify()

In [9]: t_arr = np.arange(0,1000,1)

In [10]: f_mean = smp.lambdify((x, v, t), C_fms, ['scipy', {'C_i': lambda e: 0.8}])

In [11]: try2 = f_mean(10, 0.1, t_arr)
Traceback (most recent call last):

  File "/var/folders/rd/wzfh_5h110l121rmlxn61v440000gn/T/ipykernel_3164/3786931540.py", line 1, in <module>
    try2 = f_mean(10, 0.1, t_arr)

  File "<lambdifygenerated-1>", line 2, in _lambdifygenerated
    return 0.1*v*quad(lambda t0: C_i(t0), t - x/v, t)[0]

  File "/opt/anaconda3/lib/python3.9/site-packages/scipy/integrate/quadpack.py", line 348, in quad
    flip, a, b = b < a, min(a, b), max(a, b)

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

Equation 34

In [12]: C_i = smp.Function('C_i')

In [13]: t, tao, x, v = smp.symbols('t, tao, x, v', positive=True)

In [14]: I2 = v*smp.Integral((C_i(t-tao))**2, (tao, 0, t))

In [15]: I2.doit()
Out[15]: 
v*Integral(C_i(t - tao)**2, (tao, 0, t))

In [16]: I2.doit().simplify()
Out[16]: 
v*Integral(C_i(t - tao)**2, (tao, 0, t))

In [17]: I2_s = I2.doit().simplify()

In [18]: tao_arr = np.arange(0,1000,1)

In [19]: I2_sf = smp.lambdify((v, tao), I2_s, ['scipy', {'C_i': lambda e: 0.8}])

In [20]: try2 = I2_sf(0.1, tao_arr)
Traceback (most recent call last):

  File "/var/folders/rd/wzfh_5h110l121rmlxn61v440000gn/T/ipykernel_3164/4262383171.py", line 1, in <module>
    try2 = I2_sf(0.1, tao_arr)

  File "<lambdifygenerated-2>", line 2, in _lambdifygenerated
    return v*quad(lambda tao: C_i(t - tao)**2, 0, t)[0]

  File "/opt/anaconda3/lib/python3.9/site-packages/scipy/integrate/quadpack.py", line 351, in quad
    retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,

  File "/opt/anaconda3/lib/python3.9/site-packages/scipy/integrate/quadpack.py", line 463, in _quad
    return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)

  File "/opt/anaconda3/lib/python3.9/site-packages/sympy/core/expr.py", line 345, in __float__
    raise TypeError("Cannot convert expression to float")

TypeError: Cannot convert expression to float

Answer 1

So you are passing an unevaluated Integrate to lambdify, which in turn translates it call to scipy.integrate.quad.

Looks like the integrals can't be evaluated even with doit and simplify calls. Have you actually looked at C_fms and I2_s? That's one of the first things I'd do when running this code!

I've never looked at this approach. I have seen people lambdify the objective expression, and then try to use that in quad directly.

quad has specific requirements (check the docs!). The objective function must return a single number, and the bounds must also be numbers.

In the first error, you are passing array t_arr as the t bound, and it got the usual ambiguity error when checking where it is bigger than the other bound, 0. That's that b < a test. quad cannot use arrays as bounds.

I not sure why the second case gets avoids this problem - bounds must be coming from somewhere else. But the error comes when quad calls the objective function, and expects a float return. Instead the function returns a sympy expression which sympy can't convert to float. My guess there's some variable in the expression that's still a sympy.symbol.

In diagnosing lambdify problems, it's a good idea to look at the generated code. One way is with help on the function, help(I2_sf). But with that you need to be able to read and understand python, including any numpy and scipy functions. That's not always easy.

Have you tried to use sympy's own numeric integrator? Trying to combine sympy and numpy/scipy often has problems.