quadrature

Function summary

gauss_kronrod_abscissae_and_weights(n) Gauss-Kronrod quadrature abscissae and weights
gauss_legendre_abscissae_and_weights(n) Gauss-Legendre quadrature abscissae and weights
gk_quad(f, a, b[, args, n, sum_intervals]) Integration by Gauss-Kronrod quadrature between intervals
gl_quad(f, a, b[, args, n, shanks_ind, …]) Integration by Gauss-Legendre quadrature with subdivided interval
shanks(seq[, ind]) Iterated Shanks transformation to accelerate series convergence
shanks_table(seq[, table, randomized]) Copied from sympy.mpmath.mpmath.calculus.extrapolation.py

Module listing

Numerical integration by quadrature

geotecha.mathematics.quadrature.gauss_kronrod_abscissae_and_weights(n)[source]

Gauss-Kronrod quadrature abscissae and weights

Coarse integral = Sum(f(xi) * wi1)

Fine integral = Sum(f(xi) * wi2)

For the coarse integral the unused weights are set to zero

Parameters:
n : [2-20, 32, 64, 100]

number of integration points for the Gauss points. Number of Kronrod points will automatically be 2 * n + 1.
Returns:
xi : 1d array

Abscissae for the quadrature points.

wi1 : 1d array

Weights for the coarse integral.

wi2 : 1d array

Weights for the fine integral

References

[2]Holoborodko, Pavel. 2011. ‘Gauss-Kronrod Quadrature Nodes and Weights. November 7. http://www.advanpix.com/2011/11/07/gauss-kronrod-quadrature-nodes-weights/#Tabulated_Gauss-Kronrod_weights_and_abscissae
geotecha.mathematics.quadrature.gauss_legendre_abscissae_and_weights(n)[source]

Gauss-Legendre quadrature abscissae and weights

Integral = Sum(f(xi) * wi)

Parameters:
n : [2-20, 32, 64, 100]

Number of integration points
Returns:
xi, wi : 1d array of len(n)

Abscissae and weights for numericla integration

References

[1]Holoborodko, Pavel. 2014. ‘Numerical Integration’. Accessed April 24. http://www.holoborodko.com/pavel/numerical-methods/numerical-integration/.
geotecha.mathematics.quadrature.gk_quad(f, a, b, args=(), n=10, sum_intervals=False)[source]

Integration by Gauss-Kronrod quadrature between intervals

Parameters:
f : function or method

Function to integrate.

a, b : 1d array

Limits of integration. Must have len(a)==len(b).

args : tuple, optional

args will be passed to f using f(x, *args). Default args=().

n : [7,10,15,20,25,30], optional

Number of gauss quadrature evaluation points. Default n=10. There will be 2*n+1 Kronrod quadrature points. sum_intervals : [False, True]

If sum_intervals=True the integral for each a and b, will be summed.

Otherwise each interval integration will be returned. The sum of the error estimates will also be summed.
Returns:
igral : ndarray

Integral of f between a and b. If sum_intervals=False then shape of igral will be (len(a), …) where … corresponds to however many dimensions are returned from f with scalar arguments. Each value in igral corresponds to the corresponding a-b interval. If sum_intervals=True then igral will have shape (…).

err_estimate : ndarray same size as igal

Estimate of the error in the integral. i.e. absolute value of fine integral minus coarse integral.
geotecha.mathematics.quadrature.gl_quad(f, a, b, args=(), n=10, shanks_ind=False, sum_intervals=False)[source]

Integration by Gauss-Legendre quadrature with subdivided interval

Parameters:
f : function or method

function to integrate. Must accept vector aguments for x. Might need to use numpy.vecotrize.

a, b : 1d array

limits of integration

args : tuple, optional

args will be passed to f using f(x, *args). default=()

n : [2-20, 32, 64, 100], optional

number of quadrature evaluation points. default=10

sum_intervals : [False, True]

If sum_intervals=True the integral for each a and b, will be summed. Otherwise each interval integration will be returned.
Returns:
igral : ndarray

Integral of f between a and b. If sum_intervals=False then shape of igral will be (len(a), …) where … corresponds to however many dimensions are returned from f with scalar arguments. Each value in igral corresponds to the corresponding a-b interval. If sum_intervals=True then igral will have shape (…).

Notes

Be careful when using large values of n.There may be precision issues. If f returns an ndarray when x is scalar. igral will have additonal dimensions corresponding to those of the f-with-scalar-x output.

geotecha.mathematics.quadrature.shanks(seq, ind=0)[source]

Iterated Shanks transformation to accelerate series convergence

Though normally applied to a 1d array, shanks will actually operate on the last dimension of seq which allows for multi-dimensional arrays. e.g. for 2d data each row of sequence whould be a separate sequence

