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:
numpy.trapz function in your code calculates the integral over
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
It is like you define
F(z)=z and sum
z in [0, 1]. It is going to explode if you increase the the number of
z you chose.