scipy.integrate.

nsum#

scipy.integrate.nsum(f, a, b, *, step=1, args=(), log=False, maxterms=1048576, tolerances=None)[source]#

Evaluate a convergent finite or infinite series.

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:
fcallable

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, bfloat 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.

stepfloat 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.

argstuple 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.

logbool, 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.

maxtermsint, 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, rtolfloat, 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.)

successbool

True when the algorithm terminated successfully (status 0); False otherwise.

statusint array

An integer representing the exit status of the algorithm.

  • 0 : The algorithm converged to the specified tolerances.

  • -1 : Element(s) of a, b, or step are invalid

  • -2 : Numerical integration reached its iteration limit; the sum may be divergent.

  • -3 : A non-finite value was encountered.

  • -4 : The magnitude of the last term of the partial sum exceeds the tolerances, so the error estimate exceeds the tolerances. Consider increasing maxterms or loosening tolerances. Alternatively, the callable may not be unimodal, or the limits of summation may be too far from the function maximum. Consider increasing maxterms or breaking the sum into pieces.

sumfloat array

An estimate of the sum.

errorfloat array

An estimate of the absolute error, assuming all terms are non-negative, the function is computed exactly, and direct summation is accurate to the precision of the result dtype.

nfevint array

The number of points at which f was evaluated.

See also

mpmath.nsum

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) dx \leq \sum_{k=u}^\infty f(k) \leq \int_u^\infty f(x) dx + f(u)\]

Let \(a\) represent a, \(n\) represent maxterms, \(\epsilon_a\) represent atol, and \(\epsilon_r\) represent rtol. The implementation first evaluates the integral \(S_l=\int_a^\infty f(x) dx\) as a lower bound of the infinite sum. Then, it seeks a value \(c > a\) such that \(f(c) < \epsilon_a + S_l \epsilon_r\), if it exists; otherwise, let \(c = a + n\). Then the infinite sum is approximated as

\[\sum_{k=a}^{c-1} f(k) + \int_c^\infty f(x) dx + 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.

References

[1]

Wikipedia. “Integral test for convergence.” https://en.wikipedia.org/wiki/Integral_test_for_convergence

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
np.float64(7.448762306416137e-09)
>>> (res.sum - ref)/ref  # true error
np.float64(-1.839871898894426e-13)
>>> res.nfev  # number of points at which callable was evaluated
np.int32(8561)

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)
True

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
(np.float64(0.6931471805598691), np.float64(-7.616129948928574e-14))
On this page