Issues numerically integrating in python

Question

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?

Answer 1

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.