Found an interesting simple approximation for the integral of a function

[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