Accurracy of 10-point Newton-Cotes formula

[Bernard]

Hi everyone,

my problem is not diectly related to MKL software and I think that this foum is more appropriate for asking question about the numerical integration.I 'am developing a numerical integration library which is based on well known multipoint Newton-Cotes formulas taken directly from Mathworld article.I have implemented already most of the formulas and ran simple accurracy related tests which gave a satisfactional results.In  case I have encountered a problem which is probably related to catastrophic cancelleation because of alterating values(points) beign used to achieve the convergence of Integral.Can anyone help me?I would like to know how to minimize the error which is always ~0.2.

Thanks in advance.

[mecej4]

As a test case, use F(x)=1 and see if the integral comes out equal to b-a. Next, try F(x) = x, F(x) = 1-x, etc and see if the result is exact for F(x) a polynomial of sufficiently low degree.

I suspect that there is at least one error in the integration formulae as given in the code above, since with F(x) = 1, a = -2, b = +2 the code returns the value of the integral as 3.631944 instead of 4.

[Dmitry_B_Intel]

To reduce the effect of cancellation it makes sense to try Kahan summation.

[Bernard]

Hi mecej4,

What do you mean by using F(x) = 1 , F(x) = x and F(x) = 1-x ?

The function beign integrated is normal distribution density function which is smooth at small interval [0-3] and I do not try to integrate over a larger domain because of very fast function roll off.The problem lies in 9-point quadrature formula which has negative and positive points.There is no any problem with the other quadrature.

[Bernard]

@mecej4

If you are interested I will upload a properly functioning formula 7-point quadrature.As I said earlier other quadratures are working properly.

[mecej4]

Most approximate quadrature rules give exact results when applied to integrands that are polynomials of sufficiently low degree. Even the trapezoidal rule gives exact results for these polynomials: (i) F(x) =1 (degree 0) (ii) F(x) = x (degree 1).

Simpson's rule is exact for polynomials of degree 3 or less.

Your 7-point routine is also in error. When invoked with these lines of code

[cpp] double fun(double x){ return 1.0; } int main(){ double a=-2.0, b=2.0; printf("%e",Quadrature7Points(fun,a,b,16)); } [/cpp]

it gives 3.789474 instead of the expected value of 4.

[Bernard]

I must admit that I only tested trapezoidal ,various simpson and multi-point quadratures on normal density distribution function and the results were correct.Can you test it with normal density distribution function on interval -3,3 and 0,3? I think that area of F(x) = 1 and F(x) = x would be better approximated by trapezoidal quadratures and area under some polynomial curve would be better approximated by the simpson and higher order quadratures.Maybe usage of high order quadrature is an "overkill" for the example given by you?I mean the usage of parabolas to approximate rectangular or trapezoidal areas.

[mecej4]

There is no point in testing a quadrature routine on a complicated function such as the normal probability density function when the routine is already known to fail for simpler functions such as F(x) = 1.

While it may be overkill to apply, for example, Simpson's rule, to a linear or quadratic function, such functions are special cases of a cubic function and, therefore, the integration routine should handle such special cases correctly.

You need to debug your programs and ensure that the C codes represent the Newton-Cotes composite rules correctly. Here, for example, is a corrected version of your 7-point code, with n=number of panels into which the interval is divided:

[cpp] if(func == NULL || n == 0)return 0.0; h = (b-a)/(6*n); sum = 41.0*(-func(a) + func(b)); for(i = 0; i < n; i++){ x = a+6*h*i; sum += 82.0*func(x) + 216.0*(func(x+h) + func(x+5*h)) + 27.0*(func(x+2*h) + func(x+4*h)) + 272.0*func(x+3*h); } return h*sum/140; [/cpp]

[Bernard]

