CODEWITHSUNDEEP
pythonscipynumerical-methods
OverLordGoldDragon
OverLordGoldDragon

Reputation: 19816

How does `points` work in `integrate.quad`?

I have an intricate function whose integration behavior I fail to make sense of:

  1. Does points specify intervals of possible singularities, or points?
  2. Whatever answer to above, how to explain 2.1 working but not 2? (also 1, 5, ...) Note that the only subject of singularity is at 0, but it's not even the case as without / w the function is 0 there and nearby.
  3. Behavior persists even with / w removed.
import numpy as np
from scipy.integrate import quad

def fn(w, mu=11.316582914572866):
    return (np.exp(-1 / (1 - ((w - mu) * (np.abs(w - mu) < .999))**2))
            ) * (np.abs(w - mu) < .999) / w

for pt2 in (1.0, 2.0, 2.1, 5.0):
    print(pt2, '--', quad(fn, 1e-8, 40, points=(0, pt2))[0])

1.0 -- 0.039282478570060606
2.0 -- 0.0
2.1 -- 0.03928247857037831
5.0 -- 0.03928247859275705

Upvotes: 1

Views: 1752

Answers (2)

Lutz Lehmann
Lutz Lehmann

Reputation: 26040

points is, per documentation, used to indicate "interesting" regions in the integration interval. The points given are used in the initial subdivision. So giving points=(mu,) will subdivide the integration on [0,40] to the intervals [0,mu] and [mu,40]. The the midpoint values will not lie on the line through the end-point values, so the error estimate will be non-zero and trigger further subdivision towards an appropriate sampling of the function. One could also give all points=(mu-1,mu,mu+1) so that the outer intervals are treated in the first step. This all needs not be very precise, it is sufficient to get subdivision points or intervals of feature size where the features of the function are.

for pt2 in (1.0, 2.0, 2.1, 5.0):
    print(pt2, '--', quad(fn, 1e-8, 40, points=(11,))[0])

1.0 -- 0.03928247857032486
2.0 -- 0.03928247857032486
2.1 -- 0.03928247857032486
5.0 -- 0.03928247857032486

or

for pt2 in (1.0, 2.0, 2.1, 5.0):
    print(pt2, '--', quad(fn, 1e-8, 40, points=(9,13))[0])

1.0 -- 0.039282478570323645
2.0 -- 0.039282478570323645
2.1 -- 0.039282478570323645
5.0 -- 0.039282478570323645

Upvotes: 2

Peter Meisrimel
Peter Meisrimel

Reputation: 402

  1. Your points should be within the integration interval you specified, i.e., 0 is not valid.
  2. quad is based on numerical integration, which consists of sampling and weighting function evaluations. Adding a print(w) inside your function allows you to see where this sampling occurs.
  3. From the output it appears that quad is run seperately on e.g. [1e-8, 2] and [2, 40].
  4. The output also shows there is no sample taken sufficiently close to mu. Thus, quad only sees a function that is almost constantly zero, hence the result.

Upvotes: 1

Related Questions