Overlapped Newton-Cotes Quadratures

The goal of this paper is to build numerical integration quadratures for equidistant data which use function values outside interval of integration. 

One of the main applications for such formulas is to improve accuracy of compound rules without increase in number of function evaluations. In particular we propose refined composite algorithm which deliver higher precision than Simpson’s 3/8 using the same number of function values.

Reference formulas and their properties are included.

Simpson’s 3/8 composite rule is perhaps the most used method for numerical integration of equally-spaced data.
It has fifth order accuracy and simple computational structure:

O(h^5):

(1)   \begin{equation*} \int_{x_0}^{x_{3m}} f(x) \mathrm{d}x \approx \frac{3\,h}{8}\sum_{k=1}^{m}\,\left[f_{3k-3}+3\,(f_{3k-2}+f_{3k-1})+f_{3k}\right] \end{equation*}

Composite scheme divides whole integration interval by non-overlapped m smaller intervals and uses simple Newton-Cotes rule for each of them.
Although this technique gives good results for appropriately chosen h it seems that composite rule doesn’t use full potential of numerous interior function samples. They could contribute more to the overall precision of the algorithm.

Our goal is to improve accuracy of Simpson’s 3/8 composite rule without increase of function evaluations.

With this purpose we have build special quadratures which use function values beyond integration interval. Lets denote total number of nodes in quadrature as N and outer function values as L. Then we have

    \[ \int_{x_L}^{x_{N-1-L}} f(x) \mathrm{d}x \approx \sum_{i=0}^{N-1}c_i\,f_i \]

Here are several quadratures of such type for L=\frac{N}{2}-2:

    \[ \begin{tabular}{|r|c|l|c|} \hline $N$& Overlapped Newton-Cotes Quadratures & Error Term & $L_1$ Norm\\ \hline &&&\\ 6&$\displaystyle{\frac{3\,h}{160}\,\left( -f_{{0}}+23\,f_{{1}}+58\,f_{{2}}+58\,f_{{3}}+23\,f_{{4}}-f_{{5}} \right)}$&$\displaystyle{\frac{13\,{h}^{7}}{2240}\,f^{(6)}(\xi)}$&1.03 \\ &&&\\ 8&$\displaystyle{\frac{h}{4480}\, \left[ \begin{matrix} 13\,(f_{{0}}+f_{{7}})-149\,(f_{{1}}+f_{{6}})+\\ 2049\,(f_{{2}}+f_{{5}})+4807\,(f_{{3}}+f_{{4}}) \end{matrix} \right] }$&$\displaystyle{\frac{7\,{h}^{9}}{6400}\,f^{(8)}(\xi)}$&1.04\\ &&&\\ 10&$\displaystyle{\frac{h}{89600}\, \left[ \begin{matrix} -49\,(f_{{0}}+f_{{9}})+603\,(f_{{1}}+f_{{8}})\\ -3960\,(f_{{2}}+f_{{7}})+42352\,(f_{{3}}+f_{{6}})+\\ 95454\,(f_{{4}}+f_{{5}}) \end{matrix} \right] }$&$\displaystyle{\frac{443\,{h}^{11}}{1971200}\,f^{(10)}(\xi)}$&1.06\\ &&&\\ \hline \end{tabular} \]

Now, in composite scheme we can apply overlapped quadratures for approximating integral on k=2,\dots,m-1 intervals instead of lower-order Simpson’s 3/8 rule.

Composite rules with enhanced precision without increase in function evaluations:

O(h^7):

(2)    \begin{align*} \int_{x_0}^{x_{3m}} f(x) \mathrm{d}x &\approx \frac{3\,h}{8}\sum_{k=1,\,m}\,\left[f_{3k-3}+3\,(f_{3k-2}+f_{3k-1})+f_{3k}\right]\\ &+\frac{3\,h}{160}\sum_{k=2}^{m-1}\,\left[-f_{3k-4}+23\,f_{3k-3}+58\,f_{3k-2}+58\,f_{3k-1}+23\,f_{3k}-f_{3k+1}\right] \end{align*}

O(h^9):

(3)    \begin{align*} \int_{x_0}^{x_{3m}} f(x) \mathrm{d}x &\approx \frac{3\,h}{8}\sum_{k=1,\,m}\,\left[f_{3k-3}+3\,(f_{3k-2}+f_{3k-1})+f_{3k}\right]\\ &+\frac{h}{4480}\sum_{k=2}^{m-1}\,{ \left[ \begin{matrix} 13\,(f_{{3k-5}}+f_{{3k+2}})-149\,(f_{{3k-4}}+f_{{3k+1}})+\\ 2049\,(f_{{3k-3}}+f_{{3k}})+4807\,(f_{{3k-2}}+f_{{3k-1}}) \end{matrix} \right]} \end{align*}

O(h^{11}):

(4)    \begin{align*} \int_{x_0}^{x_{3m}} f(x) \mathrm{d}x &\approx \frac{3\,h}{8}\sum_{k=1,\,m}\,\left[f_{3k-3}+3\,(f_{3k-2}+f_{3k-1})+f_{3k}\right]\\ &+\frac{h}{89600}\sum_{k=2}^{m-1}\,{ \left[ \begin{matrix} -49\,(f_{{3k-6}}+f_{{3k+3}})+603\,(f_{{3k-5}}+f_{{3k+2}})\\ -3960\,(f_{{3k-4}}+f_{{3k+1}})+42352\,(f_{{3k-3}}+f_{{3k}})+\\ 95454\,(f_{{3k-2}}+f_{{3k-1}}) \end{matrix} \right]} \end{align*}

