台形公式で積分

数値計算を勉強したことがないのがバレバレですねwww。

my.integrate <- function(f,a,b,n){
  ak <- seq(from=a,to=b,length.out=n)
  sum <- 0
  for(i in 2:n){
    sum <- sum + (ak[i] - ak[i-1]) * (f(ak[i-1]) + f(ak[i])) / 2
  }
  return(sum)
}

よしよし。

my.integrate(f=function(x){x^2},a=-1,b=1,n=10000)
[1] 0.6666667