having only discrete values sampled at equidistant points with step :
where coefficients are derived from interpolating polynomial fitted for points :
Goal of this simple idea is to build interpolating polynomial of the highest degree allowed by number of nodes (=degrees of freedom in the system) and assume that . It is implied that polynomial of high degree is able to resemble behavior of more complex functions and minimizes estimation error as a result.
However this implication is not true in general case. Quality of polynomial approximation of a function depends drastically on how interpolation nodes are chosen [Trefethens, Rivling]. It is well known (see Runge’s phenomenon) that in case of equally-spaced data approximation error could in fact become worse as we increase number of interpolation nodes and use polynomial of higher order.
This makes impossible to guarantee convergence of for with the exception of very specific type of functions [Krylov 145]. Which means that high order Newton-Cotes rules tend to give inferior results compared to lower order formulas.
Spectral properties of Newton-Cotes formulas
Let’s study properties of Newton-Cotes formulas in frequency domain (or spectral properties) using DSP tools and terminology. NC rules are symmetric digital filters with finite impulse response. One of the most important characteristic of digital filter is magnitude/frequency response – function which shows how much filter damps or amplifies magnitude of particular frequency contained in input data.
Response functions of the first several NC formulas are presented below:
As we can see classical NC rules tend to suppress high frequencies (noise) in the input data and leave lowest harmonics (pure signal) almost intact. In that sense they can be called smoothing or low-pass filters.
Simpson 3/8 rule has the strongest suppression of noise (most stable). Plain Simpson has the weakest. Boole’s formula damps mid-range harmonics but gives somewhat larger response on the noisy high frequencies. Overall they can be considered as stable and acceptable for practical applications from DSP point of view.
Let’s take a look on magnitude responses of high-order Newton-Cotes rules where numerical instability shows up.
Here we have much more interesting situation. Rules are not low-pass anymore, actually they can be called “super” noise amplifiers as they strengthen high frequencies by great extend. This makes formulas numerically unstable, as rounding errors (or other “noisy” components) are amplified and might deteriorate the accuracy.
Although it is well known that high-order Newton-Cotes formulas are not suitable for applications, frequency domain gives us clear illustration for one of the “unsuitability” reasons – noise amplification. Obviously it is closely related to the characteristic of error resistance/propagation, , which is commonly used to measure stability of numerical methods in classical analysis.
To overcome this issue usually researcher switches from equidistant integration rules to one of the more advanced quadrature of Gauss type (Legendre, Laguerre, Jacobi or else). There nodes are distributed more densely towards interval boundaries which allows preventing “wild” behavior of interpolating polynomial and leads to stable computations.
However in many practical applications only uniformly spaced function values are available and stable Newton-Cotes type formulas are needed. In the following sections we propose novel approach on how to build them using least square polynomial fit instead of interpolation.
Low Noise Newton-Cotes rules
There are several ways on how to improve high-order Newton-Cotes formulas. The most obvious is to re-target some of the degrees of freedom in the system from contributing to highest approximating order to regularization and stronger noise suppression. Natural way of doing this is to use least squares approximation instead of interpolation.
On the contrary to interpolation, least squares make possible to derive several integration filters of the same approximation order. Here we consider the most required and formulae.
First two rules are classics – Simpson’s and Simpson’s 3/8. Others have much stronger noise suppression in cost of larger constant in error term. Magnitude responses are below for comparison:
Formulas of higher approximation order can be derived analogously:
All regularized Newton-Cotes rules show improved resistance to noise (leading to elevated numerical stability) in comparison with classic counterparts. Proposed formulas can be used as a basement for composite rules, adaptive integration algorithms and even more exotic overlapped compound rules.
Stable Newton-Cotes formulas of High Order
Least squares can also be used to derive numerically stable rules of high orders. Here we consider filters of . Usual NC formula for the case is highly unstable and cannot be used in practice – it amplifies noise and its norm is high . Our approach allows to overcome both issues.
Stable Newton-Cotes quadrature of needs additional six function evaluations to achieve :
Classic NC formula is marked red, proposed by green. As we can see new rule is perfectly stable.
Actually this approach allows building the whole family of stabilized rules for any particular . Case of is just the minimum which required to achieve for . All filters of are also stable (although with larger constant in error term) and usable.
If you are interested in coefficients of the stabilized filters please contact me by email/or write comment below. I omitted them since they would clatter the page.
If you find these ideas interesting please do not hesitate to contact me for further discussion.