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?

  • $\begingroup$ So, SimpsonIntegral[2, {x, 0, 1}, 4] is not expected to evaluate? $\endgroup$ Commented Sep 20, 2020 at 4:50
  • $\begingroup$ It is not clear what role play a, b in your function? $\endgroup$ Commented Sep 20, 2020 at 11:46

2 Answers 2


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, 

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


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}}
  • $\begingroup$ 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)|$. $\endgroup$ Commented Sep 21, 2020 at 23:50
  • $\begingroup$ Do you mean formula or numerical result? For numerical result I gave some update. $\endgroup$ Commented Sep 22, 2020 at 10:31
  • $\begingroup$ 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. $\endgroup$ Commented Sep 24, 2020 at 8:49

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.

  • $\begingroup$ Thanks, I believe this should work! $\endgroup$
    – Carlos M.
    Commented Sep 19, 2020 at 23:44
  • $\begingroup$ @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. $\endgroup$
    – m_goldberg
    Commented Sep 20, 2020 at 6:36
  • $\begingroup$ @m_goldberg I am clicking (+1) for you solution but the definition of SimpsonIntegral[] should be reconsidered. $\endgroup$ Commented Sep 20, 2020 at 11:54
  • $\begingroup$ @AlexTrounev. I am in full agreement with you. $\endgroup$
    – m_goldberg
    Commented Sep 20, 2020 at 21:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.