I have an integral of a Bessel function that I want to calculate numerically.

I want to find the numeric integral of a function of z and w, across w, so that my result is an array across z. That is, I want

F(z) = integral f(z,w)dw

I thought I managed to accomplish this efficiently with numpy using the following code:

```
def readinKernel(wdummy, z, Ec, Ep, kval=1, ic = 1, od = 10000):
return (Ec * kval * special.jv(0, 2* Ec * kval * np.sqrt(np.outer(z, (1 - wdummy))))
* np.sqrt(ic)*Ep )
def steep_sigmoid(x, k=50):
return 1.0 / (1.0 + np.exp(-k*x))
def spinwave_calculation(z_values, Ec, Ep, numpoints = 100):
wdummy_values = np.linspace(1e-10, 1-1e-10, numpoints)
readin_values = readinKernel(wdummy_values, z_values, Ec, Ep)
readin_integrals = np.trapz(readin_values, wdummy_values, axis=1)
return readin_integrals
numpoints = 1000
z_values = np.linspace(1e-10, 1-1e-10, numpoints)
E_c_val = 1
E_p_val = 1
outputspinwave = spinwave_calculation(z_values, E_c_val, E_p_val, numpoints)
output = np.sum(np.square(outputspinwave))
print(output)
```

But it appears as though this solution appears to "blow up" with the number of points that I use. So if I increase numpoints by 10, then my output integral is 10 times larger.

How do I get an accurate estimate of this integral, and is it possible to keep both accuracy and speed?

### >Solution :

This feels more like a math problem:

The `numpy.trapz`

function in your code calculates the integral over `w`

. Your `numpoints`

is the number of `z`

points to get the array of `F(z)`

. The sum of **squared** `F(z)`

would definitely increase when you increase the number of `z`

points.

It is like you define `F(z)=z`

and sum `F(z)^2`

over `z`

in [0, 1]. It is going to explode if you increase the the number of `z`

you chose.