@mecej4 Thanks for your help I really appreciate it.Does the quadrature corrected by you properly calculate simple areas like F(x) = 1 and F(x) = x? Yes I did not perform an extensive testing and debugging.I saw the corrected version and you removed the if-else block and the sum is initialized to the crudest calculation of the integrand area.My question is (because I'm not able to test it now) does the improper initialization of the sum before the main loop could lead to the wrong result? Thanks for your help

[mecej4]

Does the quadrature corrected by you properly calculate simple areas like F(x) = 1 and F(x) = x?

Yes. It is trivial to check that for yourself, as well.

the sum is initialized to the crudest calculation of the integrand area

That is your gratuitous interpretation, and an incorrect one, as well: note the minus sign preceding func(a).That is to compensate for F(a) being counted twice in the for loop for i = 0, as if x = a was also a point at which two panels met.

To see if a particular code is correct, write down the composite Newton-Cotes of your desired order with, say, two panels. With n = 2 in the code, write down the expression for the integral in terms of the function values and coefficients. Compare the two.

[Bernard]

An extensive  testing of various qudratures is planned and I will post the results.I need to solve the issue when the multipoint quadratureworked for the complicated function and failed for thesimple area.I still think that higher order polynomial are not appropriate for calculation of the simple rectangle and trapezoidal areas.I simply needto test it.

[mecej4]

I still think that higher order polynomial are not appropriate for calculation of the simple rectangle and trapezoidal areas.

From a computational point of view, or an aesthetic one, yes. From a mathematical point of view, however, any polynomial of degree M is a particular instance of polynomials of degree N greater than M, and any code that implements Newton-Cotes formula that is exact for polynomials of degree N should also work exactly for any lower degree, assuming that functions contained in the integrand can be evaluated with an error of no more than a few ULPs.

To see this, look up the error term at the Mathworld page that you cited, for the 7-point formula. The error term is O(h⁸) F⁽⁸⁾ (ξ), for ξ between a and b. If F is a polynomial of any degree from 0 to 7, the 8th derivative is zero and the Newton-Cotes formula is exact.

[Bernard]

@mecej4

Thank you for the explanation.

[Bernard]

Hi  mecej4,

I have an update.As you advised I tested my newest quadrature function(alternating points) based on 11-point Newton-Cotes formula and I was able to achieve an accuracy of 5 decimal places when F(x) = 1 [-2,2] was evaluated , albeit at very high price of more than 1000 panels.I suppose that approximation of rectangle area by higher order polynomial and rounding errors contributed to some loss of accuracy.I must say that when quadrature is called with less than 100 panels the error is not acceptable ~0.20.

[mecej4]

Once again you are blaming the quadrature rule itself when the error is in your implementation. I suggest that you select, say, the 3 or 5-point Newton-Cotes rule to keep things simple, and perform the calculation by hand using a single panel (i.e., n = 1). For F(x) = 1, or F(x) = x, or F(x) = x², you should obtain the exact result. Then, code the same calculation and run the code with various small values of n. Your code should also give the exact result. If not, fix the code first. It is a mistake to attempt to make the error go away by using ridiculously large values of n or by using higher order rules such as the 11-point rule.

The attached file demonstrates that the correctly implemented 11-point rule matches the exact result for two cases: F(x)=1; F(x)=normal p.d.f. For both cases, a single panel is used (n = 1).

[Bernard]

HI mecej4,

thanks for your answer.I must admit that I was slightly confused when started to implement multipoint quadratures,because for complicated functions they were more accurate than simpson and/or trapezoid ,but failed for simple integrals .3-point quadratures ie simpson composite and n-point simpson 3/8 were also tested and the results were qiute accurate for simple integrals,so as you stated the problem lies in my implementation and i was lazy to do integration by hand.I saw your corrected version and the main difference lies in initialization of the sum variable and initialization of int i = 0.

Btw. I also implemented 3-point quadrature with single panel and it gave satisfactionary results.I suppose that my error was due to inacurrate implementation when the interval was subdivided into hundreds and thousands of panels.

I will corrrect m implementation exactly as you did and run more extensive tests.

Thank you for your help and comments.