Found an interesting simple approximation for the integral of a function
06-15-2021, 09:16 PM (This post was last modified: 06-18-2021 03:13 AM by Namir.)
Post: #1
 Namir Senior Member Posts: 1,032 Joined: Dec 2013
Found an interesting simple approximation for the integral of a function
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

Namir
06-18-2021, 12:01 AM
Post: #2
 ttw Member Posts: 275 Joined: Jun 2014
RE: Found an interesting simple approximation for the integral of a function
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-18-2021, 12:52 AM (This post was last modified: 06-18-2021 11:53 AM by Albert Chan.)
Post: #3
 Albert Chan Senior Member Posts: 2,559 Joined: Jul 2018
RE: Found an interesting simple approximation for the integral of a function
(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
06-18-2021, 05:13 AM (This post was last modified: 06-18-2021 05:15 AM by Namir.)
Post: #4
 Namir Senior Member Posts: 1,032 Joined: Dec 2013
RE: Found an interesting simple approximation for the integral of a function
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, 12:08 PM
Post: #5
 Albert Chan Senior Member Posts: 2,559 Joined: Jul 2018
RE: Found an interesting simple approximation for the integral of a function
(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, 01:29 PM (This post was last modified: 06-18-2021 05:12 PM by Namir.)
Post: #6
 Namir Senior Member Posts: 1,032 Joined: Dec 2013
RE: Found an interesting simple approximation for the integral of a function
(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
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)