Численное интегрирование главного значения Коши

Я использую функцию интегрирования R, но подынтегральное выражение доставляет мне некоторые проблемы (у него есть особенность). Учти это,

distrib <- function(x){
  dnorm(x, mean=700,sd=50)
}

phi <- function(nu, epsilon=20, maxi=Inf){

  fun <- function(nup)  {
    lambdap <- 1e7/nup
    distrib(lambdap) / (nup*(nup - nu))
  }
  # first part: from epsilon to nu - epsilon
  I1 <- integrate(fun, epsilon, nu-epsilon, rel.tol = 1e-8,
                  subdivisions = 200L)
  # second part: from nu + epsilon to Infty
  I2 <- integrate(fun, nu+epsilon, maxi, rel.tol = 1e-8, 
                  subdivisions = 200L)

  I1$value + I2$value 

}


x <- seq(1e7/500, 1e7/1000, length=200)

phi1 <- sapply(x, phi)    
phi1 <- sapply(x, phi, epsilon=5)
#Error in integrate(fun, epsilon, nu - epsilon, rel.tol = 1e-08, subdivisions = #200L) : 
#  the integral is probably divergent

Есть ли лучший способ вычисления таких основных значений в R, чем играть с параметром отсечки и настраивать его вручную, пока он не даст результат (не обязательно точный)?


r numerical-integration
person baptiste    schedule 06.03.2013    source источник
comment
Можно ли переформулировать проблему: если вы распределяете значения, при которых вы оцениваете функцию, как 1/fun, то вы в конечном итоге оцениваете во многих точках, близких к сингулярности. Другой способ взглянуть на это — взять обратную функцию и интегрировать ее дважды — по обе стороны от сингулярности. Теперь ваш x простирается до бесконечности, но ваша сингулярность исчезла...   -  person Floris    schedule 07.03.2013
comment
@baptiste, вы можете ослабить параметр допуска, поскольку точность не требуется. например, с rel.tol = 1e-6 он сходится.   -  person agstudy    schedule 07.03.2013
comment
@Floris, вы предлагаете изменить переменные, верно? Я не знаю подходящего для этого типа подынтегрального выражения, но я думаю, что это один из вариантов.   -  person baptiste    schedule 07.03.2013
comment
@agstudy да, я мог бы настроить параметры интеграции, но я пытаюсь посмотреть, придумали ли люди более надежные и общие методы для работы с такими неправильными интегралами.   -  person baptiste    schedule 07.03.2013
comment
Есть ли у вашего подынтегрального выражения две особенности — в 0 и nu? Я не уверен, что правильно читаю ваш код... /(nup*(nup-nu))   -  person Floris    schedule 07.03.2013
comment
ага, 0 и ню, извините непонятно.   -  person baptiste    schedule 07.03.2013

Ответы (1)


arrow_upward
1
arrow_downward

Если вы рассматриваете интеграл на интервале с центром в сингулярности, вы можете использовать замену переменной для симметризации подынтегральной функции: сингулярность исчезает, как подробно описано в статье Численная оценка обобщенного главного значения Коши (А. Ньири и Л. Баяни, 1999 г.).

# Principal value of \( \int_{x_0-a}^{x_0+a} \dfrac{ f(x) }{ h(x) - h(x0) } dx \)
pv_integrate <- function( f, h, x0, a, epsilon = 1e-6 ) {
  # Estimate the derivatives of f and h
  f1 <- ( f(x0+epsilon) - f(x0-epsilon) ) / ( 2 * epsilon )
  h1 <- ( h(x0+epsilon) - h(x0-epsilon) ) / ( 2 * epsilon )
  h2 <- ( h(x0+epsilon) + h(x0-epsilon) - 2 * h(x0) ) / epsilon^2 
  # Function to integrate. 
  # We just "symmetrize" the integrand, to get rid of the singularity,
  # but, to be able to integrate it numerically, 
  # we have to provide a value at the singularity.
  # Details: http://mat76.mat.uni-miskolc.hu/~mnotes/downloader.php?article=mmn_16.pdf
  g <- function(u) 
    ifelse( u != 0,
      f( x0 - u ) / ( h(x0 - u) - h(x0) ) + f( x0 + u ) / ( h(x0+u) - h(x0) ),
      2 * f1 / h1 - f(x0) * h2 / h1^2
    )
  integrate( g, 0, a )
}

# Compare with the numeric results in the article
pv_integrate( function(x) 1,      function(x) x,    0, 1  ) # 0
pv_integrate( function(x) 1,      function(x) x^3,  1, .5 ) # -0.3425633
pv_integrate( function(x) exp(x), function(x) x,    0, 1  ) # 2.114502
pv_integrate( function(x) x^2,    function(x) x^4,  1, .5 ) # 0.1318667
person Vincent Zoonekynd    schedule 07.03.2013
comment
выглядит очень хорошо, спасибо! Я попробую со своим подынтегральным выражением и вернусь к вам позже. - person baptiste; 07.03.2013