Почему мое приближение слишком велико с использованием правила составного Симпсона в R (численное интегрирование)?

Я пытаюсь аппроксимировать следующий интеграл, используя численное интегрирование в R:

Подлежащий вычислению интеграл,

где функция mu определяется этой формулой:

введите описание изображения здесь

Для этого я реализовал составное правило Симпсона как функцию в R, которая принимает в качестве параметров функцию (интегрировать), интервал интегрирования ([a, b] ) и желаемое количество подинтервалов (n).

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

Мой метод заключался в том, чтобы сначала определить внутренний интеграл в терминах его приближения Составного Симпсона как функции t в R. Затем снова использовать правило Составного Симпсона, чтобы вычислить внешний интеграл, просмотрев внутреннее приближение как интегрируемая функция.

При этом внутреннее приближение является правильным, если рассчитывается само по себе, как и ожидалось, но приближение всего выражения становится слишком большим, и я не могу понять, почему.

Я сравниваю приближения с данными Maple; внутреннее выражение, вычисленное само по себе с использованием t = 20, должно дать 0,8157191, а все выражение должно быть 12,837. R правильно вычисляет 0,8157191, но дает 32,9285 для всего выражения.

Я пробовал упростить использование множества различных математических функций и сделать функции независимыми от t в R, но, похоже, все они приводят к одной и той же ошибке. Итак, чтобы подвести итог, мой вопрос: почему неправильно аппроксимируется только внешний интеграл?

Я был бы очень признателен за любые подсказки или указатели - я включил сюда свой код, иллюстрирующий проблему:

compositesimpson <- function(integrand, a, b, n) {
  
  h<- (b-a)/n #THE DEFINITE INTERVAL IS SCALED BY
  #THE DESIRED NUMBER OF SUBINTERVALS
  
  xi<- seq.int(a, b, length.out = n+1) #DIVIDES THE DEFINITE INTERVAL INTO THE
  xi<- xi[-1]                          #DESIRED NUMBER OF SUBINTERVALS, 
  xi<- xi[-length(xi)]                 #EXCLUDING a AND b
  
  #THE APPROXIMATION ITSELF
  approks<- (h/3)*(integrand(a) + 2*sum(integrand(xi[seq.int(2, length(xi), 2)])) + 
                     4*sum(integrand(xi[seq.int(1, length(xi), 2)])) + integrand(b))
  return(approks)
  
}

# SHOULD YIELD -826.5755 BY Maple, SO THE FUNCTION IS WORKING HERE
ftest<- function(x) {
  return(exp(2*x)*sin(3*x))
}
compositesimpson(ftest, -4, 4, 100000)


# MU FUNCTION FOR TESTING
mu.01.kvinde<- function(x){ 0.000500 + 10^(5.728 + 0.038*(x+48) -10)}


#INNER INTEGRAL AS A FUNCTION OF ITS APPROXIMATION
indreintegrale.person1<- function(t){
  indre<- exp(-compositesimpson(mu.01.kvinde, 0, t, 100000))
  return(indre)
}

indreintegrale.person1(20) #YIELDS 0.8157191, WHICH IS CORRECT

compositesimpson(indreintegrale.person1, 20, 72, 100000) #YIELDS 32.9285,
#BUT SHOULD BE 12.837 ACCORDING TO MAPLE

r numerical-methods numerical-integration simpsons-rule
person VHS    schedule 21.02.2021    source источник

Ответы (1)


arrow_upward
0
arrow_downward

Это как-то связано с попыткой использовать векторизацию на двух уровнях рекурсии, и она не делает того, что вы хотите. Например. сравнивать

indreintegrale.person1(20) 
#> [1] 0.8157191
indreintegrale.person1(c(20, 72))
#> [1] 0.8157191 0.4801160
indreintegrale.person1(72) 
#> [1] 2.336346e-10

Я думаю, что средний ответ неверен, но два других верны.

Самое быстрое исправление, сделайте эту замену:

indreintegrale.person1 <- function(t){
  sapply(t, function(t2) exp(-compositesimpson(mu.01.kvinde, 0, t2, 100000)))
}

и теперь он дает ожидаемый ответ (но для расчета требуется немного больше времени!).

person pseudospin    schedule 21.02.2021
comment
Спасибо за ваше объяснение и решение! Если бы я использовал это для вычисления других интегралов того же типа, вы бы порекомендовали просто продолжать использовать sapply для внутренних выражений, или вы бы внесли изменения в compositesimpson Сама функция? - person VHS; 23.02.2021
comment
Я думаю, вы должны относиться к функции compositesimpson() как к невекторизованной (вызовы seq в ней затруднили бы ее векторизацию). Однако любая функция, которую вы добавляете как integrand, должна быть векторизована, чтобы compositiesimpson работала. sapply - это один из способов превратить невекторизованную функцию в векторизованную. - person pseudospin; 23.02.2021