Bir sayısal analiz uygulaması olarak Gauss integralinin usulüne uygun kestirilmesinde integrasyon aralığının tespiti

Matematiksel istatistik ve kuantum harmonik sarkaç başta olmak üzere, fen bilimleri ve mühendisliğin değişik sahalarında çan eğrisi fonksiyonunun integrali mutlaka karşımıza çıkar. \begin{equation*} I := \int\limits_{-\infty}^{\infty} \exp (-x^{2}) dx = \sqrt{\pi} \approx 1,77245385091 \end{equation*} Tek değişkenli bir fonksiyonun integrali olmasına rağmen, ne yazık ki bu integral tek değişkenli analizin imkanları çerçevesinde hesaplanılamaz; çok değişkenli analize başvurmak gerekir. Bu postada öyle yapmayacağız ve usulüne uygun bir biçimde $I$ integralini sayısal yöntemlerle kestirmeye çalışacağız. Usulüne uygun olmaktan kastımız, yaptığımız kestirme için bir hata sınırı tayin etmekten veyahut verilen bir hata sınırına uyan bir kestirmenin asgari sayıda hesaplama ile yapılmasından ibarettir. Ayrıca işlem sayısında da cimri davranıp, yaptığımız herbir çarpma ve toplama işlemini sayarak bunları da azaltmanın yollarını sorgulayacağız.

Hesaplamaya başlamadan önce ilk olarak çan eğrisi fonksiyonunun paritesinin çift olduğunu gözlüyoruz. Diğer bir deyişle $f(x) := e^{-x^{2}}$ için $f(-x)=f(x)$. O zaman $I$ integralini hesaplamak için çan eğrisini tüm gerçel sayılarda hesaplamak yerine sadece negatif olmayan gerçel sayılarda hesaplamamız yeterlidir. Bu gözlem yapmamız gereken işlem sayısını hemen yarı yarıya düşürür. \begin{equation*} I = 2 \int \limits_{^0}^{\infty} \exp(-x^{2}) dx \end{equation*}

İntegrali hesaplarken Riemann'ın fikrinden faydalanacak ve fonksiyonun altındaki alanı küçük dikdörtgenlere bölüp, bu dikdörtgenlerin toplamıyla kestireceğiz. Ancak burada bir sorun var: integralimizde üst limit $\infty$! Hiç kimse sonsuza kadar olan bir bölgeyi küçük dörtgenlere bölüp de burada sonsuz tane hesaplama yapmadı. Bu yüzden $I$ integralini kestirirken, integralin üst limitini $1 < a < \infty$ ile değiştireceğiz. Kuşkusuz buradan bize bir hata ($\varepsilon > 0$) gelecek.

Farkında mısınız, bir önceki paragrafta calculus'taki ve genel olarak analizdeki limit süreçlerinin tersinden gittik: önce $a$ değerini verip daha sonra bunun doğurduğu hatayı hesaplamayı önerdik. Halbuki analizde, önce hata ($\varepsilon > 0$) verilir, daha sonra bu hataya uyan $a=a(\varepsilon)$ değeri bulunur.

İntegrali $(0,a)$ aralığında hesaplayacağımıza göre, hesaplamadığımız kısım bizim "hata" dediğimiz nicelik olacaktır. Yani \begin{equation*} \varepsilon := 2 \int \limits _{a}^{\infty} \exp(-x^{2}) dx \end{equation*} olur. Ne yazık ki bu integral analitik yöntemlerle hesaplanamıyor. Hatta bu integral bir sabite kadar literatürde ${\rm erfc}(x)$ yani bütünleyen hata fonksiyonu (complementary error function) olarak bilinir ve istatistikte merkezi bir öneme sahiptir. Verdiğimiz hataya üst sınırlar bulmak için birkaç eşitsizlikten faydalanacağız. Herbiri hesaplama maliyetimizi bir miktar düşürecek.

Lemma 1: $x \geq 1$ ise, o zaman $\exp (-x^{2}) \leq \exp (-x)$.
İspat: $x \in \mathbb{R}$ için üstel fonksiyonu artandır. (Basitçe birinci türevinin her zaman pozitif olduğunu gözleyiniz.) Ayrıca $x \geq 1$ için $x^{2} \geq x$ diğer bir ifadeyle $-x^{2} \leq -x$ olur. Üstel fonksiyonun artanlığı kullanılarak ispat tamamlanır.

