bundles / scipy latest / scipy / integrate / _tanhsinh / nsum

function

`scipy.integrate._tanhsinh:nsum`

source: /scipy/integrate/_tanhsinh.py :962

Signature

   def     nsum (    f    ,    a    ,    b    ,  * ,     step    =  1   ,    args    =  ()   ,    log    =  False   ,    maxterms    =  1048576   ,    tolerances    =  None     )

Summary

Evaluate a convergent finite or infinite series.

Extended Summary

For finite a and b, this evaluates

f(a + np.arange(n)*step).sum()

where n = int((b - a) / step) + 1, where f is smooth, positive, and unimodal. The number of terms in the sum may be very large or infinite, in which case a partial sum is evaluated directly and the remainder is approximated using integration.

Parameters

f : callable

The function that evaluates terms to be summed. The signature must be

f(x: ndarray, *args) -> ndarray

where each element of x is a finite real and args is a tuple, which may contain an arbitrary number of arrays that are broadcastable with x.

f must be an elementwise function: each element f(x)[i] must equal f(x[i]) for all indices i. It must not mutate the array x or the arrays in args, and it must return NaN where the argument is NaN.

f must represent a smooth, positive, unimodal function of x defined at all reals between a and b.

a, b : float array_like

Real lower and upper limits of summed terms. Must be broadcastable. Each element of a must be less than the corresponding element in b.

step : float array_like

Finite, positive, real step between summed terms. Must be broadcastable with a and b. Note that the number of terms included in the sum will be floor((b - a) / step) + 1; adjust b accordingly to ensure that f(b) is included if intended.

args : tuple of array_like, optional

Additional positional arguments to be passed to f. Must be arrays broadcastable with a, b, and step. If the callable to be summed requires arguments that are not broadcastable with a, b, and step, wrap that callable with f such that f accepts only x and broadcastable *args. See Examples.

log : bool, default: False

Setting to True indicates that f returns the log of the terms and that atol and rtol are expressed as the logs of the absolute and relative errors. In this case, the result object will contain the log of the sum and error. This is useful for summands for which numerical underflow or overflow would lead to inaccuracies.

maxterms : int, default: 2**20

The maximum number of terms to evaluate for direct summation. Additional function evaluations may be performed for input validation and integral evaluation.

atol, rtol : float, optional

Absolute termination tolerance (default: 0) and relative termination tolerance (default: eps**0.5, where eps is the precision of the result dtype), respectively. Must be non-negative and finite if log is False, and must be expressed as the log of a non-negative and finite number if log is True.

Returns

res : _RichResult

An object similar to an instance of scipy.optimize.OptimizeResult with the following attributes. (The descriptions are written as though the values will be scalars; however, if f returns an array, the outputs will be arrays of the same shape.)

success: success
status: status
sum: sum
error: error
nfev: nfev

Notes

The method implemented for infinite summation is related to the integral test for convergence of an infinite series: assuming step size 1 for simplicity of exposition, the sum of a monotone decreasing function is bounded by

\int_{u}^{\infty} f (x) d x \leq k = u \sum \infty f (k) \leq \int_{u}^{\infty} f (x) d x + f (u)

Let $a$ represent a, $n$ represent maxterms, $ϵ_{a}$ represent atol, and $ϵ_{r}$ represent rtol. The implementation first evaluates the integral $S_{l} = \int_{a}^{\infty} f (x) d x$ as a lower bound of the infinite sum. Then, it seeks a value $c > a$ such that $f (c) < ϵ_{a} + S_{l} ϵ_{r}$ , if it exists; otherwise, let $c = a + n$ . Then the infinite sum is approximated as

k = a \sum c - 1 f (k) + \int_{c}^{\infty} f (x) d x + f (c) /2

and the reported error is $f (c) /2$ plus the error estimate of numerical integration. Note that the integral approximations may require evaluation of the function at points besides those that appear in the sum, so f must be a continuous and monotonically decreasing function defined for all reals within the integration interval. However, due to the nature of the integral approximation, the shape of the function between points that appear in the sum has little effect. If there is not a natural extension of the function to all reals, consider using linear interpolation, which is easy to evaluate and preserves monotonicity.

The approach described above is generalized for non-unit step and finite b that is too large for direct evaluation of the sum, i.e. b - a + 1 > maxterms. It is further generalized to unimodal functions by directly summing terms surrounding the maximum. This strategy may fail:

If the left limit is finite and the maximum is far from it.
If the right limit is finite and the maximum is far from it.
If both limits are finite and the maximum is far from the origin.

In these cases, accuracy may be poor, and nsum may return status code 4.

Although the callable f must be non-negative and unimodal, nsum can be used to evaluate more general forms of series. For instance, to evaluate an alternating series, pass a callable that returns the difference between pairs of adjacent terms, and adjust step accordingly. See Examples.

Array API Standard Support

nsum has experimental support for Python Array API Standard compatible backends in addition to NumPy. Please consider testing these features by setting an environment variable SCIPY_ARRAY_API=1 and providing CuPy, PyTorch, JAX, or Dask arrays as array arguments. The following combinations of backend and device (or other capability) are supported.

====================  ====================  ====================
Library               CPU                   GPU
====================  ====================  ====================
NumPy                 ✅                     n/a                 
CuPy                  n/a                   ✅                   
PyTorch               ⛔                     ⛔                   
JAX                   ⛔                     ⛔                   
Dask                  ⛔                     n/a                 
====================  ====================  ====================

See dev-arrayapi for more information.

Examples

Compute the infinite sum of the reciprocals of squared integers.

import numpy as np
from scipy.integrate import nsum
res = nsum(lambda k: 1/k**2, 1, np.inf)
ref = np.pi**2/6  # true value
res.error  # estimated error
(res.sum - ref)/ref  # true error
res.nfev  # number of points at which callable was evaluated

✓

Compute the infinite sums of the reciprocals of integers raised to powers ``p``, where ``p`` is an array.

from scipy import special
p = np.arange(3, 10)
res = nsum(lambda k, p: 1/k**p, 1, np.inf, maxterms=1e3, args=(p,))
ref = special.zeta(p, 1)
np.allclose(res.sum, ref)

✓

Evaluate the alternating harmonic series.

res = nsum(lambda x: 1/x - 1/(x+1), 1, np.inf, step=2)

✓

res.sum, res.sum - np.log(2)  # result, difference vs analytical sum

✗

Aliases

scipy.integrate.nsum

Referenced by

This package

release:1.15.0-notes