柯西主值數值積分

我使用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，使得主值的一個更好的辦法，比切斷參數玩耍和手動調節它，直到它給出了一個結果（不一定準確）？

來源

2013-03-06 baptiste

是否有可能重新說明問題 - 如果您將評估函數的值分配爲1/fun，那麼您最終會在接近奇點的多個點進行評估。另一種觀察方式是採用反函數並將其整合兩次 - 在奇點的任一側。現在你的x會擴展到無窮大，但你的奇點已經消失... – Floris 2013-03-06 23:16:22

@baptiste你可以放寬公差參數，因爲精度不是必須的。例如用'rel.tol = 1e-6'來收斂。 – agstudy 2013-03-06 23:42:20

@弗洛裏你建議改變變量，對吧？我不知道這種類型的被積函數是好的，但我想這是一種選擇。 – baptiste 2013-03-07 04:39:27

如果考慮整體上集中於單一，可以使用變量的變化的對稱化積的時間間隔：奇點消失，在紙張Numerical evaluation of a generalized Cauchy principal value (A. Nyiri and L. Bayanyi, 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

來源

2013-03-07 17:48:09

看起來真的很好，謝謝！我會隨着我的積分一起去，然後再回到你身邊。 – baptiste 2013-03-07 19:09:49

柯西主值數值積分

回答

相關問題