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$ – J. M.'s discontentment Sep 20 at 4:50
  • $\begingroup$ It is not clear what role play a, b in your function? $\endgroup$ – Alex Trounev Sep 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.

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

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}}
| improve this answer | |
  • $\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$ – A little mouse on the pampas Sep 21 at 23:50
  • $\begingroup$ Do you mean formula or numerical result? For numerical result I gave some update. $\endgroup$ – Alex Trounev Sep 22 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$ – A little mouse on the pampas Sep 24 at 8:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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