Integrating special function and plotting

Question

I am trying to plot an integration of a special (eg. Bessel) function and my minimal code is the following.

#!/usr/bin/env python
import matplotlib.pyplot as plt
import numpy as np
import scipy.integrate as integrate
import scipy.special as sp
from scipy.special import jn


#x = np.arange(0.0, 10.0, 0.1)
U = np.linspace(0,10,1000)



#Delta = U**2
#Delta = U-4+8*integrate.quad(lambda x: sp.jv(1,x)/(x*(1.0+np.exp(U*x*0.5))), 0, 100)
Delta = U-4+8*integrate.quad(lambda x: jn(1,x)/(x*(1.0+np.exp(U*x*0.5))), 0.1, 1000)

plt.plot(U,Delta)

plt.xlabel('U')
plt.ylabel('$\Delta$')
plt.show()

However, this gives me several an error messages saying quadpack.error: Supplied function does not return a valid float whereas the function gets easily plotted in Mathematica. Do Python's Bessel's functions have limitations?

I have used this documentation for my plotting.

Your integral over [0,4.5] is a single number, not an array. What are you trying to plot? — Stephen Terry, Jul 09 '17 at 14:59
@StephenTerry sorry posted a wrong question, I have edited the post. — hbaromega, Jul 09 '17 at 16:04
*"I am trying to **plot an integration** of a special (eg. Bessel) function"*. I do not know what this is... Please explain in more details what exactly are you trying to plot. — AGN Gazer, Jul 09 '17 at 16:16

AGN Gazer · Answer 1 · 2017-07-10T02:22:53.507

It is difficult to provide an answer that solves the problem before understanding what exactly you are trying to do. However, let me list a number of issues and provide an example that may not achieve what you are trying to do but at least it will provide a path forward.

Because your lambda function multiplies x by an array U, it returns an array instead of a number. A function that needs to be integrated should return a single number. You could fix this, for example, by replacing U by u:
```
f = lambda x, u: jn(1,x)/(x*(1.0+np.exp(u*x*0.5)))
```
Make Delta a function of u AND make quad pass additional argument u to f (defined in the previous point) AND extract only the value of the integral from the returned tuple from quad (quad returns a tuple of several values: the integral, the error, etc.):
```
 Delta = lambda u: -4+8*integrate.quad(f, 0.1, 1000, args=(u,))[0]
```
Compute Delta for each u:
```
 deltas = np.array(map(Delta, U))
```
plot the data:
```
 plt.plot(U, deltas)
```

Integrating special function and plotting

1 Answers1