New Trapedoizal Tail integration method

[Namir]

Hi All,

I am sharing a new "trapezoidal tail" method. The method samples two function values in each interval at f(x+h/2-h/3) and fx(x+h/2+h/3) where h is the size of the integration interval. If you remove the h/2 part, the algorthm degrades in the results it return.

Given function fx to integrate from x=A to x=B in N even steps.

Code:

h = (B - A) / N
h3 = h / 3
X = A + h / 2
Sum = 0
Do While X < B
  Sum = Sum + fx(X - h3) + fx(X + h3)
  X = X + h
Loop
Integral = h * Sum / 2

The above algorithm performs better than similar trapezoidal integration methods, like the mid-point method or regular trapezoidal method.

How does the new algorithm compare with Simpson's rule? In general the latter methods performs better. I noticed that the absolute ratios of % errors of the new algorithm to Simpson's rule rises quadratically with the number of interval divisions. This relation seems to hold for the few different functions I tested and the curve fitted (using 5 points for each case) were perfect quadratic fits!

Namir

[Paul Dale]

I'm failing to see any advantage here. Two function evaluations per interval and first order convergence. Simpson's rule has the same number of evaluations with second order convergence.

Pauli

[Namir]

(05-13-2021 06:17 AM)Paul Dale Wrote:  I'm failing to see any advantage here. Two function evaluations per interval and first order convergence. Simpson's rule has the same number of evaluations with second order convergence.

Pauli

The algorithm is an "open" type and does not calculate values at the beginning or end of the interval. As such it can avoid runtime errors in cases where the bounds of the interval generate runtime error (i.e Simpson rule will bomb!)

Namir

[Albert Chan]

(05-12-2021 06:05 PM)Namir Wrote:  The method samples two function values in each interval at f(x+h/2-h/3) and fx(x+h/2+h/3)

You might also try f(x+h/2 ± h/(2√3))
This is the points and weights for 2 points Gaussian quadrature.

Because 2√3 ≈ 3.464 relatively close to 3, this explained the good result you observed.

[Maximilian Hohmann]

Hello!

(05-13-2021 08:28 AM)Namir Wrote:  The algorithm is an "open" type and does not calculate values at the beginning or end of the interval. As such it can avoid runtime errors in cases where the bounds of the interval generate runtime error (i.e Simpson rule will bomb!)

But what if the runtime errors occur at those 1/3 points? I guess that the probability for a runtime error will be the same for any point.

Regards
Max

[Albert Chan]

(05-13-2021 12:12 PM)Maximilian Hohmann Wrote:  But what if the runtime errors occur at those 1/3 points?
I guess that the probability for a runtime error will be the same for any point.

If singularity occurs inside integration limits, we *have* to break it up. (runtime errors or not)

This is an integral made famous by Richard Feynman.
(I saw this from Nahin's book "Inside Interesting Integrals", page 18)

\(\displaystyle \int _0 ^1 {dx \over [a\,x + b\,(1-x)]^2} = {1 \over a\,b}\)

If no singularity within the limits, integral is correct.
But, what if denominator goes zero ?

a*x + b*(1-x) = b + (a-b)*x = 0
x = b / (b-a)       // assumed a ≠ b, ok since a=b made integrand constant: 1/b^2 = 1/(ab)

If singularity occured inside integral limits:

0 < x < 1
0 < b < b-a

a < 0 and b > 0

Without break up integral to pieces, we have LHS > 0, RHS < 0, which make no sense.

[Namir]

(05-13-2021 12:05 PM)Albert Chan Wrote:  

(05-12-2021 06:05 PM)Namir Wrote:  The method samples two function values in each interval at f(x+h/2-h/3) and fx(x+h/2+h/3)

You might also try f(x+h/2 ± h/(2√3))
This is the points and weights for 2 points Gaussian quadrature.

Because 2√3 ≈ 3.464 relatively close to 3, this explained the good result you observed.

Thanks Albert! The 2*sqrt(3) makes sense. I tried several points on either side of the half interval with varying degrees of success.

Namir