- Weighted functions can be realized as products of the f(x)w(x) kind for some smooth function f(x) with a non-negative weight function w(x) containing singularities.
- An illustrative example is given by cos(πx/2)/√x. We could regard this case as the product of cos(πx/2) with w(x)=1/√x. The weight presents a single singularity of x=0, and is smooth otherwise.
- The usual way to treat these integrals is by means of weighted Gaussian quadrature formulas. For example, to perform principal value integrals of functions of the form f(x)/(x-c), we issue quad with the weight='cauchy' and wvar=c arguments. This calls the routine QAWC from QUADPACK.
Let's experiment with the Fresnel-type sine integral of g(x) = sin(x)/x on the [-1,1] interval and compare it with romberg:
In [56]: value, abs_error = quad(f, -1, 1, weight='cauchy',wvar=0);
....: print (value)
1.89216614073
In [57]: romberg(g, -1, 1)
Out[57]: 2.35040238729
- In the case of integrals of functions with weights (x-a)α(b-x)β, where a and b are the endpoints of the domain of integration and both alpha and beta are greater than -1, we issue quad with the weight='alg' and wvar=(alpha, beta) arguments. This calls the QAWS routine from QUADPACK.
- Let's experiment with the Fresnel-type cosine integral of g(x)=cos(πx/2)/√x on the [0,1] interval and compare it with quadrature:
In [58]: def f(x): return np.cos(np.pi * x * 0.5)
In [59]: def g(x): return np.cos(np.pi * x * 0.5) / np.sqrt(x)
In [60]: value, abs_error = quad(f, 0, 1, weight='alg',
wvar=(-0.5,0));
....: print (value)
1.55978680075
In [61]: quadrature(g, 0, 1)
quadrature.py:178: AccuracyWarning: maxiter (50) exceeded. Latest
difference = 3.483182e-04
Out[61]: (1.5425452942607543, 0.00034831815190772275)
- If the weight has the form w(x)=(x-a)α(b-x)βln(x-a), w(x)=(x-a)α(b-x)βln(b-x), or w(x)=(x-a)α(b-x)βln(x-a)ln(b-x), we issue quad with the weight='alg-loga', weight='alg-logb', or weight='alg-log' argument respectively, and in each case, wvar=(alpha, beta). For example, for the function g(x)=x2ln(x) on the interval [0,1], we can issue the following:
In [62]: def f(x): return x**2
In [63]: def g(x): return x**2 * np.log(x)
In [64]: value, abs_error = quad(f, 0, 1, weight='alg-loga',
wvar=(0,0));
....: print (value)
-0.111111111111
The actual value of this integral is -1/9.