How to do it...

  1. 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.
  2. 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.
  3. 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
  1. 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.
  2. 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)
  1. 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.

..................Content has been hidden....................

You can't read the all page of ebook, please click here login for view all page.
Reset