scijit.integrate.newton_cotes

scijit.integrate.newton_cotes(rn, equal=0)

Weights and error coefficient of a Newton-Cotes rule.

Takes NO callback of either style: it returns the rule, it does not apply it. For an equally-spaced rule, integrate with dx * sum(an * f(a + arange(N+1) * dx)) where dx = (b - a) / N.

Parameters:

rnint, float or 1-D array_like

Either the order N, the number of intervals, so the rule uses N + 1 equally-spaced points; or the N + 1 sample positions, which must start at 0 and end at N. A float order is accepted and used by value.

equalint, optional

Set to 1 to treat the samples as equally spaced whatever rn holds, which also replaces the positions with 0..N. Default 0, and then equal spacing is detected from rn itself when every gap is exactly 1, in which case the positions are kept.

Returns:

anfloat64 array, shape (N+1,)

Weights, scaled so the integral is dx * sum(an * f_i).

Bfloat

Error coefficient. The error term is B * dx**(N+2) * f**(N+1)(xi), and B * dx**(N+3) * f**(N+2)(xi) when the samples are equally spaced and N is even.

Raises:

IndexError

An empty rn.

ValueError

The sample positions do not start at 0 and end at N; N < 1; rn of rank 2 or more.

See also

scipy.integrate.newton_cotes

The scipy routine this mirrors.

Notes

For N = 1..14 equally spaced the weights match scipy.integrate. newton_cotes to rounding and the error coefficient B is exact.

Above 14, and for unequally spaced samples, the weights are built through the Vandermonde path and lose accuracy as the system grows ill-conditioned. High-order Newton-Cotes weights alternate in sign and grow, so rules above about N = 8 are numerically unusable regardless of who computes them. Prefer composite low-order rules, romb(), or quad().

N < 1 raises here. scipy reaches rn / float(N), emits a numpy RuntimeWarning: invalid value encountered in divide, and then raises ValueError: math domain error from math.log(0).

An rn of rank 2 or more raises here. scipy reads its order from len(rn), which is the first axis alone, so the rest of the array selects nothing and the rule returned describes an order the caller did not ask for.

Pure @njit, no state, so prange-safe.

Examples

>>> import numpy as np
>>> from numba import njit
>>> import scijit.integrate as si
>>> @njit
... def run():
...     return si.newton_cotes(4, 0)
>>> an, B = run()
>>> an
array([0.31111111, 1.42222222, 0.53333333, 1.42222222, 0.31111111])
>>> B
-0.008465608465608466