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?
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.