Calculating Middle Points when using Simpson's Rule

Question

I've been learning Excel-VBA and trying to implement some basic functions using numerical analysis techniques. One of the things I'm working on implementing is Simpson's Rule for numerical Integrals: $\int_a^b f(x)\approx \frac{b-a}{3}\left(f(a) +4f(\frac{a+b}{2})+f(b)\right)$. My code for implementing this is:

Option Explicit

Function SimpsonIntegration(KnownXs As Variant, KnownYs As Variant) As Variant

    'Calculates the area under a curve using the Simpson rule.
    'KnownXs and KnownYs are the known (x,y) points of the curve.

    Dim i As Integer

    'Check if the X values are range.
    If Not TypeName(KnownXs) = "Range" Then
        SimpsonIntegration = "Xs range is not valid"
        Exit Function
    End If

    'Check if the Y values are range.
    If Not TypeName(KnownYs) = "Range" Then
        SimpsonIntegration = "Ys range is not valid"
        Exit Function
    End If

    'Check if the number of X values are equal to the number of Y values.
    If KnownXs.Rows.Count <> KnownYs.Rows.Count Then
        SimpsonIntegration = "Number of Xs <> Number of Ys"
        Exit Function
    End If

    SimpsonIntegration = 0

    'Apply the Simpson Rule: (x(i+2)-x(i))/6 * (y(i+2) + 4 * y(i+1) + y(i))
    For i = 1 To KnownXs.Rows.Count - 2 Step 2
        SimpsonIntegration = SimpsonIntegration + (KnownXs.Cells(i + 2, 1) - KnownXs.Cells(i, 1)) / 6 _
        * (KnownYs.Cells(i + 2, 1) + 4 * KnownYs.Cells(i + 1, 1) + KnownYs.Cells(i, 1))
    Next i

End Function

If you don't know Excel-VBA then pretty much all I'm doing is taking a list with all the equally spaced time points ($\textbf{X} = \{x_1, x_2, ..., x_n\}$) and a list with all the measured values ($\textbf{Y} = \{y_1, y_2, ..., y_n\}$). Then I'm initializing the integral to be zero and taking a running sum of the integral for every three point window. So for the first three points the integral is $\frac{x_3 - x_1}{3}\left(x_1 + 4 x_2 + x_3\right)$. This will be the integral of $\int_{x_1}^{x_3}\bf{Y}$. Then I do the same thing for $\int_{x_3}^{x_5}\bf{Y}, ...., \int_{x_{n-2}}^{x_{n}}\bf{Y}$. My issue is that I'm trying to plot the running integral against $\bf{X}$ on the same plot with $\bf{Y}$, but I only have points of the integral for every other point on $\bf{X}$, making the integral plot a bit choppy and only good for odd indexes of $\bf{X}$. I was looking for an overlapping version of Simpson's Rule but could not find anything. Is there any way to calculate these middle points or any overlapping version of Simpson's Rule?

Answer 1

Since you did not include a plot, I can only guess what you want. And I believe, that you actually want a good way to plot what the Simpson's Rule does.

The Simpson's Rule represents $f$ with a quadratic polynomial $p$ on each interval. Then the integral is approximated by $$\int_If(x)\mathrm{d}x \approx ∫_Ip(x)\mathrm{d}x.$$

You currently plot function values and their integral value, which (IMO) does not show anything.
Just imagine a plot with $f(x)=x$ and $F(x)=\frac{1}{2}x^2$.
What is the intention of that plot? How do you interpret it?

Therefore it is more intuitive to plot $f$ and the polynomial $p$. Since the integral basically is the area between $f$ and the $x$-axis, you can directly see the connection of the real Integral $∫f\mathrm{d}x$ and the approximation $∫p\mathrm{d}x$.

---

Overlapping formula:
I have never seen that, but on paper you can do that. What you have to do, is dig into the derivation of the Simpson's Rule. I don't know if you will get a nice formulation, but here is the rough idea.

On the interval $[x_0,x_2]$, with mid-point $x_1$, the polynomial reads: \begin{align*}P_{012}(x) =&\ f(x_0)\underbrace{\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}}_{a_0}\\ &+f(x_1)\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}\\ &+f(x_2)\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)} \end{align*}

That is just the quadratic Lagrange-Interpolation of $f$ using the nodal points $(x_0,x_1,x_2)$. Integrating this polynomial over $[x_0,x_2]$ yields the known Simpson's Rule:

$$I^{(2)}(f)=\frac{x_2-x_0}{6}\left(f(x_0) + 4f(x_1) + f(x_2)\right).$$

The integration can be done quite elegantly by transforming from the interval $[x_0,x_2]$ to $[0,2]$. For example when you integrate the first summand of $P_{012}$, you will get: \begin{align*}∫_{x_0}^{x_2}a_0\mathrm{d}x &= H\int_0^2\frac{x-1}{0-1}\frac{x-2}{0-2}\mathrm{d}x \\ &=\frac{H}{2}\int_{0}^{2}(x^2-3x+2)\mathrm{d}x = \frac{1}{3}H, \end{align*} with $H=(x_2-x_0)/2$.

We just calculated the first summand of the Simpson's Rule: $$I^{(2)}(f)=\frac{H}{3}f(x_0) + …$$

And everything we just did is basically the approach of the so-called Newton-Cotes formulas.

In your case, you have an overlapping integration of the two functions $P_{012}$ and $P_{123}$ on the interval $[x_0,x_3]$. Here you could do the integration like:

$$\int_{x_0}^{x_3}f(x)\mathrm{d}x\approx \int_{x_0}^{x_1}P_{012}(x)\mathrm{d}x + \frac{1}{2}\int_{x_1}^{x_2}P_{012}(x)+P_{123}(x)\mathrm{d}x +\int_{x_2}^{x_3}P_{123}(x)\mathrm{d}x $$

---

Edit: Maybe just a typo in your question: The Simpson rule has a $\frac{1}{6}$ not a $\frac{1}{3}$ factor.