Notice that we use extra nodes from neighboring intervals for estimation on every step k=2,\dots,m-1. This allows us to improve precision on interior intervals (first and last intervals still have O(h^5)). Most importantly, proposed high precision composite rules use the same number of nodes as Simpson’s 3/8 does. They can be plugged in place of Simpson’s 3/8 rules seamlessly in any application.

Special attention should be payed to the stability of the proposed formulas. As you see overlapped quadratures have negative coefficients and L_1 norm is bigger than 1.0 as a result. This usually indicates possible instability of the quadratures in some cases. However their magnitude responses show that they are stable:

Noisy high frequencies (usual indicators of instability) are suppressed, there are no other abnormalities. Besides we have tested proposed & Simpson’s 3/8 composite rules on 100+ test functions for integration from [1]. Overlapped quadratures gave more accurate result in 75% of cases without special tuning of the step size h.

This combined with the fact that actually norm is only slightly higher than 1.0 allows us to recommend refined composite schemes for applications.

References

I have used many excellent sources, but only one for the moment:
[1] H. Engels, Numerical Quadrature and Cubature, Academic Press, London, 1980.


1 Star2 Stars3 Stars4 Stars5 Stars (4 votes, average: 5.00)
Loading...

7 Comments

  1. renato
    Posted May 25, 2011 at 5:51 pm | #

    Pavel, thank you for replying! As you suggested, I’ll try overlapped formulae also. Just a question : in (2)-(4) is the first summation running for k=0,m or perhaps for k=1,m?
    If k=0,m, wouldn’t it require 3 additional evaluations (f(-3),f(-2),f(-1)) ?

    thank you,

    renato

  2. Bernard
    Posted April 9, 2013 at 8:25 pm | #

    Hi Pavel,

    I implemented overlapped Newton Cotes formula(3),but after extensive testing only for 1000 points I’m getting satisfactionary results.I suppose either wrong initialization of the accumulated sum or wrong indices of for loop contributed to this.
    Here is my code:

    double QuadratureComposite8Points(double(*func)(double),double a,double b,int n)
    {
    if(func == NULL || n <= 0)
    return 0.0;

    double term,x1,x2,sum,h,k;
    int i,ii;
    h = (b-a)/(3*n-1);

    sum = FUNC(a) + 3.0*(FUNC(a+h) + FUNC(a+2*h)) + FUNC(a+3*h);
    //sum = FUNC(-a) + FUNC(b);
    printf("Sum = %.17f\n",sum);
    for(i = 1;i < n;i++){
    x1 = a+3*h*i;
    sum += FUNC(x1) +3.0*(FUNC(x1+h) + FUNC(x1+2*h)) + FUNC(x1+3*h);
    }

    term = 3*h/8;

    term*sum;
    //sum += 13.0*(FUNC(a) + FUNC(a+7*h)) – 149.0*(FUNC(a+h) + FUNC(a+6*h)) + 2049.0*(FUNC(a+2*h) + FUNC(a+5*h)) + 4807.0*(FUNC(a+3*h) + FUNC(a+4*h));
    for(ii = 2;ii< n-1;ii++){
    x2 = a + 3*h*ii;
    sum += 13.0*(FUNC(x2) + FUNC(x2+7*h)) – 149.0*(FUNC(x2+h) + FUNC(x2+6*h)) + 2049.0*(FUNC(x2+2*h) + FUNC(x2+5*h)) + 4807.0*(FUNC(x2+3*h) + FUNC(x2+4*h));

    }

    return (h*sum)/4480;

    }

    Thanks in advance:)

    • Posted April 10, 2013 at 4:40 pm | #

      Hi Bernard,

      You don’t need the first loop. By k=1,m I meant that k has only two values 1 and m, not a range.
      So you don’t need first loop.

  3. Bernard
    Posted April 10, 2013 at 6:28 pm | #

    Hi Pavel,

    I really appreciate your help thank you.
    Bernard

  4. Bernard
    Posted April 10, 2013 at 6:34 pm | #

    Pavel ,

    so proper implementation should initialize the sum variable to this FUNC(a) + 3.0*(FUNC(a+h) + FUNC(a+2*h)) + FUNC(a+3*h); and calculate only the second loop which should start from k = 2 to m-1 ?

    P.S
    I was confused by the sigma in the 3/8 composite rule.

  5. Bernard
    Posted April 10, 2013 at 9:25 pm | #

    Hi Pavel,

    Can you post the results obtained by you when you tested Composite Rule (4) on probability density function.Unfortunately my results differ by large margin when compared to composite simpson rule.Is this possible to post also your code implementation.

    Thanks in advance
    Bernard

Post a Comment

Your email is never published nor shared.

Use native LaTeX syntax to include formulas: $ ... $, \[ ... \], etc. Do not forget to preview comment before posting.

Also you may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>

Subscribe without commenting