Scipy's quad fails on "pretty smooth" function

Question

The data generating functions are quite messy, so I'll try to be as clear as I can be. Please comment if not good enough at doing so.

I'm having a function in several variables and try to integrate it, holding the other one constant. However, for different values of that secondary variable, the integration process yields completely different outcomes - despite the function not altering much!

Here's my sample code (unfort. not reproducible):

fig, ax = plt.subplots()
tPrimeGrid = [0.16, 0.18, 0.22]
from scipy import integrate
for ttPrime in tPrimeGrid:
    # compute grid
    lowerBar = continuousTime.getPLowerBarAnalytical(ttPrime, Param.S, Param)
    upperBar = Param.r

    # integrate function
    h = integrate.quad(myFunc, lowerBar, upperBar, args=(0.6, ttPrime, Param), full_output=True)
    print('tPrime {} value {} erorr {}'.format(ttPrime, h[0], h[1]))
    # plot function to be evaluated
    pGrid = np.linspace(lowerBar, upperBar, 100)
    plt.plot(pGrid, [myFunc(pp, 0.6, ttPrime, Param) for pp in pGrid], label='t {} (h: {})'.format(ttPrime, h[0]))
plt.legend()

The output of quad is

tPrime 0.16 value 6.63310536371 erorr 0.000345564616621
tPrime 0.18 value 5.93645658492 erorr 0.00045956820487
tPrime 0.22 value 0.359208009237 erorr 3.98801002485e-15

Errors are all well within my level of tolerance, but the value jump a lot. Here's the graph:

So these functions have very similar shape. Looking at it with the naked eye, the difference in integrals doesn't make sense indeed. So then I computed them "manually":

(np.array([myFunc(pp, 0.6, ttPrime, Param) for pp in pGrid])).sum()*(pGrid[1]-pGrid[0])
0.35538925924926557

Implying that the small value actually is the correct one. Therefore, I'll also add here the full output of quad() for one of the bad ones:

