integracja numeryczna z możliwym dzieleniem przez „zero”

Próbuję zintegrować

\int_{0}^{1} t^{2 n + 2} \exp (\frac{α r_{0}}{t}) d t

$\int^1_0 t^{2n+2}\exp\left({\frac{\alpha r_0}{t}}\right)dt$

która jest prostą transformacją

\int_{1}^{\infty} x^{2 n} \exp (- α r_{0} x) d x

$\int^{\infty}_1 x^{2n}\exp(-\alpha r_0 x)dx$

za pomocą $t = \frac1{x}$ ponieważ trudno jest liczbowo aproksymować nieprawidłowe całki. Prowadzi to jednak do problemu oceny nowego całka blisko zera. Bardzo łatwo będzie uzyskać odpowiednią liczbę węzłów kwadraturowych, ponieważ przedział ma tylko długość 1 (więc porównywalny $dt$ może być bardzo mały), ale jakie uwagi należy wziąć pod uwagę przy integracji w pobliżu zera?

Na pewnym poziomie myślę, że to po prostu biorąc $\int^1_\epsilon t^{2n+2}\exp({\frac{\alpha r_0} {t}})dt$ jest dobrym pomysłem gdzie $\epsilon$ to jakaś niewielka liczba. Jednak jaki numer powinienem wybrać? Czy powinien to być epsilon maszynowy? Czy podział przez epsilon maszynowy jest liczbą dobrze skwantyfikowaną? Co więcej, jeśli podział mojej maszyny epsilon (lub blisko niej) daje niewiarygodnie dużą liczbę, to biorąc $\exp(\frac{1}{\epsilon})$ stanie się jeszcze większy.

Jak mam to uwzględnić? Czy istnieje sposób na dobrze zdefiniowaną całkę numeryczną tej funkcji? Jeśli nie, jaki jest najlepszy sposób integracji funkcji?

quadrature accuracy

— drjrm3
źródło

Czy zastanawiałeś się nad użyciem Monte Carlo?

— Faheem Mitha

Mam wrażenie, że nie rozwiązałoby to problemu. Integracja Monte Carlo jest często zarezerwowana dla całek wysokowymiarowych. Miałem takie same problemy z Monte Carlo, po prostu miałem mniej kontroli nad tym, gdzie moja funkcja jest oceniana.

— drjrm3

Możesz mieć rację.

— Faheem Mitha

Myślę, że dobrze byłoby mieć odpowiedź (być może na osobne, bardziej ogólne pytanie) wyjaśniającą, w jaki sposób dokonuje się integracji numerycznej, gdy funkcja jest rozbieżna na jednym limicie, w ogólnym przypadku, w którym nie jest możliwe wykonanie całki analitycznej. Z drugiej strony, równie dobrze można to znaleźć w przepisach numerycznych ...

— David Z

@Faheem: „Monte Carlo to bardzo zła metoda; należy jej używać tylko wtedy, gdy wszystkie metody alternatywne są gorsze”. - Alan Sokal

— JM

Odpowiedzi:

Można tego dokonać poprzez integrację części:

\int_{1}^{\infty} x e^{- a x} = \frac{- 1}{a} x e^{- a x} ∣_{1}^{\infty} - \frac{- 1}{a} \int_{1}^{\infty} e^{- a x} = \frac{e^{- a}}{a} + \frac{e^{- a}}{a^{2}} = \frac{a + 1}{a^{2}} e^{- a}

$\int^\infty_1 x e^{-ax} = \frac{-1}{a} x e^{-ax}\mid^\infty_1 - \frac{-1}{a} \int^\infty_1 e^{-ax} = \frac{e^{-a}}{a} + \frac{e^{-a}}{a^2} = \frac{a+1}{a^2} e^{-a}$ i kontynuowanie przez indukcję

\int_{1}^{\infty} x^{k} e^{- a x} = \frac{- 1}{a} x^{k} e^{- a x} ∣_{1}^{\infty} - \frac{- k}{a} \int_{1}^{\infty} x^{k - 1} e^{- a x} = \frac{e^{- a}}{a} + \frac{k}{a} \int_{1}^{\infty} x^{k - 1} e^{- a x}

$\int^\infty_1 x^k e^{-ax} = \frac{-1}{a} x^k e^{-ax}\mid^\infty_1 - \frac{-k}{a} \int^\infty_1 x^{k-1} e^{-ax} = \frac{e^{-a}}{a} + \frac{k}{a} \int^\infty_1 x^{k-1} e^{-ax}$ po to aby

I (k) = \frac{e^{- a}}{a} + \frac{k}{a} I (k - 1)

$I(k) = \frac{e^{-a}}{a} + \frac{k}{a} I(k-1)$ i

I (0) = \frac{e^{- a}}{a}

$I(0) = \frac{e^{-a}}{a}$ .

— Matt Knepley
źródło

absolutnie nie mam pojęcia, jak to przeoczyłem. Dziękuję Ci.

— drjrm3

Sprytne zamiany i integracja przez części zawsze powinny być jedną z pierwszych rzeczy, które robisz z niesfornymi całkami.

— JM,

Często dobrym pomysłem jest zapytać system algebry komputerowej, gdy masz taką całkę. Klon ocenia „

\int_{1}^{\infty} x^{2 n} \exp (- α x) d x assuming n :: nonnegint, α > 0

$\int_1^\infty x^{2n} \exp(-\alpha x) \mathrm{d}x \text{ assuming }n :: \text{nonnegint}, \alpha > 0$ „natychmiast w

Γ (2 n + 1, α) α^{- 2 n - 1}

$\Gamma(2n+1, \alpha) \alpha^{-2n-1}$ ; Jestem pewien, że Mathematica robi to samo. (Nadal dobrym pomysłem jest to weryfikacja numeryczna, oczywiście, co ci faceci zazwyczaj mogą również zrobić.)

— Erik P.

W rzeczywistości Mathematica wybiera reprezentację odpowiedzi jako ExpIntegralE [-2 n, ar]. Jeśli uruchomisz na niej FunctionExpand, to daje taką samą odpowiedź jak Maple.

— Searke

Spójrz na QUADPACK . Ma procedury integracji z (pół-) nieskończonymi domenami.

— GertVdE
źródło