L_SIMPSON

=L_SIMPSON(y, [x], [dx])

Description

L_SIMPSON is an implementation of the composite Simpson's 1/3 Rule for numerical integration of irregularly spaced data. You would use it when you want higher accuracy than the trapezoidal rule, but you still only have a set of (x,y) points rather than the ability to evaluate the function at more points.

The following image shows a set of 8 points. The black line is the actual curve which for this particular example is the function \(f(x)=0.5x^3+5x^2-10x+10\). The composite Simpson's 1/3 Rule estimates the area under the curve between each subsequent set of 3 points (The "1/3" in the name of the rule comes from the resulting formula, not the fact that it uses 3 points).

Simpson's 1/3 Rule can be derived from a quadratic polynomial fit over 3 points and integrates the quadratic polynomial to find the area under the curve for the interval \([x_i,x_{i+2}]\). In the above example, the first area is shown as light green over the interval \([0,2.5]\), and the quadratic polynomial fit is \(y=3.25x^2-8.75x+10\). The green area is estimated as \(\int_0^{2.5}3.25x^2-8.75x+10dx=14.58333\). It is an estimate because the quadratic fit may not exactly fit the true curve over this interval, even though we can exactly integrate a quadratic polynomial.

The red and purple areas for intervals \([2.5,4.2]\) and \([4.2,6.5]\) are calculated in turn.

If there is an even number of points like this example, the final area (shown as light brown) is calculated by using a quadratic fit of the last 3 points but integrating only over the final interval \([x_{N-1},x_N]=[6.5,7]\).

The L_SIMPSON function uses the formulas provided in references [1-3] rather than using L_POLYFIT and L_POLYINT. Aside from machine precision error, the results are the same as using L_POLYFIT and L_POLYINT, but with greater efficiency. Higher-order Newton-Cotes formulas (the Trapezoidal Rule being 1st-order and Simpson's 1/3 Rule being 2nd-order) can be implemented using L_POLYFIT and L_POLYINT.

Uniformly Spaced Points

If the x values are uniformly spaced, you can omit the x parameter and specify the value for dx=Δx where \(h={\Delta}x\) is constant. With a constant h spacing, the formulas for Simpson's 1/3 Rule are simplified to:

The formula for the final area if there are an even number of points simplifies to:

Cumulative Integration

I have not implemented a cumulative option for L_SIMPSON like with L_TRAPZ or L_PPINT because the function does not produce a result for each corresponding value of y (the integration occurs over 3 points rather than just 2).

With some limited experimentation, I have found that integrating a cubic spline using L_PPINT and L_CSPLINE can produce results with approximately the same order of error as L_SIMPSON when L_PDIFF(x,y,1,3) is used to estimate the slopes for L_CSPLINE. More rigor in the derivation of the error would be needed to further test this.

Lambda Formula

This code for using L_SIMPSON in Excel is provided under the License as part of the LAMBDA Library, but to use just this function, you may copy the following code directly into your spreadsheet.

/**
* Numerical Integration using Composite Simpson's 1/3 Rule
*/
L_SIMPSON = LAMBDA(y,[x],[dx],
LET(doc,"https://www.vertex42.com/lambda/simpson.html",
    dx,IF(ISOMITTED(dx),1,dx),
    n,ROWS(y),
    integral,REDUCE(0,SEQUENCE(INT((n-1)/2),1,1,2),LAMBDA(acc,i,
        IF(ISOMITTED(x),
            acc+(1/3*dx*(INDEX(y,i)+4*INDEX(y,i+1)+INDEX(y,i+2))),
            LET(
                h_1,INDEX(x,i+1)-INDEX(x,i),
                h_2,INDEX(x,i+2)-INDEX(x,i+1),
                area,(h_1+h_2)/6*(
                    (2-h_2/h_1)*INDEX(y,i) + 
                    (h_1+h_2)^2/(h_1*h_2)*INDEX(y,i+1) + 
                    (2-h_1/h_2)*INDEX(y,i+2)),
                acc+area
            )
        )
    )),
    lastinterval,IF(ISODD(n),0,
        IF(ISOMITTED(x),
            dx*(5/12*INDEX(y,n)+2/3*INDEX(y,n-1)-1/12*(INDEX(y,n-2))),
            LET(
                h_1,INDEX(x,n)-INDEX(x,n-1),
                h_2,INDEX(x,n-1)-INDEX(x,n-2),
                alpha,(2*h_1^2+3*h_1*h_2)/(6*(h_2+h_1)),
                beta,(h_1^2+3*h_1*h_2)/(6*h_2),
                nu,(h_1^3)/(6*h_2*(h_2+h_1)),
                alpha*INDEX(y,n)+beta*INDEX(y,n-1)-nu*INDEX(y,n-2)
            )
        )
    ),
    integral+lastinterval
));

Name: L_SIMPSON
Description: Numerical Integration using Simpson's 1/3 Rule
Arguments: y, x, dx
Function:
[in the works]

L_SIMPSON Examples

Example 1

The L_PPINT function used the example where \(f(x)=x\sin(x)/4+4\) was integrated between 0 and 6π with only 31 points. The exact value is (70.6858347...).

The following is the result of using L_SIMPSON for the same 31 points as in the L_PPINT example:

=LET(
    xs, L_LINSPACE(0,6*PI(),31),
    ys, xs/4*SIN(xs)+4,
    L_SIMPSON(ys,xs)
)

Result: 70.681554
Error:  -0.004280

Note: The trapezoidal rule using L_TRAPZ for these same 31 points results in 70.8419 with an error of 0.156.