The easiest numerical integral with Simpson's rule

The easiest numerical integral would be one with Simpson's rule. Let's follow the rule and don's stop on tracks!

Composite Simpson's rule

The formula of Simpson's rule is: \begin{align} \int_a^b f(x) dx \approx \frac{1}{3} \left[ f(x_0) + 2 \sum_{j=1}^{n/2-1} f(x_{2j}) + 4 \sum_{j=1}^{n/2} f(x_{2j-1}) + f(x_n) \right]. \nonumber \end{align} Consider the function \(f(x) = 4x + 3x^2+ 2x^3+ 1x^4\),

f = function(x) {4*x + 3*x^2 + 2*x^3 + x^4}
curve(f,2,10)

and its integral \begin{align} \int_2^{10} f(x) dx &= \int_2^{10} 4x + 3x^2+ 2x^3+ 1x^4 dx \nonumber \\ &= \left[ 2 x^2 + x^3 + \frac{1}{2} x^4 + \frac{1}{5} x^5 \right]_2^{10}. \nonumber \end{align}

Analytical solution

The analytical solution is

F = function(x){2*x^2+x^3+1/2*x^4+1/5*x^5}
F(10)-F(2)

\begin{align} \int_2^{10} f(x) dx =26169.6. \nonumber \end{align}

Approximated solution by composite Simpson's rule

The numbers of numerical integral obtained by composite Simpson's rule are:

simp = function(f,n,a,b){
    h = (b-a)/n
    x = seq(a,b,h)
    y = f(x)
    coef = c(1,rep(c(4,2),(n/2-1)),4,1)
    res = sum(h/3*y*coef)
    return(res)
}

c(
  simp(f,10,2,10),
  simp(f,50,2,10),
  simp(f,100,2,10),
  simp(f,1000,2,10)
)

\begin{align} \mbox{simp}(f,10,2,10)&= 26170.0369066667 \nonumber \\ \mbox{simp}(f,50,2,10)&= 26169.6006990507 \nonumber \\ \mbox{simp}(f,100,2,10)&= 26169.6000436907 \nonumber \\ \mbox{simp}(f,1000,2,10)&= 26169.6000000044 \nonumber \end{align} Obviously, it looks not efficient, however it would be useful for development stages, because the possibility for making mistake in the code is quite low, hopefully.