I'm trying to integrate a complicated combination of functions using scipy.integration.quadrature, and it's throwing accuracy warnings, and the 'Latest difference' values are (gulp) nan:
'C:\Program Files\Anaconda3\lib\site-packages\scipy\integrate\quadrature.py:199: AccuracyWarning: maxiter (1000) exceeded. Latest difference = nan'
The problem seems similar to this post, but using a different integration method, and the parameters are different.
SciPy Quad Integration: Accuracy Warning
Rather than post the original set of functions, here is an example of a simpler function that also throws the warning, although the 'Latest difference' values here are actually defined numbers.
def func(phi):
return phi**3
def func2(phi):
return 1/(phi)
def int(phi):
return func(phi)/abs(2/func2(phi)**5)
res, err = integrate.quadrature(int, 0, 1, maxiter=10)
print("The numerical result is {:f} (+-{:g})"
.format(res, err))
Question: why is this behaviour (warning and values) occurring?
I note that raising the value of the maxiter value here (eg by powers of 10) drastically changes the result, but the Latest difference values increase - suggests the integral is diverging from a result.
Interestingly, use of scipy.integration.quad with the same inputs gives the warning 'The integral is probably divergent, or slowly convergent. warnings.warn(msg, IntegrationWarning)'. So is this just a case of a bad choice of function for the integration? Note this is not the actual function, but one that seems to give similar (but not identical) behaviour.
The warning basically means that the number of iterations is not enough to find a solution within the desired tolerance. So the first step to try would be to increase the number of iterations. However, this does not work here because the function does not want to be (numerically) integrated.
What's wrong with the function? Here is an even simpler function with the same behavior:
def fint(phi):
return 1/phi
This may make the problem clearer to see: there is a pole (function returns inf) at phi=0. No matter how many iterations we use, there will always be increasingly bigger values the closer it gets to 0.
If we do not include the pole in the integration range there is no warning:
res, err = integrate.quadrature(fint, 0.5, 1, maxiter=100)
print("The numerical result is {:f} (+-{:g})".format(res, err))
The numerical result is 0.693147 (+-6.53485e-10)
(See how small the error is?)
There is one more subtlety regarding the original warning about "Latest difference = nan". This indicates that your complicated functions produce nan values. That is often a result of calculating something like 0/0, inf/inf, inf - inf, ...
We can reproduce this with another simple function:
def fint(phi):
return phi/phi
That is basically a constant value with a tiny ugliness at phi = 0.
res, err = integrate.quadrature(fint, -1, 1, maxiter=100)
print("The numerical result is {:f} (+-{:g})".format(res, err))
The numerical result is 2.000000 (+-nan)
RuntimeWarning: invalid value encountered in true_divide return phi/phi
AccuracyWarning: maxiter (100) exceeded. Latest difference = nan
To wrap it up, not every function can be integrated over any range. And not every function that can be integrated is suitable for numerical integration.