Definite integral

This small program consists of a function that can compute the definite integral of a function over a given interval. Full description with examples can be found here and the code is here

integral_definida: approximates the definite integral of a function over [a,b], given a number of partitions.

• Arguments:

  • funcion: Mathematical function that we want to integrate.

  • a: Lower limit of integration.

  • b: Upper limit of integration.

  • n: Number of partitions

  • aprox: Chosen integration technique . The definite integral can be approximated by rectangles (‘rect’), the trapezoid rule (‘trap’) or by Simpson’s rule (‘simp’). The latter gives the best approximations, but when n is an even number.

    Approximation by rectangles: $\small\displaystyle \int_a^b f(x)dx = \lim_{n\rightarrow \infty} I_n=\lim_{n\rightarrow \infty}S_n$.

    $\small In=\sum_{i=1}^n f(m_i)\cdot (\frac{b-a}{n})$ ; $\small Sn=\sum_{i=1}^n f(M_i)\cdot (\frac{b-a}{n})$

    Trapezoid rule: $\small\displaystyle \int_a^b f(x) dx \approx \frac{b-a}{2n}[f(x_0) + 2f(x_1) + 2f(x_2) + 2f(x_3) + … + 2f(x_{n-2}) + 2f(x_{n-1}) + f(x_n)]$

    Simpson’s rule: $\small\displaystyle \int_a^b f(x) dx \approx \frac{b-a}{3n}[f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + … + 2f(x_{n-2}) + 4f(x_{n-1}) + f(x_n)]$

• Output:

  • If aprox='rect’: list with two elements representing the lower sum (In) and the upper sum (Sn).

  • If aprox='trap’ o aprox='simp’: vector of length one representing the estimate of the definite integral over a and b.

integral_definida <- function(funcion,a,b,n,aprox='rect'){

  interv <- (b-a)/n
  
  seg <- seq(from=a,to=b,by=interv)
  
  if (aprox=='trap'){
    
    strt_end <- c(funcion(seg[1]),funcion(seg[n+1]))
    mid <- 2*funcion(seg[2:n])
    
    return(( (b-a)/(2*n) )*sum(strt_end,mid))
  }
  else if (aprox=='simp'){
    
    strt_end <- c(funcion(seg[1]),funcion(seg[n+1]))
    seg_cut <- seg[2:n]
    
    indices <- 1:length(seg_cut)
    mid <- c(4*funcion(seg_cut[indices%%2 != 0]),
             2*funcion(seg_cut[indices%%2 == 0]))
    
    return(( (b-a)/(3*n) )*sum(strt_end,mid))
  }
  else{
    
    In <- double(1)
    Sn <- double(1)
  
    for (i in 1:n){
      fmi <- min(funcion(seg[i]:seg[i+1]))
      In <- In+fmi*interv
      fMi <- max(funcion(seg[i]:seg[i+1]))
      Sn <- Sn+fMi*interv
    }
    return(list(In=In,Sn=Sn))
  }
}