How can I fix my Simpson's rule integration code?

Question

I was tasked with writing a function on Mathematica similar to NIntegrate, but that uses the Composite Simpson's Rule as its method for doing the calculation. The function should take as arguments, the function to be integrated (aka the integrand), the integration variable (as well as its upper and lower limits) and how many equally spaced subdivisions it should use to do the calculation.

The function is only required to integrate in regards to one variable; however, if presented with a multivariable function it must be able to integrate over only one designed variable and ignore the others.

I was able to handle that following another question here. And came up with the following code:

    SetAttributes[SimpsonIntegral, HoldAll]
    SimpsonIntegral[f_[a___, var_, b___], {var_, xmin_, xmax_}, steps_] :=
    (xmax - xmin)/(3 steps) Sum[f[a, xmin + (xmax - xmin)/steps (2*y - 2), b] + 
    4*f[a, xmin + (xmax - xmin)/steps (2*y - 1), b] + 
    f[a, xmin + (xmax - xmin)/steps (2*y), b], {y, 1, steps/2}];

It worked perfectly fine with single variable functions (like Sin[x]) and even multivariable functions (like BesselJ[n,z]). However, it fails when the input function takes an expression as its argument. For example:

SimpsonIntegral[Sin[x - 1], {x, 0, 1}, 6]

would simply return unevaluated:

SimpsonIntegral[Sin[x - 1], {x, 0, 1}, 6]

instead of the various terms of the sum.

How can I fix this?

Answer 1

I can recommend more practical for numerical calculations next definition:

SetAttributes[SimpsonIntegral, HoldAll]
SimpsonIntegral[f_, x_, xmin_, xmax_, 
   steps_] := (xmax - xmin)/(3 steps) (Sum[
      f /. {x -> xmin + (xmax - xmin)/steps (2*y - 2)}, {y, 1, 
       steps/2}] + 
     4*Sum[f /. {x -> xmin + (xmax - xmin)/steps (2*y - 1)}, {y, 1, 
        steps/2}] + 
     Sum[f /. {x -> xmin + (xmax - xmin)/steps (2*y)}, {y, 1, 
       steps/2}]);

Then we have for BesselJ[]

SimpsonIntegral[BesselJ[2, 3 x + 2], x, 0, 1, 8] // N

Out[]= 0.366069 

Let compare with NIntegrate

NIntegrate[BesselJ[2, 3 x + 2], {x, 0, 1}]

Out[]= 0.3660498384281397

As expected the error for unit interval is bounded as $\frac{h^4}{180}max|f^{(4)}(x)|$, so with h=1/8 it gives $1.35634\times 10^{-6}\times 20.496=2.78 \times 10^{-5} $, and we have from above error of $1.9\times 10^{-5}$. From the other side we can calculate exactly

i = Integrate[BesselJ[2, 3 x + 2], {x, 0, 1}]

Out[]= 1/72 (125 HypergeometricPFQ[{3/2}, {5/2, 3}, -(25/4)] - 
   8 HypergeometricPFQ[{3/2}, {5/2, 3}, -1]) 

Then we can numerically compute

i // N

0.36604983842813943

Let compare results given by NIntegrate and Integrate to support error bound verification. So in this case NIntegrate really working with MachinePrecision. To calculate $max|f^{(4)}(x)|$ we use

FindMaximum[D[BesselJ[2, 3 x + 1], {x, 4}], {x, .8}]

Out[]= {20.496, {x -> 0.752259}}

Answer 2

The function MyPlot, which is the model for your implementation, is specialized to a certain argument pattern and won't evaluate when given Sin[x + 1] as its argument, either. It requires the variable of interest to appear in its argument sequence "naked". It can not be a factor in an expression. The simple solution is to define a new function which provides the required isolation. Like so:

f[x_] := Sin[x - 1]
SimpsonIntegral[f[x], {x, 0, 1}, 6]

1/18 (-4 Sin[1/6] - 2 Sin[1/3] - 4 Sin[1/2] - 2 Sin[2/3] - 4 Sin[5/6] - Sin[1])

If that isn't an acceptable solution, then you need find a different approach.