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:
:
(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*}](http://www.holoborodko.com/pavel/wp-content/ql-cache/quicklatex.com-5605e794ac90baac35e5ca0d71b174e4_l3.png "Rendered by QuickLaTeX.com")
Composite scheme divides whole integration interval by non-overlapped 
smaller intervals and uses simple Newton-Cotes rule for each of them.
Although this technique gives good results for appropriately chosen 
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 
and outer function values as 
. 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 \]](http://www.holoborodko.com/pavel/wp-content/ql-cache/quicklatex.com-dcd86f943cae7466342db1621c5a7212_l3.png "Rendered by QuickLaTeX.com")
Here are several quadratures of such type for 
:
![\[ \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} \]](http://www.holoborodko.com/pavel/wp-content/ql-cache/quicklatex.com-70dcf8d55696a9bb399555b0d627d793_l3.png "Rendered by QuickLaTeX.com")
Now, in composite scheme we can apply overlapped quadratures for approximating integral on 
intervals instead of lower-order Simpson’s 3/8 rule.
Composite rules with enhanced precision without increase in function evaluations:
: (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*}](http://www.holoborodko.com/pavel/wp-content/ql-cache/quicklatex.com-11e239bf144473f9c6e2d7080bb7f730_l3.png "Rendered by QuickLaTeX.com")
: (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*}](http://www.holoborodko.com/pavel/wp-content/ql-cache/quicklatex.com-94c641248f6ed9860c21b4e9b65b2a4e_l3.png "Rendered by QuickLaTeX.com")
: (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*}](http://www.holoborodko.com/pavel/wp-content/ql-cache/quicklatex.com-77ad45025db6aaa59af1d947734a8d76_l3.png "Rendered by QuickLaTeX.com")
Notice that we use extra nodes from neighboring intervals for estimation on every step . This allows us to improve precision on interior intervals (first and last intervals still have ). 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 
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 .
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.
(4 votes, average: 5.00)
Loading...
7 Comments
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
Yep, this is mistake – already fixed. Thank you.
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:)
Hi Bernard,
You don’t need the first loop. By
I meant that 
has only two values 
and 
, not a range.
So you don’t need first loop.
Hi Pavel,
I really appreciate your help thank you.
Bernard
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.
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