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
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.
(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
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.
(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
(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