Вычисление определенного интеграла в R

integrate() функция вычисляет определенный интеграл. 
синтаксис: integrate(f, lower, upper)

Параметры:

f:  интегрируемая функция

lower: нижний предел

upper: верхний предел

# Create function 
f<-function(x) sin(x)
 # Integrate from 0 to 1 
integrate(f, 0, pi) 

> # Create function
> f<-function(x) sin(x)
> # Integrate from 0 to 1
> integrate(f, 0, pi)
2 with absolute error < 2.2e-14

# Create function
f<-function(x) 1/(1+x^2)
# Integrate from 0 to 1
integrate(f, 0, 1)
1-(0.78539816339744839002*4-pi)

> # Create function
> f<-function(x) 1/(1+x^2)
> # Integrate from 0 to 1
> integrate(f, 0, 1)
0.78539816339744839002 with absolute error < 8.7e-15
> 1-(0.78539816339744839002*4-pi)
[1] 0.99999999999999955591

Так же возможно установить дополнительный пакет и вычислять определенный интеграл методом Симпсона с заданной точностью. 

#install.packages("fda.usc")
library(fda.usc)
options(digits = 20)
x<-seq(0,pi,length = 20)
fx<-fdata(sin(x),x)
int.simpson(fx,method ="ESR" )

# Posible values for method are:
# "TRAPZ": Trapezoid rule integration.
# "CSR": Composite Simpson's rule integration.
# "ESR": Extended Simpson's rule integration.
# If method=NULL (default), the value of par.fda.usc$int.method is used.
# "TRAPZ": Интеграция правил трапеции.
# "CSR": Комплексная интеграция правил Симпсона.
# "ESR": Расширенная интеграция правил Симпсона.

x = seq(0, 1, length = 10000)
fx = fdata(1/(1+x^2),x)
i = int.simpson(fx,method ="ESR" )
1 - abs(i*4-pi)
i

> x = seq(0, 1, length = 10000)
> fx = fdata(1/(1+x^2),x)
> i = int.simpson(fx,method ="ESR" )
> 1 - abs(i*4-pi)
[1] 1
> i
[1] 0.785398163397448279

Суть метода парабол сводится к более сложной фор муле значение функции на промежутку. Симпсон вывел приближённую формулу вычисление интеграла на промежутке из предположения, что функция на промежутке является параболой. Использование этой функции даёт лучший результат в сравнении с методом прямоугольника и трапеции. Формула Симпсона представлена на рисунке

# Вычисление интеграла 
# методом парабол (Симпсона)
options(digits=18) # Установка максимального количества используемых цифр
fx = function(x) 1/(1+x^2)
x1 = 0
x2 = 1
n = 100000
#
dx = integer = (x2 - x1)/n
for (i in  1:n) {
  integer = integer + (dx/6)*(fx(x1) + 4*fx((2*x1+dx)/2)+
                                fx(x1+dx))
  x1 = x1 + dx
}
integer

Хотя вычислить интеграл точно по формуле Симпсона возможно, воспользовавшись конструкцией приведенной ниже.

# Вычисление интеграла 
# методом парабол (Симпсона) по точной формуле

options(digits=20)

fx = function(x) 1/(1+x^2)
a = 0
b = 1
n = 100000
i = 0
stepx = (b-a)/(n*3)
x = seq(a, b, (b-a)/n)
y = fx(x)
i = i + y[1] + y[length(y)]
for (k in 2:(length(y)-1)) {
  if (k%%2==0) i = i + y[k]*4
  if (k%%2!=0) i = i + y[k]*2
  
}
I = stepx*i
I

> I
[1] 0.7853981633974453924196