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

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?

• So, SimpsonIntegral[2, {x, 0, 1}, 4] is not expected to evaluate? – J. M.'s ennui Sep 20 '20 at 4:50
• It is not clear what role play a, b in your function? – Alex Trounev Sep 20 '20 at 11:46

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.

• Thanks, I believe this should work! – Carlos M. Sep 19 '20 at 23:44
• @CarlosM. I'm glad you find my answer useful. Please consider accepting it. You can do that by clicking on the check mark that appears on the left of the answer below the down arrow. – m_goldberg Sep 20 '20 at 6:36
• @m_goldberg I am clicking (+1) for you solution but the definition of SimpsonIntegral[] should be reconsidered. – Alex Trounev Sep 20 '20 at 11:54
• @AlexTrounev. I am in full agreement with you. – m_goldberg Sep 20 '20 at 21:07

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}}

• This is a great answer. I wonder if you can use MMA to get the truncation error $\frac{h^4}{180}max|f^{(4)}(x)|$. – A little mouse on the pampas Sep 21 '20 at 23:50
• Do you mean formula or numerical result? For numerical result I gave some update. – Alex Trounev Sep 22 '20 at 10:31
• If you can, I hope you can write a function as general as possible to solve the representation of truncation error in numerical analysis. Thank you. – A little mouse on the pampas Sep 24 '20 at 8:49