Я хочу выполнить численное интегрирование в одном измерении, где подынтегральная функция имеет векторное значение. integrate()
допускает только скалярные интегранты, поэтому мне нужно будет вызывать его несколько раз. Пакет cubature
кажется хорошо подходящим, но он плохо работает с одномерными интегралами. Рассмотрим следующий пример (интегральное выражение со скалярным знаком и одномерное интегрирование):
library(cubature)
integrand <- function(x, a=0.01) exp(-x^2/a^2)*cos(x)
Nmax <- 1e3
tolerance <- 1e-4
# using cubature's adaptIntegrate
time1 <- system.time(replicate(1e3, {
a <<- adaptIntegrate(integrand, -1, 1, tol=tolerance, fDim=1, maxEval=Nmax)
}) )
# using integrate
time2 <- system.time(replicate(1e3, {
b <<- integrate(integrand, -1, 1, rel.tol=tolerance, subdivisions=Nmax)
}) )
time1
user system elapsed
2.398 0.004 2.403
time2
user system elapsed
0.204 0.004 0.208
a$integral
> [1] 0.0177241
b$value
> [1] 0.0177241
a$functionEvaluations
> [1] 345
b$subdivisions
> [1] 10
Каким-то образом adaptIntegrate
, похоже, использует гораздо больше вычислений функций для аналогичной точности. Оба метода, по-видимому, используют квадратуру Гаусса-Кронрода (одномерный случай: квадратурное правило Гаусса с 15 точками), хотя ?integrate
добавляет «алгоритм Эпсилон Винна». Объясняет ли это большую разницу во времени?
Я открыт для предложений альтернативных способов работы с векторными интегралами, такими как
integrand <- function(x, a = 0.01) c(exp(-x^2/a^2), cos(x))
adaptIntegrate(integrand, -1, 1, tol=tolerance, fDim=2, maxEval=Nmax)
$integral
[1] 0.01772454 1.68294197
$error
[1] 2.034608e-08 1.868441e-14
$functionEvaluations
[1] 345
Спасибо.
fDim=2
(последний пример, тоже 345 оценок), сравнение - это просто случай двойного вызоваintegrate
,str(lapply(c(integrand1, integrand2), integrate, -1,1, rel.tol=tolerance, subdivisions=Nmax))
дает 10+1 = 11 оценок. Я хочу сказать, что да,adaptIntegrate
нацелен на многомерную интеграцию и, возможно, на векторные интегралы, но случай одномерной интеграции гораздо менее эффективен, чем повторный вызовintegrate
, но с большим запасом (здесь ~ 30 раз). - person baptiste   schedule 17.03.2013