integrate.quad(myFunc, lowerBar, upperBar, args=(0.6, 0.16, Param), full_output=True)
Out[19]: 
(6.634157704675579,
 0.004721148834250418,
 {'alist': array([ 1.99994895,  1.78826738,  1.86060326,  1.79090489,  1.93030163,
          1.72120652,  1.96515082,  1.75605571,  1.98257541,  1.7734803 ,
          1.9912877 ,  1.7821926 ,  1.99564385,  1.78872682,  1.99782193,
          1.78654874,  1.99940018,  1.78763778,  1.99945548,  1.99891096,
          1.78845456,  1.99972774,  1.99918322,  1.78831843,  1.99988939,
          1.99931935,  1.7881823 ,  1.99999149,  1.99942145,  1.7882844 ,
          1.9998979 ,  1.9999447 ,  1.99940443,  1.78825036,  1.99996597,
          1.99986387,  1.99991492,  1.99995746,  1.99938742,  1.78827589,
          1.99993194,  1.99990641,  1.99998298,  1.99988089,  1.99995321,
          1.99939592,  1.78827164,  1.99994044,  1.99995108,  1.99940231]),
  'blist': array([ 1.99995108,  1.78827164,  1.93030163,  1.86060326,  1.96515082,
          1.75605571,  1.98257541,  1.7734803 ,  1.9912877 ,  1.7821926 ,
          1.99564385,  1.78654874,  1.99782193,  1.79090489,  1.99891096,
          1.78763778,  1.99940231,  1.7881823 ,  1.99972774,  1.99918322,
          1.78872682,  1.99986387,  1.99931935,  1.78845456,  1.9998979 ,
          1.99938742,  1.78825036,  2.        ,  1.99945548,  1.78831843,
          1.99990641,  1.99994895,  1.99942145,  1.78826738,  1.99998298,
          1.99988089,  1.99993194,  1.99996597,  1.99939592,  1.7882844 ,
          1.99994044,  1.99991492,  1.99999149,  1.99988939,  1.99995746,
          1.99940018,  1.78827589,  1.9999447 ,  1.99995321,  1.99940443]),
  'elist': array([  4.61134496e-03,   4.14135579e-06,   1.27804393e-15,
           1.55808958e-15,   5.54429458e-16,   0.00000000e+00,
           2.90617202e-16,   0.00000000e+00,   2.27720143e-15,
           0.00000000e+00,   3.91514838e-15,   0.00000000e+00,
           2.15987477e-14,   5.40749179e-17,   1.19682144e-12,
           0.00000000e+00,   1.03521880e-04,   0.00000000e+00,
           2.48275100e-08,   6.85735811e-13,   6.78454062e-18,
           8.65204618e-09,   1.64268148e-13,   3.39437556e-18,
           4.65487056e-08,   1.12644353e-13,   0.00000000e+00,
           6.19443638e-08,   4.08066769e-09,   8.48813317e-19,
           5.99147851e-07,   1.21478039e-06,   9.39741081e-10,
           0.00000000e+00,   6.32087778e-14,   2.19912539e-10,
           6.97448944e-09,   1.85114515e-13,   1.28558440e-13,
           2.12217047e-19,   8.89451362e-08,   8.45167083e-09,
           1.25621961e-14,   2.59393721e-08,   2.45738052e-15,
           1.19106919e-12,   1.06110581e-19,   4.91812027e-08,
           2.72831861e-15,   3.06819516e-12]),
  'iord': array([ 1, 17, 50, 49, 48, 47, 46, 45, 44, 43, 42, 29, 40, 39, 38, 23, 35,
         11, 34,  3,  7, 21, 30, 18, 12,  8,  6,  0,  0,  0,  0,  0,  0,  0,
          0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0], dtype=int32),
  'last': 50,
  'neval': 2079,
  'rlist': array([  1.04205921e-01,   6.52426361e-06,   1.15115963e-01,
           1.40340233e-01,   4.99385660e-02,   0.00000000e+00,
           2.33230318e-02,   0.00000000e+00,   1.12779305e-02,
           0.00000000e+00,   5.54633015e-03,   0.00000000e+00,
           2.75040018e-03,   4.87063561e-03,   1.36955727e-03,
           0.00000000e+00,   8.85474806e-03,   0.00000000e+00,
           1.73731989e+00,   3.41803845e-04,   6.11097092e-04,
           1.73678061e+00,   1.70814248e-04,   3.05738170e-04,
           2.00530044e-01,   8.53852175e-05,   0.00000000e+00,
           5.29681716e-06,   1.51991445e-01,   7.64543067e-05,
           2.17985628e-01,   2.00512965e-01,   7.26773425e-02,
           0.00000000e+00,   1.06564581e-05,   3.34545761e-01,
           5.59012296e-01,   5.32838374e-06,   1.06721256e-05,
           1.91148122e-05,   3.34512038e-01,   2.38773007e-01,
           5.32807434e-06,   1.85664300e-01,   2.66423055e-06,
           5.33597722e-06,   9.55759147e-06,   1.85647516e-01,
           1.33212494e-06,   8.93844378e-03])},
 'The maximum number of subdivisions (50) has been achieved.\n  If increasing the limit yields no improvement it is advised to analyze \n  the integrand in order to determine the difficulties.  If the position of a \n  local difficulty can be determined (singularity, discontinuity) one will \n  probably gain from splitting up the interval and calling the integrator \n  on the subranges.  Perhaps a special-purpose integrator should be used.')

So here are my related questions:

1. How can it be that with such similar shape, some parameters work and some don't? How can I "predict" these failures in the future?
2. The full output (below) tells me that maximum number of subdivisions has been achieved. It does not tell me that the error is potentially miss-estimated. That must be implied? 5.9-0.0004 is not close to 0.3.
3. As increasing limit did not help (see below), what are potential alternatives? How can I integrate this function?

I tried increasing limit, but that only gave me the following (not better) output:

integrate.quad(expectedUtility, lowerBar, upperBar, args=(0.6, 0.16, Param), full_output=True, limit=500)
Out[24]: 
(6.633112814769514,
 4.74743687572826e-06,
[...]
 'The occurrence of roundoff error is detected, which prevents \n  the requested tolerance from being achieved.  The error may be \n  underestimated.')

Answer 1

quad has a somewhat funky convention for the return value when full_output is True. If quad doesn't detect any problems, it returns y, abserr, infodict. If quad detects a failure of some form, it returns y, abserr, infodict, message. infodict doesn't contain a simple status field to indicate a failure, so you have to check for the existence of the fourth return value. If it is there, it is a string that describes the problem. (In your code, you might add something like if len(h) > 3: ...handle the failure....) In your full output of a bad case, you can see that message is the string containing this:

The maximum number of subdivisions (50) has been achieved.
If increasing the limit yields no improvement it is advised to analyze
the integrand in order to determine the difficulties.  If the position of a
local difficulty can be determined (singularity, discontinuity) one will
probably gain from splitting up the interval and calling the integrator
on the subranges.  Perhaps a special-purpose integrator should be used.

This means the numerical integration has failed. Without knowing more about the integrand myFunc, it is hard to say why it failed.