# Issues numerically integrating in python

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)

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.