I need to integrate a function defined in 2Dims (
z and radius
r), for which I don't have an expression.
I can just query the function at any position
(z,r) and get the returned value.
I have the integration range across z:
[-z_range, z_range], which I partition into
N_z points as:
z_is = -z_range + (np.arange(N_z) + 0.5) * (2.*z_range/N_z)
For each value in
z_is, the range to integrate across
r_thresh_at_this_z is obtained from the value of
def get_r_thresh(z_i): return expression_of_z_i # returns a positive float
So the range on the radial integral is dependent on the value on
I have the function
def f(r,z): return interpolator_for_f(r,z)
I want to use the
quadpy package in the most efficient way as it has been created to be able to be used in this way.
I was thinking to use a for loop to loop through the
z_is and to perform a gaussian quadrature across
[0, r_thres_for_that_z], for each value in
I could use:
results = np.zeros((N_z)) errs = np.zeros((N_z)) for i in range(N_z): def f(r): return interpolator_for_f(r,z_is[i]) results[i], errs[i] = quadpy.quad(f, [0, r_thresh_at_this_z])
But I feel a for loop is not the most efficient way to use quadpy.
Can you tell me what I am missing in doing this integral only with fast numpy arrays, so no for-loops?
[I have read about the shapes of the input x to the function f. In my case d = 1 because it's a line integral across r. n = N_z because I want to perform N_z such line integrals which I will then add up to obtain 1 single scalar, the result of the (whole) double integral.
p = 1000 because say I want 1000 integration points across r, for each value of z.
So I will need to sample the function at N_z * 1000 points.
Function f shall return an array shaped (N_z, 1000)
Is the identification of these parameters helpful in any way?]