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 |