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)

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

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

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:

(2)

(3)

(4)

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.

## 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 $k=1,m$ I meant that $k$ has only two values $1$ and $m$, 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