CODEWITHSUNDEEP
pythonscipynumerical-methodsnumerical-integration
MaviPranav
MaviPranav

Reputation: 335

Different results while integrating Chebyshev weight function using different quadratures

Can someone explain the following behaviour I observe while integrating the Chebyshev weight function using 3 different routines and 2 different representations of the exponent? The expected answer is Pi in each case:

from scipy.integrate import quadrature, quad, fixed_quad

print fixed_quad(lambda x: 1/(1 - x**2)**(1/2), -1, 1)
print fixed_quad(lambda x: 1/(1 - x**2)**(0.5), -1, 1)

print quadrature(lambda x: 1/(1 - x**2)**(1/2), -1, 1) 
print quadrature(lambda x: 1/(1 - x**2)**(0.5), -1, 1)

print quad(lambda x: 1/(1 - x**2)**(1/2), -1, 1)
print quad(lambda x: 1/(1 - x**2)**(0.5), -1, 1)

which gives

(2.0000000000000009, None)
(2.8254100794589787, None)

(1.9999999999999996, 4.4408920985006262e-16)
/usr/lib/python2.7/dist-packages/scipy/integrate/quadrature.py:168: AccuracyWarning:    maxiter (50) exceeded. Latest difference = 6.965869e-04
 AccuracyWarning)
(3.107110439388189, 0.00069658693569163432)

(2.0, 2.220446049250313e-14)
(3.141592653589564, 6.200200353134733e-10)

Firstly, as may be seen, the answer depends on whether 1/2 or 0.5 is given in the exponent: why is that so?

Secondly, the result depends on the routine chosen; can someone explain why FORTRAN's QUADPACK gets the right answer but the Gaussian quadratures miss the result completely?

Note: A similar question is addressed here but does not address the two specific questions above.

Upvotes: 0

Views: 456

Answers (2)

Warren Weckesser
Warren Weckesser

Reputation: 114911

As @gjdanis points out, in python 2.7, 1/2 is 0 (unless you include from __future__ import division in your code).

Your integrand has singularities at 1 and -1. fixed_quad and quadrature perform Gaussian quadrature with a weighting function w(x) = 1, so those singularities are not handled well.

fixed_quad is not adaptive (hence the name). The default order is 5. You'll have to increase the order (by a lot) to get a reasonable approximation:

In [179]: print fixed_quad(lambda x: 1/(1 - x**2)**(0.5), -1, 1, n=100)
(3.124265558250825, None)

In [180]: print fixed_quad(lambda x: 1/(1 - x**2)**(0.5), -1, 1, n=2000)
(3.1407221810853478, None)

quadrature simply calls fixed_quad with increasing order (up to a maximum given by the maxiter argument) until the difference between successive integral estimates is sufficiently small. The warning that is printed tells you that the maximum order was reached without meeting the desired error tolerance. The default for maxiter is 50; you'll need to increase maxiter to get better results. For example, here's the result with maxiter=200:

In [2]: print quadrature(lambda x: 1/(1 - x**2)**(0.5), -1, 1, maxiter=200)
/Users/warren/local_scipy/lib/python2.7/site-packages/scipy/integrate/quadrature.py:183: AccuracyWarning: maxiter (200) exceeded. Latest difference = 4.353464e-05
  AccuracyWarning)
(3.1329074742380407, 4.3534643496823122e-05)

If you use maxiter, you should also make judicious use of miniter. quadrature naively starts at miniter and increases the order by 1 until the error estimate is small enough or maxiter is reached.

To learn more about how fixed_quad and quadrature work, take a look at the source code: https://github.com/scipy/scipy/blob/master/scipy/integrate/quadrature.py

As you point out, quad is a wrapper of the Fortran library QUADPACK. This code is far more sophisticated than the simple Gaussian quadrature used by fixed_quad and quadrature.

Upvotes: 1

rookie
rookie

Reputation: 2923

In python 2.7 (which you're using) integer division is default. This means 1/2 will evaluate to 0. If you want to use floating point division, add from __future__ import division to the top of your code.

Upvotes: 3

Related Questions