Skip to main content
`Ctrl`+`K`
scipy.integrate.

# nsum

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

  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:
    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.)

      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.

  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 Support for the array API standard for more information.

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