nsum#
- scipy.integrate.nsum(f, a, b, *, step=1, args=(), log=False, maxterms=1048576, tolerances=None)[source]#
Evaluate a convergent, monotonically decreasing finite or infinite series.
For finite b, this evaluates:
f(a + np.arange(n)*step).sum()
where
n = int((b - a) / step) + 1
, where f is smooth, positive, and monotone decreasing. b 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 andargs
is a tuple, which may contain an arbitrary number of arrays that are broadcastable withx
.f must be an elementwise function: each element
f(x)[i]
must equalf(x[i])
for all indicesi
. It must not mutate the arrayx
or the arrays inargs
, and it must return NaN where the argument is NaN.f must represent a smooth, positive, and monotone decreasing function of x defined at all reals between a and b.
nsum
performs no checks to verify that these conditions are met and may return erroneous results if they are violated.- a, bfloat array_like
Real lower and upper limits of summed terms. Must be broadcastable. Each element of a must be finite and less than the corresponding element in b, but elements of b may be infinite.
- 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 thatf(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
, whereeps
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 attributes. (The descriptions are written as though the values will be scalars; however, if func returns an array, the outputs will be arrays of the same shape.)- successbool
True
when the algorithm terminated successfully (status0
);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.
- 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 func was evaluated.
See also
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
.Although the callable f must be non-negative and monotonically decreasing,
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, maxterms=1e3) >>> ref = np.pi**2/6 # true value >>> res.error # estimated error 4.990061275730517e-07 >>> (res.sum - ref)/ref # true error -1.0104163408712734e-10 >>> res.nfev # number of points at which callable was evaluated 1209
Compute the infinite sums of the reciprocals of integers raised to powers
p
, wherep
is an array.>>> from scipy import special >>> p = np.arange(2, 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 (0.6931471805598691, -7.616129948928574e-14)