Found an interesting simple approximation for the integral of a function

[Namir]

Hi,

I found in the 9th edition of "Numerical Analysis" by Burden & Faires, a question (number 16) on page 202 about estimating an integral using:

Integral of f(x) from a to b = 9/4*h*f(x1) + 3/4*h*f(x2)

Where h = (b-a)/3, x1 = a + h, and x2 = b.

You can divide your interval (A, B) into a series of small (a, b) and reuse the above equation. Here is a Python implementation that uses the above equation is a series of small intervals:

Code:

from math import *

def fx(x):
  return 1/x
  
def Integral(fx,A,B,h):
  a = A
  b = a + h
  s = 0
  h2=h/3
  t=3/4*h2
  while a < B:
    s += t*(3*fx(a+h2)+fx(b))
    a = b
    b += h
  return s
  
def Exact(a,b):
  return log(b/a)
  
a = 1
b = 100
h = 0.01
integ =Integral(fx,a,b,h)
print("Area=",integ)
ex = Exact(a,b)
print("%error=", 100*(ex-integ)/ex)

The output is:

Code:

Area= 4.605170176738395
%error= 2.008546184997127e-07

Does anyone know about this rather unusual algorithm?

Namir

[ttw]

The geometry seems to fit a trapezoid so I'd guess the rule is quadratic. I did plug in various values for general polynomials (X^k is sufficient) and used (the interval 0 to 1 is sufficient) and tried various values of k. If I did it correctly (not guaranteed, of course), X^2 would be integrated correctly but not X^3. The midpoint rule and composite trapezoidal rules do the same for fewer function evaluations.

[Albert Chan]

(06-15-2021 09:16 PM)Namir Wrote:  Integral of f(x) from a to b = 9/4*h*f(x1) + 3/4*h*f(x2)

Where h = (b-a)/3, x1 = a + h, and x2 - b.

If we extend for more intervals, we have weight of

0301301301301 ...

Flip the limits, integrating from right to left, weight also flipped

1031031031031 ...

∫(f(x), x=a .. b) = ∫(f(x) + f(a+b-x), x=a .. b)/2

RHS folded integral, we add the weights:

1332332332332 ...

Folded integral is Simpson's 3/8 rule

[Namir]

I was curious about the name of the algorithm, since the book authors just tossed the equation in one of their questions without giving it a name. I do recognize that it is a variant of teh trapezoidal integration. I was curious that it calculated fx(x) at the end of the interval b and at x= a +(b-a)/3, giving that point more weight.

In the absence of a name, I am calling this type of algorithm, the Quasi Trapezoidal method.

I reversed the weights calculating fx(x) at a and at b-(b-a)/3. Here is the updated Micro Python code, for bother versions, that I ran on a TI NSpire emulator:

Code:

from math import *

def fx(x):
  return 1/x
  
def QuasiTrapezoidalInteg1(fx,A,B,h):
  a = A
  b = a + h
  s = 0
  h2=h/3
  t=3/4*h2
  while a < B:
    s += t*(3*fx(a+h2)+fx(b))
    a = b
    b += h
  return s
  
def QuasiTrapezoidalInteg2(fx,A,B,h):
  a = A
  b = a + h
  s = 0
  h2=h/3
  t=3/4*h2
  while a < B:
    s += t*(fx(a)+3*fx(b-h2))
    a = b
    b += h
  return s
    
def Exact(a,b):
  return log(b/a)
  
a = 1
b = 100
h = 0.01
integ =QuasiTrapezoidalInteg1(fx,a,b,h)
print("Area=",integ)
ex = Exact(a,b)
print("%error=", 100*(ex-integ)/ex)
integ =QuasiTrapezoidalInteg2(fx,a,b,h)
print("Area=",integ)
print("%error=", 100*(ex-integ)/ex)

The output is:

Code:


Area= 4.605170176738395
%error= 2.008546184997127e-07
Area= 4.605170195256284
%error= -2.012562416029829e-07

Shows that both versions, as expected, give close results.

[Albert Chan]

(06-18-2021 05:13 AM)Namir Wrote:  In the absence of a name, I am calling this type of algorithm, the Quasi Trapezoidal method.

On my previous post (now corrected), I messed up sum of weight.
Folded integral is really Simpson 3/8 rule.

Perhaps we should name this alogirhtm Quasi Simpson 3/8 method

[Namir]

(06-18-2021 12:08 PM)Albert Chan Wrote:  

(06-18-2021 05:13 AM)Namir Wrote:  In the absence of a name, I am calling this type of algorithm, the Quasi Trapezoidal method.

On my previous post (now corrected), I messed up sum of weight.
Folded integral is really Simpson 3/8 rule.

Perhaps we should name this alogirhtm Quasi Simpson 3/8 method

Sounds like you are talking about combining the two versions of the Quasi-Trapezoidal method to end up with a four-point Quasi Simpson 3/8 method[.

I went ahead and combined the two versions Quasi Trapezoidal methods and indeed came out with something that looks like SImpsin's 3/8 rule. Here is the Mini Python code for the whole thing:

Code:

from math import *

def fx(x):
  return 1/x
  
def QuasiTrapezoidalInteg1(fx,A,B,h):
  a = A
  b = a + h
  s = 0
  h2=h/3
  t=3/4*h2
  while a < B:
    s += t*(3*fx(a+h2)+fx(b))
    a = b
    b += h
  return s
  
def QuasiTrapezoidalInteg2(fx,A,B,h):
  a = A
  b = a + h
  s = 0
  h2=h/3
  t=3/4*h2
  while a < B:
    s += t*(fx(a)+3*fx(b-h2))
    a = b
    b += h
  return s
    
def QuasiSimpsom38(fx,A,B,h):
  a = A
  b = a + h
  s = 0
  h2=h/3
  t=3/4*h2
  while a < B:
    s += t*(fx(a)+3*(fx(a+h2)+fx(b-h2))+fx(b))
    a = b
    b += h
  return s/2    
    
def Exact(a,b):
  return log(b/a)
  
a = 1
b = 100
h = 0.01
integ =QuasiTrapezoidalInteg1(fx,a,b,h)
print("Area=",integ)
ex = Exact(a,b)
print("%error=", 100*(ex-integ)/ex)
integ =QuasiTrapezoidalInteg2(fx,a,b,h)
print("Area=",integ)
print("%error=", 100*(ex-integ)/ex)
integ =QuasiSimpsom38(fx,a,b,h)
print("Area=",integ)
print("%error=", 100*(ex-integ)/ex)

The output is:

Code:


Area= 4.605170176738395
%error= 2.008546184997127e-07
Area= 4.605170195256284
%error= -2.012562416029829e-07
Area= 4.605170185997313
%error= -2.002329551551245e-10

Showing that the Quasi Simpson 3/8 does better than the other two methods, as expected.

Namir