Lemma 1'deki eşitsizliğin her iki tarafında da intagral alırsak, $a \geq 1$ için \begin{equation*} \varepsilon = 2 \int \limits _{a}^{\infty} \exp(-x^{2}) dx \leq 2 \int \limits _{a}^{\infty} \exp(-x) dx = 2 e^{-a} \end{equation*} sınırını tayin ederiz. Diyelim ki C programlama dilinin float veri yapısına uygun olarak $\varepsilon \leq 10^{-8}$ olsun istedik. O zaman $2e^{-a} = 10^{-8}$ denklemi uyarınca $a = -\log 0,5 \times 10^{-8} = 19,11382792$ buluruz. Bu sonuç iyi bir integrasyon algoritmasıyla, kabaca $(0,20)$ aralığında çan eğrisinin integralini alırsak, istediğimiz hata çerçevesinde Gauss integralini kestirebileceğimizi söylüyor.

Lemma 2: $x \geq 1$ ise, o zaman $\exp (-x^{2}) \leq x \exp (-x^{2})$.
İspat: Bariz.

Şimdi, bu son derece aşikar eşitsizlikte her iki tarafında integralini alırsak, $a \geq 1$ için \begin{equation*} \varepsilon = 2 \int \limits _{a}^{\infty} \exp(-x^{2}) dx \leq \int \limits _{a}^{\infty} 2 x \exp(-x^{2}) dx = -\int \limits _{a}^{\infty} d \exp(-x^{2}) = e^{-a^{2}} \end{equation*} sonucunu elde ederiz. Tek değişkenli analizin imkanlarıyla bir tam diferansiyel haline getirilemeyen çan eğrisi fonksiyonunu tam diferansiyel haline getirecek şekilde $x$ ile çarptığımızı ama buna mukabil eşitlik yerine eşitsizlikle uğraşmak zorunda kaldığımızı gözleyiniz. İntegralin üst limitini sonsuz yerine $a \geq 1$ ile değiştirdiğimizde gelen hataya bir üst sınır tayin ettik. Bize verilen hatayı (veya daha azını) yapma garantisi veren asgari $a$ değeri de böylece bulunmuş oluyor: $a(\varepsilon) = \sqrt{- \log \varepsilon}$. Karşılaştırma yapabilmek için bir önceki yöntemde kullandığımız hata payını kullanırsak, o zaman $a = \sqrt{8 \log 10} \approx 4,29193205258$ kullanmamız gerektiği ortaya çıkar. Sayısal integrallerde hesaplama maliyeti integrasyon aralığı ile doğru orantılı olarak artar. Bir önceki yöntemdeki integrasyon aralığını kabaca dörtte birine düşürdüğümüze göre, yapacağımız işlem sayısını da dörtte birine düşürdüğümüzü rahatlıkla iddia edebiliriz.

Lemma 3: (Laplace.) $0 \notin (a,b)$, $a \neq 0$ ve $b \neq 0$ için $F^{\prime}(x) = xf(x)$ ise, o zaman \begin{equation*} \int \limits _{a}^{b} f(x) dx = \frac{F(b)}{b} - \frac{F(a)}{a} + \int \limits _{a}^{b} \frac{F(x)}{x^{2}}dx \end{equation*} İspat: $f(x) = \frac{1}{x} xf(x)$ yazıp, kısmi integrasyon tekniğinde $u = \frac{1}{x}$ ve $dv = dF = xf(x)dx$ konarak bu teorem kolayca ispatlanabilir.

Şimdi Lemma 3'te $f(x)= 2\exp(-x^{2})$ ve $F^{\prime}(x) = 2x\exp(-x^{2})$ yazdığımızda $F(x) = -\exp(-x^{2})$ olur. O zaman \begin{equation*} \varepsilon = \int \limits _{a}^{\infty} 2 \exp(-x^{2}) dx = \frac{e^{-a^{2}}}{a} - \int \limits _{a}^{\infty} \frac{\exp(-x^{2})}{x^{2}}dx \leq \frac{e^{-a^{2}}}{a} \end{equation*} elde ederiz. Son basamakta integrand $(a,\infty)$ aralığında pozitif olduğu için, integralin de pozitif olduğunu kullandık. Daha önceki analizlerimizle kıyasla yapabilmek için $\varepsilon = \frac{e^{-a^{2}}}{a}$ denklemini $\varepsilon = 10^{-8}$ için çözeceğiz. Bu denklem $e^{-a^{2}} - a \varepsilon = 0$ şeklinde yeniden düzenlenebilir. Bu aşkın (transendental) denklemi örneğin Raphson-Newton algoritmasıyla yaklaşık olarak çözebiliriz. Raphson-Newton denklemi konumuzun dışında kaldığı için detaylarına girmeyeceğiz ama başlanıç tahmini için $a = \sqrt{-\log \varepsilon}$ kullanıldığında sadece 8 iterasyon sonra sonucun $a = 4,12358554$ değerine yakınsadığını belirtelim.

Öyle görünüyor ki kullandığımız yöntemlerde integrasyon aralığını tayin etme verimliliğinin coğrafi sınırlarına ulaştık. Burası bu postayı bitirmek için iyi bir nokta. Bir sonraki postada integral hesabına kaldığımız yerden devam edeceğiz.