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

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

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

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.

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

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

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

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.

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

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

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.

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

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.

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

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]

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

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

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.

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

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.

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

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.

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

@mecej4

Thank you for the explanation.

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

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.

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

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

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

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