Source code for Inelastica.math.gausskronrod

from __future__ import print_function

import numpy as N


[docs] def GaussKronrod(n, tol=1e-10): r""" Computes the (:math:`2n+1`)-point Gauss-Kronrod quadrature abscissa :math:`\{x_i\}` and weights :math:`\{w_{1,i}\}`, :math:`\{w_{2,i}\}`. Kronrod adds :math:`n+1` points to an :math:`n`-point Gaussian rule. Definition: :math:`\int_{-1/2}^{1/2} f(x) dx \sim \sum_i w_{1,i} f(x_i)`. The integration range corresponds to the first Brillouin zone in Inelastica. The weights :math:`\{w_{2,i}\}` are used for error estimation. Parameters ---------- n : int Points for the Gaussian rule. tol : float The requested absolute accuracy of the abscissas. Returns ------- x : ndarray Abscissa for which the function is to be evaluated. w1 : ndarray Integration weights. w2 : ndarray Weights for error estimation. """ # Rescale to [-0.5,0.5] and expand to get all points. x, w1, w2 = _kronrod(n, tol) indx = N.argsort(x) x, w1, w2 = x[indx], w1[indx], w2[indx] if N.abs(x[0]) < 1e-7: x, w1, w2 = N.concatenate((-x, x[1:])), N.concatenate((w1, w1[1:])), N.concatenate((w2, w2[1:])) else: x, w1, w2 = N.concatenate((-x, x)), N.concatenate((w1, w1)), N.concatenate((w2, w2)) indx = N.argsort(x) x, w1, w2 = x[indx], w1[indx], w2[indx] return x/2, w1/2, w2/2
def _abwe1(n, m, tol, coef2, even, b, x): #*****************************************************************************80 ## ABWE1 calculates a Kronrod abscissa and weight. # Licensing: # This code is distributed under the GNU LGPL license. # Modified: # 30 April 2013 # Author: # Original FORTRAN77 version by Robert Piessens, Maria Branders. # Python version by John Burkardt. # Reference: # Robert Piessens, Maria Branders, # A Note on the Optimal Addition of Abscissas to Quadrature Formulas # of Gauss and Lobatto, # Mathematics of Computation, # Volume 28, Number 125, January 1974, pages 135-139. if x == 0.0: ka = 1 else: ka = 0 for iter in range(1, 51): b1 = 0.0 b2 = b[m+1-1] yy = 4.0 * x * x - 2.0 d1 = 0.0 if even: ai = m + m + 1 d2 = ai * b[m+1-1] dif = 2.0 else: ai = m + 1 d2 = 0.0 dif = 1.0 for k in range(1, m + 1): ai = ai - dif i = m - k + 1 b0 = b1 b1 = b2 d0 = d1 d1 = d2 b2 = yy * b1 - b0 + b[i-1] if not even: i = i + 1 d2 = yy * d1 - d0 + ai * b[i-1] if even: f = x * (b2 - b1) fd = d2 + d1 else: f = 0.5 * (b2 - b0) fd = 4.0 * x * d2 # # Newton correction. # delta = f / fd x = x - delta if ka == 1: break if abs(delta) <= tol: ka = 1 # # Catch non-convergence. # if ka != 1: print('') print('ABWE1 - Fatal error!') print(' Iteration limit reached.') print(' Last DELTA was %e' % (delta)) sys.exit('ABWE1 - Fatal error!') # # Computation of the weight. # d0 = 1.0 d1 = x ai = 0.0 for k in range(2, n + 1): ai = ai + 1.0 d2 = ((ai + ai + 1.0) * x * d1 - ai * d0) / (ai + 1.0) d0 = d1 d1 = d2 w = coef2 / (fd * d2) return x, w def _abwe2(n, m, tol, coef2, even, b, x): #*****************************************************************************80 ## ABWE2 calculates a Gaussian abscissa and two weights. # Licensing: # This code is distributed under the GNU LGPL license. # Modified: # 30 April 2013 # Author: # Original FORTRAN77 version by Robert Piessens, Maria Branders. # PYTHON version by John Burkardt. # Reference: # Robert Piessens, Maria Branders, # A Note on the Optimal Addition of Abscissas to Quadrature Formulas # of Gauss and Lobatto, # Mathematics of Computation, # Volume 28, Number 125, January 1974, pages 135-139. if x == 0.0: ka = 1 else: ka = 0 for iter in range(1, 51): p0 = 1.0 p1 = x pd0 = 0.0 pd1 = 1.0 if n <= 1: if x != 0.0: p2 = (3.0 * x * x - 1.0) / 2.0 pd2 = 3.0 * x else: p2 = 3.0 * x pd2 = 3.0 ai = 0.0 for k in range(2, n + 1): ai = ai + 1.0 p2 = ((ai + ai + 1.0) * x * p1 - ai * p0) / (ai + 1.0) pd2 = ((ai + ai + 1.0) * (p1 + x * pd1) - ai * pd0) / (ai + 1.0) p0 = p1 p1 = p2 pd0 = pd1 pd1 = pd2 delta = p2 / pd2 x = x - delta if ka == 1: break if abs(delta) <= tol: ka = 1 if ka != 1: print('') print('ABWE2 - Fatal error!') print(' Iteration limit reached.') print(' Last DELTA was %e' % delta) sys.exit('ABWE2 - Fatal error!') an = n w2 = 2.0 / (an * pd2 * p0) p1 = 0.0 p2 = b[m+1-1] yy = 4.0 * x * x - 2.0 for k in range(1, m + 1): i = m - k + 1 p0 = p1 p1 = p2 p2 = yy * p1 - p0 + b[i-1] if even: w1 = w2 + coef2 / (pd2 * x * (p2 - p1)) else: w1 = w2 + 2.0 * coef2 / (pd2 * (p2 - p0)) return x, w1, w2 def _kronrod(n, tol): #*****************************************************************************80 # ## KRONROD adds N+1 points to an N-point Gaussian rule. # # Discussion: # # This subroutine calculates the abscissas and weights of the 2N+1 # point Gauss Kronrod quadrature formula which is obtained from the # N point Gauss quadrature formula by the optimal addition of N+1 points. # # The optimally added points are called Kronrod abscissas. The # abscissas and weights for both the Gauss and Gauss Kronrod rules # are calculated for integration over the interval [-1,+1]. # # Since the quadrature formula is symmetric with respect to the origin, # only the nonnegative abscissas are calculated. # # Note that the code published in Mathematics of Computation # omitted the definition of the variable which is here called COEF2. # # Storage: # # Given N, let M = ( N + 1 ) / 2. # # The Gauss-Kronrod rule will include 2*N+1 points. However, by symmetry, # only N + 1 of them need to be listed. # # The arrays X, W1 and W2 contain the nonnegative abscissas in decreasing # order, and the weights of each abscissa in the Gauss-Kronrod and # Gauss rules respectively. This means that about half the entries # in W2 are zero. # # For instance, if N = 3, the output is: # # I X W1 W2 # # 1 0.960491 0.104656 0.000000 # 2 0.774597 0.268488 0.555556 # 3 0.434244 0.401397 0.000000 # 4 0.000000 0.450917 0.888889 # # and if N = 4, (notice that 0 is now a Kronrod abscissa) # the output is # # I X W1 W2 # # 1 0.976560 0.062977 0.000000 # 2 0.861136 0.170054 0.347855 # 3 0.640286 0.266798 0.000000 # 4 0.339981 0.326949 0.652145 # 5 0.000000 0.346443 0.000000 # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 30 April 2013 # # Author: # # Original FORTRAN77 version by Robert Piessens, Maria Branders. # PYTHON version by John Burkardt. # # Reference: # # Robert Piessens, Maria Branders, # A Note on the Optimal Addition of Abscissas to Quadrature Formulas # of Gauss and Lobatto, # Mathematics of Computation, # Volume 28, Number 125, January 1974, pages 135-139. # # Parameters: # # Input, integer N, the order of the Gauss rule. # # Input, real TOL, the requested absolute accuracy of the # abscissas. # # Output, real X(N+1), the abscissas. # # Output, real W1(N+1), the weights for the Gauss-Kronrod rule. # # Output, real W2(N+1), the weights for the Gauss rule. # from numpy import pi from numpy import zeros from math import floor from math import sin from math import sqrt tau_dim = (n + 1) / 2 tau_dim = int(floor(tau_dim)) b = zeros(((n+1)/2)+1) tau = zeros(tau_dim) w1 = zeros(n + 1) w2 = zeros(n + 1) x = zeros(n + 1) m = int(floor((n + 1) / 2)) even = (2 * m == n) d = 2.0 an = 0.0 for k in range(1, n + 1): an = an + 1.0 d = d * an / (an + 0.5) tau[1-1] = (an + 2.0) / (an + an + 3.0) b[m-1] = tau[1-1] - 1.0 ak = an for l in range(1, m): ak = ak + 2.0 tau[l+1-1] = ((ak - 1.0) * ak \ - an * (an + 1.0)) * (ak + 2.0) * tau[l-1] \ / (ak * ((ak + 3.0) * (ak + 2.0) \ - an * (an + 1.0))) b[m-l-1] = tau[l+1-1] for ll in range(1, l + 1): b[m-l-1] = b[m-l-1] + tau[ll-1] * b[m-l+ll-1] b[m+1-1] = 1.0 bb = sin(0.5 * pi / (an + an + 1.0)) x1 = sqrt(1.0 - bb * bb) s = 2.0 * bb * x1 c = sqrt(1.0 - s * s) coef = 1.0 - (1.0 - 1.0 / an) / (8.0 * an * an) xx = coef * x1 coef2 = 2.0 / (2 * n + 1) for i in range(1, n + 1): coef2 = coef2 * 4.0 * i / (n + i) for k in range(1, n + 1, 2): [xx, w1[k-1]] = _abwe1(n, m, tol, coef2, even, b, xx) w2[k-1] = 0.0 x[k-1] = xx y = x1 x1 = y * c - bb * s bb = y * s + bb * c if k == n: xx = 0.0 else: xx = coef * x1 [xx, w1[k+1-1], w2[k+1-1]] = _abwe2(n, m, tol, coef2, even, b, xx) x[k+1-1] = xx y = x1 x1 = y * c - bb * s bb = y * s + bb * c xx = coef * x1 if even: xx = 0.0 [xx, w1[n+1-1]] = _abwe1(n, m, tol, coef2, even, b, xx) w2[n+1-1] = 0.0 x[n+1-1] = xx return x, w1, w2