"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "mpmath/calculus/quadrature.py" between
mpmath-1.0.0.tar.gz and mpmath-1.1.0.tar.gz

About: mpmath is a Python library for arbitrary-precision floating-point arithmetic.

quadrature.py  (mpmath-1.0.0):quadrature.py  (mpmath-1.1.0)
skipping to change at line 384 skipping to change at line 384
if verbose and k % 300 == 150: if verbose and k % 300 == 150:
# Note: the number displayed is rather arbitrary. Should # Note: the number displayed is rather arbitrary. Should
# figure out how to print something that looks more like a # figure out how to print something that looks more like a
# percentage # percentage
print("Calculating nodes:", ctx.nstr(-ctx.log(diff, 10) / prec)) print("Calculating nodes:", ctx.nstr(-ctx.log(diff, 10) / prec))
ctx.prec -= extra ctx.prec -= extra
return nodes return nodes
class GaussLegendre(QuadratureRule): class GaussLegendre(QuadratureRule):
""" r"""
This class implements Gauss-Legendre quadrature, which is This class implements Gauss-Legendre quadrature, which is
exceptionally efficient for polynomials and polynomial-like (i.e. exceptionally efficient for polynomials and polynomial-like (i.e.
very smooth) integrands. very smooth) integrands.
The abscissas and weights are given by roots and values of The abscissas and weights are given by roots and values of
Legendre polynomials, which are the orthogonal polynomials Legendre polynomials, which are the orthogonal polynomials
on `[-1, 1]` with respect to the unit weight on `[-1, 1]` with respect to the unit weight
(see :func:`~mpmath.legendre`). (see :func:`~mpmath.legendre`).
In this implementation, we take the "degree" `m` of the quadrature In this implementation, we take the "degree" `m` of the quadrature
skipping to change at line 408 skipping to change at line 408
Comparison to tanh-sinh quadrature: Comparison to tanh-sinh quadrature:
* Is faster for smooth integrands once nodes have been computed * Is faster for smooth integrands once nodes have been computed
* Initial computation of nodes is usually slower * Initial computation of nodes is usually slower
* Handles endpoint singularities worse * Handles endpoint singularities worse
* Handles infinite integration intervals worse * Handles infinite integration intervals worse
""" """
def calc_nodes(self, degree, prec, verbose=False): def calc_nodes(self, degree, prec, verbose=False):
""" r"""
Calculates the abscissas and weights for Gauss-Legendre Calculates the abscissas and weights for Gauss-Legendre
quadrature of degree of given degree (actually `3 \cdot 2^m`). quadrature of degree of given degree (actually `3 \cdot 2^m`).
""" """
ctx = self.ctx ctx = self.ctx
# It is important that the epsilon is set lower than the # It is important that the epsilon is set lower than the
# "real" epsilon # "real" epsilon
epsilon = ctx.ldexp(1, -prec-8) epsilon = ctx.ldexp(1, -prec-8)
# Fairly high precision might be required for accurate # Fairly high precision might be required for accurate
# evaluation of the roots # evaluation of the roots
orig = ctx.prec orig = ctx.prec
skipping to change at line 439 skipping to change at line 439
for j in xrange(1, upto): for j in xrange(1, upto):
# Asymptotic formula for the roots # Asymptotic formula for the roots
r = ctx.mpf(math.cos(math.pi*(j-0.25)/(n+0.5))) r = ctx.mpf(math.cos(math.pi*(j-0.25)/(n+0.5)))
# Newton iteration # Newton iteration
while 1: while 1:
t1, t2 = 1, 0 t1, t2 = 1, 0
# Evaluates the Legendre polynomial using its defining # Evaluates the Legendre polynomial using its defining
# recurrence relation # recurrence relation
for j1 in xrange(1,n+1): for j1 in xrange(1,n+1):
t3, t2, t1 = t2, t1, ((2*j1-1)*r*t1 - (j1-1)*t2)/j1 t3, t2, t1 = t2, t1, ((2*j1-1)*r*t1 - (j1-1)*t2)/j1
t4 = n*(r*t1- t2)/(r**2-1) t4 = n*(r*t1-t2)/(r**2-1)
t5 = r
a = t1/t4 a = t1/t4
r = r - a r = r - a
if abs(a) < epsilon: if abs(a) < epsilon:
break break
x = r x = r
w = 2/((1-r**2)*t4**2) w = 2/((1-r**2)*t4**2)
if verbose and j % 30 == 15: if verbose and j % 30 == 15:
print("Computing nodes (%i of %i)" % (j, upto)) print("Computing nodes (%i of %i)" % (j, upto))
nodes.append((x, w)) nodes.append((x, w))
nodes.append((-x, w)) nodes.append((-x, w))
skipping to change at line 573 skipping to change at line 572
`-1 \le x \le 1`, `y \ge 0`:: `-1 \le x \le 1`, `y \ge 0`::
>>> mp.dps = 50 >>> mp.dps = 50
>>> 2*quad(lambda x: sqrt(1-x**2), [-1, 1]) >>> 2*quad(lambda x: sqrt(1-x**2), [-1, 1])
3.1415926535897932384626433832795028841971693993751 3.1415926535897932384626433832795028841971693993751
One can just as well compute 1000 digits (output truncated):: One can just as well compute 1000 digits (output truncated)::
>>> mp.dps = 1000 >>> mp.dps = 1000
>>> 2*quad(lambda x: sqrt(1-x**2), [-1, 1]) #doctest:+ELLIPSIS >>> 2*quad(lambda x: sqrt(1-x**2), [-1, 1]) #doctest:+ELLIPSIS
3.141592653589793238462643383279502884...216420198 3.141592653589793238462643383279502884...216420199
Complex integrals are supported. The following computes Complex integrals are supported. The following computes
a residue at `z = 0` by integrating counterclockwise along the a residue at `z = 0` by integrating counterclockwise along the
diamond-shaped path from `1` to `+i` to `-1` to `-i` to `1`:: diamond-shaped path from `1` to `+i` to `-1` to `-i` to `1`::
>>> mp.dps = 15 >>> mp.dps = 15
>>> chop(quad(lambda z: 1/z, [1,j,-1,-j,1])) >>> chop(quad(lambda z: 1/z, [1,j,-1,-j,1]))
(0.0 + 6.28318530717959j) (0.0 + 6.28318530717959j)
**Examples of 2D and 3D integrals** **Examples of 2D and 3D integrals**
 End of changes. 4 change blocks. 
5 lines changed or deleted 4 lines changed or added

Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)