Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

- Intel Community
- Software
- Software Development SDKs and Libraries
- Intel® oneAPI Math Kernel Library
- Accurracy of 10-point Newton-Cotes formula

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Mute
- Printer Friendly Page

Bernard

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

04-05-2013
04:29 AM

72 Views

Accurracy of 10-point Newton-Cotes formula

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.

Link Copied

15 Replies

mecej4

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

04-05-2013
05:23 AM

72 Views

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

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

04-05-2013
05:44 AM

72 Views

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

Bernard

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

04-05-2013
06:43 AM

72 Views

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

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

04-05-2013
06:50 AM

72 Views

mecej4

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

04-05-2013
08:18 AM

72 Views

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

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

04-05-2013
08:43 AM

72 Views

mecej4

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

04-05-2013
09:42 AM

72 Views

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

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

04-05-2013
01:37 PM

72 Views

mecej4

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

04-05-2013
01:48 PM

72 Views

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

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

04-05-2013
02:06 PM

72 Views

mecej4

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

04-05-2013
02:42 PM

72 Views

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^{8}) F^{(8)} (ξ), 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

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

04-06-2013
12:13 AM

72 Views

@mecej4

Thank you for the explanation.

Bernard

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

04-08-2013
05:49 AM

72 Views

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

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

04-08-2013
07:19 AM

72 Views

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^{2}, 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

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

04-08-2013
10:49 AM

72 Views

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.

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page

For more complete information about compiler optimizations, see our Optimization Notice.