0

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
OverLordGoldDragon
  • 1
  • 9
  • 53
  • 101

2 Answers2

1
  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.
Peter Meisrimel
  • 392
  • 2
  • 10
  • Changing 0 to 1e-8 doesn't make a difference (at least here), but I did confirm 4 via 2, which helps. The `points` behavior still remains a mystery, though, and the "intervals vs points" question wasn't addressed. – OverLordGoldDragon Nov 24 '20 at 21:43
  • I do not have a guaranteed answer for the interval vs. points question, but I am leaning strongly towards intervals. E.g., using `points = [2]`, yields the intervals `[1e-8, 2]` and `[2, 40]`. The sampled points suggest `quad` is run successively on those intervals. Checking the source code did not provide me with extra cues on what exactly is going on. – Peter Meisrimel Nov 25 '20 at 09:06
  • +1 for useful info, and I agree `points` is likelier intervals, but keeping q open for a bit more definitive take. Seems to me that `points` shifts the sampling space, with a chance to entirely miss the narrow non-zero region. I ultimately went with `trapz` since key behavior is easy to predict. – OverLordGoldDragon Nov 25 '20 at 13:41
1

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
Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
  • Thanks; so `points=(10, 30)` would divide into `[0, 10]`, `[10, 30]`, and `[30, 40]`? Is each interval sampled "equally rigorously" (e.g. `linspace(0, 10, 1000)` and `linspace(10, 30, 1000)`)? And `pt2` in your demos is unfunctional, btw. – OverLordGoldDragon Dec 01 '20 at 19:53
  • Yes to the subdivision, No, to the conclusion, the error is computed by a simple halving subdivision. So if the endpoints have values almost zero and the midpoint too, the integral over the interval is taken as zero. You could of course give `points=arange(0.1,40,0.5)` or so to force a dense-enough sampling. – Lutz Lehmann Dec 01 '20 at 20:20