Parameters:
seq : list or array

If seq is a numpy array then it’s elements will be modified in-place. If seq is a list then seq will not be modified.

ind : int, optional

Start index for extrapolation. Can be negative, e.g. ind=-5 will extrapolate based on the last 5 elements of the seq. default ind=0 i.e. use all elements.
Returns:
out : array with 1 dim less than seq, or float if seq is only 1d.

Extrapolated value. If seq is a numpy array then due to in-place modification the result will also be in seq[…, -1].

See also

shanks_table
Copy of sympy.mpmath.calculus.extrapolation.shanks Provides the whole epsilon table and error estimates.
numpy.apply_along_axis
If your sequence is not in the last dimension of an array then use np.apply_along_axis to apply it along a specific axis.

Notes

I think this will also work on multi-dimensional data. The shanks extrapolation will be performed on the last dimension of the data. So for 2d data each row is a separate sequence.

For sequence:

The partial sum is first defined as:

\[A_n=\sum_{m=0}^{n}a_m\]

This forms a new sequence, the convergence of which can be sped up by repeated use of:

\[S(A_n)=\frac{A_{n+1}A_{n-1}-A_n^2}{A_{n+1}-2A_n+A_{n-1}}\]
geotecha.mathematics.quadrature.shanks_table(seq, table=None, randomized=False)[source]

Copied from sympy.mpmath.mpmath.calculus.extrapolation.py

This shanks function is taken almost verbatim (minus an initial ctx argument???) from sympy.mpmath.mpmath.calculus.extrapolation.py:

mpmath is BSD license

Notes

Given a list seq of the first N elements of a slowly convergent infinite sequence (A_k), shanks() computes the iterated Shanks transformation S(A), S(S(A)), ldots, S^{N/2}(A). The Shanks transformation often provides strong convergence acceleration, especially if the sequence is oscillating.

The iterated Shanks transformation is computed using the Wynn epsilon algorithm (see [1]). shanks() returns the full epsilon table generated by Wynn’s algorithm, which can be read off as follows:

  • The table is a list of lists forming a lower triangular matrix, where higher row and column indices correspond to more accurate values.
  • The columns with even index hold dummy entries (required for the computation) and the columns with odd index hold the actual extrapolates.
  • The last element in the last row is typically the most accurate estimate of the limit.
  • The difference to the third last element in the last row provides an estimate of the approximation error.
  • The magnitude of the second last element provides an estimate of the numerical accuracy lost to cancellation.

For convenience, so the extrapolation is stopped at an odd index so that shanks(seq)[-1][-1] always gives an estimate of the limit.

Optionally, an existing table can be passed to shanks(). This can be used to efficiently extend a previous computation after new elements have been appended to the sequence. The table will then be updated in-place.

The Shanks transformation:

The Shanks transformation is defined as follows (see [2]): given the input sequence (A_0, A_1, ldots), the transformed sequence is given by

\[S(A_k) = \frac{A_{k+1}A_{k-1}-A_k^2}{A_{k+1}+A_{k-1}-2 A_k}\]

The Shanks transformation gives the exact limit A_{infty} in a single step if A_k = A + a q^k. Note in particular that it extrapolates the exact sum of a geometric series in a single step.

Applying the Shanks transformation once often improves convergence substantially for an arbitrary sequence, but the optimal effect is obtained by applying it iteratively: S(S(A_k)), S(S(S(A_k))), ldots.

Wynn’s epsilon algorithm provides an efficient way to generate the table of iterated Shanks transformations. It reduces the computation of each element to essentially a single division, at the cost of requiring dummy elements in the table. See [1] for details.

Precision issues:

Due to cancellation effects, the sequence must be typically be computed at a much higher precision than the target accuracy of the extrapolation.

If the Shanks transformation converges to the exact limit (such as if the sequence is a geometric series), then a division by zero occurs. By default, shanks() handles this case by terminating the iteration and returning the table it has generated so far. With randomized=True, it will instead replace the zero by a pseudorandom number close to zero. (TODO: find a better solution to this problem.)

Examples (truncated from original)

We illustrate by applying Shanks transformation to the Leibniz series for pi:

>>> S = [4*sum((-1)**n/(2*n+1) for n in range(m))
... for m in range(1,30)]
>>>
>>> T = shanks_table(S[:7])
>>> for row in T:
...  print('['+', '.join(['{:.6g}'.format(v) for v in row])+']')
...
[-0.75]
[1.25, 3.16667]
[-1.75, 3.13333, -28.75]
[2.25, 3.14524, 82.25, 3.14234]
[-2.75, 3.13968, -177.75, 3.14139, -969.938]
[3.25, 3.14271, 327.25, 3.14166, 3515.06, 3.14161]