New Quadratic Integration

[Namir]

If you are familiar with the basic Simpson's rule, you know that it integrates f(x) in an interval (a, b). This algorithm calculates f(x) at a, b, and (a+b)/2. In other words, the Simpson's rule samples f(x) at midpoint of (a, b). Typically, we use the chained version of Simpson's rule, where we use multiple equidistant values of x to calculate the integral.

I decided to develop a similar algorithm that calculates the integral of f(x) for the interval (a, b). Instead of calculating f((a + b)/2), I decided to determine the value of x=m where f(m) = (f(a) + f(b)) / 2. In other words, I seek the value of x where f(x) is the average of f(a) and f(b). The first step is to perform a simple inverse Lagrangian interpolation to calculate m. Then perform an analytical integration for the Lagrangian polynomial that passes through (a, f(a), (m, f(m), and (b, f(b)).

Here is the pseudo-code for the algorithm:

Quote:Given a, b, and function f(x).
To calculate the Area of f(x) between a and b.

c = (a + b) / 2

let f(m) = (f(a) + f(b)) / 2

m = a * [(f(m) - f(c))(f(m) - f(b))]/[(f(a) - f(c))(f(a) - f(b))] +
c * [(f(m) - f(a))(f(m) - f(b))]/[(f(c) - f(a))(f(c) - f(b))] +
b * [(f(m) - f(a))(f(m) - f(b))]/[(f(b) - f(a))(f(b) - f(c))]

P(x) = f(a)*(x-m)*(x-b)/[(a-m)*(a-b)] +
f(m)*(x-a)*(x-b)/[(m-a)*(m-b)] +
f(b)*(x-a)*(x-m)/[(b-a)*(b-m)]

= Ca*(x^2-(m+b)x + m*b) +
Cm*(x^2-(a+b)x + a*b) +
Cb*(x^2-(m+a)x + m*a)

= (Ca + Cm + Cb)*x^2 -
(Ca*(m+b) + Cm*(a+b) + Cb*(m+a))*x +
(Ca*m*b + Cm*a*b + Cb*m*a)

= K1*x^2 - K2*x + K3

Area = K1*X^3/3 - K2*x^2/2 + K3*x ] from a to b
= K1*(b^3-a^3)/3 - K2*(b^2-a^2)/2 + K3*(b-a)
=[b-a]*[K1/3*(b^2+ ab + a^2) -
K2/2 * (b + a) + K3]

Where,

K1 = (Ca + Cb + Cb)
K2 = (Ca*(m+b) + Cm*(a+b) + Cb*(m+a))
K3 = (Ca*mb + Cm*ab + Cb*ma)

Ca = f(a)/[(a-m)(a-b)]
Cm = f(m)/[(m-a)(m-b)]
Cb = f(b)/[(b-a)(b-m)]

Here is an implementation of the function Area in Excel VBA:

Quote:Function f(ByVal X As Double) As Double
' Sample function. Edit the code to customize for other mathematical functions
f = 1 / X
End Function

Function Area(ByVal A As Double, ByVal B As Double, ByVal N As Integer)
Dim Sum As Double
Dim Blast As Double, Fa As Double, Fb As Double, Fc As Double, Fm As Double
Dim Ca As Double, Cb As Double, Cm As Double
Dim K1 As Double, K2 As Double, K3 As Double
Dim Incr As Double, m As Double, C As Double

Blast = B
Incr = (B - A) / N
Sum = 0
Do
B = A + Incr
C = (A + B) / 2
Fa = f(A)
Fb = f(B)
Fc = f(C)
Fm = (Fa + Fb) / 2

m = A * ((Fm - Fc) * (Fm - Fb)) / ((Fa - Fc) * (Fa - Fb)) + _
C * ((Fm - Fa) * (Fm - Fb)) / ((Fc - Fa) * (Fc - Fb)) + _
B * ((Fm - Fa) * (Fm - Fc)) / ((Fb - Fa) * (Fb - Fc))

Fm = f(m)

Ca = Fa / (A - m) / (A - B)
Cm = Fm / (m - A) / (m - B)
Cb = Fb / (B - A) / (B - m)

K1 = (Ca + Cm + Cb)
K2 = (Ca * (m + B) + Cm * (A + B) + Cb * (m + A))
K3 = (Ca * m * B + Cm * A * B + Cb * m * A)

Sum = Sum + Incr * (K1 / 3 * (B ^ 2 + A * B + A ^ 2) - K2 / 2 * (B + A) + K3)

A = B
Loop Until A >= Blast
Area = Sum

End Function

The parameters for function Area are:

1) A and B that define the integration interval.
2) N is the number of divisions within the integration intervl. Can be any positive value.

I tested the above code with f(x)=x and f(x)=1/x and obtained good results.

My rational behind this algorithm is the integration of function that experience steep changes. That i why I seek the value of x=m instead of x+(a+b)/2.

Enjoy!

Namir