11 Sayısal Türev ve İntegral
Analiz derslerinde türev ve integral, kurallar ve formüllerle hesaplanır. Uygulamada ise fonksiyon çoğu zaman bir formül olarak değil, bir ölçüm tablosu ya da başka bir programın çıktısı olarak elimizdedir. Formül olsa bile ilkel her zaman bulunamaz: \(e^{-x^2}\) ve \(\sin x / x\) fonksiyonlarının ilkelleri elementer fonksiyonlarla yazılamaz. SymPy ile Sembolik Hesap bölümünde SymPy’nin bu integraller için erf ve Si adlı özel fonksiyonlara başvurduğunu görmüştük. Böyle durumlarda türev ve integral, fonksiyonun sonlu sayıda noktadaki değerlerinden sayısal olarak hesaplanır.
Bu bölümde önce fark bölümleriyle türev alacağız ve hatanın adım boyuna bağlı davranışını inceleyeceğiz: adımı küçültmek bir noktaya kadar işe yarar, sonra sonucu bozar. Ardından ızgara verisinin türevini veren np.gradient fonksiyonunu, yamuk ve Simpson kurallarını (önce kendi kodumuzla, sonra np.trapezoid ve scipy.integrate.simpson ile), SciPy’nin hata tahmini de veren uyarlamalı quad fonksiyonunu ve katlı integraller için dblquad ile tplquad fonksiyonlarını kullanacağız. Formüllerin türetilmesi ve hata terimlerinin ispatı Nümerik Analiz notlarındadır (bkz. Nümerik Analiz ve Nümerik Analiz). Burada onları Python’da uygulayıp sınıyoruz.
Vektörizasyon ve Broadcasting bölümünde fark bölümlerini dizilerle hesaplamış ve merkezi farkın hatasının \(h^2\) ile küçüldüğünü görmüştük. Şimdi aynı formüllerin ne kadar doğru olabileceğini ve doğruluğu neyin sınırladığını soruyoruz.
11.1 Fark Bölümleriyle Türev
Türevin tanımındaki limiti küçük ama sıfır olmayan bir adımda durdurmak, akla gelen ilk yaklaşımdır.
Tanım 11.1 (İleri ve Geri Fark Bölümü) \(f\), \(x\)’i içeren bir aralıkta tanımlı ve \(h > 0\) olsun.
\[D_h^{+} f(x) = \frac{f(x + h) - f(x)}{h}, \qquad D_h^{-} f(x) = \frac{f(x) - f(x - h)}{h}\]
sayılarına \(f\)’nin \(x\)’teki ileri fark bölümü ve geri fark bölümü (forward and backward difference quotient), \(h\) sayısına da adım (step size) denir.
Yani ileri fark \(x\)’in sağındaki, geri fark ise solundaki bir noktayla kurulan kirişin eğimidir. \(h \to 0\) iken ikisi de \(f'(x)\)’e yaklaşır; türevin tanımı tam olarak budur.
Tanım 11.2 (Merkezi Fark Bölümü) \(f\), \(x\)’i içeren bir aralıkta tanımlı ve \(h > 0\) olsun.
\[D_h f(x) = \frac{f(x + h) - f(x - h)}{2h}\]
sayısına \(f\)’nin \(x\)’teki merkezi fark bölümü (central difference quotient) denir.
Yani merkezi fark, \(x\)’in iki yanında eşit uzaklıktaki noktalardan geçen kirişin eğimidir. İleri ve geri fark bölümlerinin aritmetik ortalamasına eşittir:
\[D_h f(x) = \frac{D_h^{+} f(x) + D_h^{-} f(x)}{2}.\]
İki formülü Python fonksiyonu olarak yazıp \(f(x) = e^x\) için \(x = 1\)’de deneyelim. Kesin türev \(f'(1) = e\)’dir.
import math
def forward_diff(f, x, h):
"""İleri fark bölümü: (f(x + h) - f(x)) / h."""
return (f(x + h) - f(x)) / h
def central_diff(f, x, h):
"""Merkezi fark bölümü: (f(x + h) - f(x - h)) / (2h)."""
return (f(x + h) - f(x - h)) / (2 * h)
x = 1.0
exact = math.exp(x) # (e^x)' = e^x
print(" h ileri fark hatası merkezi fark hatası")
for h in [0.1, 0.01, 0.001]:
e_fwd = abs(forward_diff(math.exp, x, h) - exact)
e_cen = abs(central_diff(math.exp, x, h) - exact)
print(f"{h:6.0e} {e_fwd:17.3e} {e_cen:19.3e}")Çıktı:
h ileri fark hatası merkezi fark hatası
1e-01 1.406e-01 4.533e-03
1e-02 1.364e-02 4.530e-05
1e-03 1.360e-03 4.530e-07Hata ileri farkta \(h\) ile, merkezi farkta \(h^2\) ile orantılı küçülüyor: \(h\) on kat küçülünce birinci sütun on kat, ikinci sütun yüz kat küçülüyor. Bunun nedeni Taylor açılımıdır.
Önerme 11.1 (Fark Bölümlerinin Kesme Hatası) \(h > 0\) olsun.
- \(f\), \([x, x + h]\) aralığında \(C^2\) sınıfındaysa \(x < \xi < x + h\) olan bir \(\xi\) için \[D_h^{+} f(x) - f'(x) = \frac{h}{2}\, f''(\xi)\] olur.
- \(f\), \([x - h, x + h]\) aralığında \(C^3\) sınıfındaysa \(x - h < \eta < x + h\) olan bir \(\eta\) için \[D_h f(x) - f'(x) = \frac{h^2}{6}\, f'''(\eta)\] olur.
Birinci eşitliğin ispatı için bkz. Nümerik Analiz. İkinci eşitlik oradaki üç-nokta orta nokta formülüdür (bkz. Nümerik Analiz).
Yani hata ifadesindeki \(h\)’nin kuvveti yaklaşımın mertebesini (order) verir: ileri fark bölümünün hatası \(O(h)\), merkezi farkınki \(O(h^2)\)’dir (bkz. Nümerik Analiz). İleri farka birinci mertebeden, merkezi farka ikinci mertebeden bir yaklaşım denir. \(f(x) = e^x\) için \(f'' = f''' = e^x\) olduğundan \(x = 1\) civarında hatalar yaklaşık \(\frac{e}{2}h\) ve \(\frac{e}{6}h^2\)’dir. \(h = 0{,}001\) için bunlar yaklaşık \(1{,}36 \cdot 10^{-3}\) ve \(4{,}53 \cdot 10^{-7}\) eder; tablodaki sayılarla uyuşur.
11.2 Adımı Seçmek: Kesme ve Yuvarlama Hatası
Hata formülleri, \(h\)’yi olabildiğince küçük seçmeyi önerir gibidir. Deneyelim: aynı türevi \(h = 10^{-1}, 10^{-2}, \dots, 10^{-15}\) adımlarıyla hesaplayalım.
import math
x = 1.0
exact = math.exp(x)
print(" h ileri fark hatası merkezi fark hatası")
for k in range(1, 16):
h = 10.0**-k
fwd = (math.exp(x + h) - math.exp(x)) / h
cen = (math.exp(x + h) - math.exp(x - h)) / (2 * h)
print(f"{h:6.0e} {abs(fwd - exact):17.3e}"
f" {abs(cen - exact):19.3e}")Çıktı:
h ileri fark hatası merkezi fark hatası
1e-01 1.406e-01 4.533e-03
1e-02 1.364e-02 4.530e-05
1e-03 1.360e-03 4.530e-07
1e-04 1.359e-04 4.531e-09
1e-05 1.359e-05 5.859e-11
1e-06 1.359e-06 1.635e-10
1e-07 1.399e-07 5.859e-11
1e-08 6.603e-09 6.603e-09
1e-09 2.154e-07 6.603e-09
1e-10 1.548e-06 6.727e-07
1e-11 3.263e-05 1.043e-05
1e-12 4.323e-04 2.103e-04
1e-13 4.559e-04 4.559e-04
1e-14 9.338e-03 9.338e-03
1e-15 3.903e-01 1.683e-01Tablo beklenmedik bir şey söylüyor. İleri farkın hatası \(h = 10^{-8}\)’e kadar küçülüyor, sonra büyümeye başlıyor; \(h = 10^{-15}\)’te türevin ilk ondalık basamağı bile yanlış. Merkezi fark en iyi sonucunu \(h = 10^{-5}\) civarında veriyor; \(h = 10^{-7}\)’de aynı küçüklükte bir hatanın çıkması ise bir rastlantıdır ve nedenini birazdan göreceğiz. Bozulmanın kaynağı paydaki çıkarmadır:
import math
h = 1e-12
a = math.exp(1.0 + h)
b = math.exp(1.0)
print(f"e^(1 + h) = {a:.17f}")
print(f"e^1 = {b:.17f}")
print(f"fark = {a - b:.6e}")
print(f"e * h = {math.e * h:.6e}")
print(f"(1 + h) - 1 = {(1.0 + h) - 1.0:.6e}")Çıktı:
e^(1 + h) = 2.71828182846176380
e^1 = 2.71828182845904509
fark = 2.718714e-12
e * h = 2.718282e-12
(1 + h) - 1 = 1.000089e-12\(e^{1 + h}\) ile \(e^1\)’in ilk on ondalık basamağı aynıdır. İki sayı da makinede yaklaşık 16 anlamlı basamakla saklandığı için farklarının yalnız ilk dört basamağı doğru çıkar: \(2{,}718282 \cdot 10^{-12}\) olması gereken fark \(2{,}718714 \cdot 10^{-12}\) bulundu. Bu, Python ile İlk Adımlar bölümünde gördüğümüz sadeleşme hatasıdır: birbirine çok yakın iki sayının farkında anlamlı basamaklar kaybolur. Adımın kendisi bile tam değildir: \(1 + h\) makinede tam gösterilemediği için \((1 + h) - 1\) farkı \(10^{-12}\) yerine \(1{,}000089 \cdot 10^{-12}\) çıkar. Bu hataların etkisi bölme sırasında \(1/h\) ile büyür.
Önerme 11.2 (Yuvarlama Hatasıyla Birlikte Toplam Hata) \(f\)’nin değerleri makinede \(\tilde f\) olarak hesaplansın ve \(x\)’in yakınındaki her \(t\) için \(|\tilde f(t) - f(t)| \le \delta\) olsun. Çıkarma ve bölmenin kendi yuvarlamasını ihmal edelim. \(\tilde D_h^{+} f(x)\) ve \(\tilde D_h f(x)\), fark bölümlerinin \(\tilde f\) değerleriyle hesaplanan karşılıkları olsun.
- \(x\)’in yakınındaki her \(t\) için \(|f''(t)| \le M_2\) ise \[\big|\tilde D_h^{+} f(x) - f'(x)\big| \le \frac{M_2\, h}{2} + \frac{2\delta}{h}\] olur. Sağ taraf en küçük değerini \(h_{+} = 2\sqrt{\delta / M_2}\) adımında alır.
- \(x\)’in yakınındaki her \(t\) için \(|f'''(t)| \le M_3\) ise \[\big|\tilde D_h f(x) - f'(x)\big| \le \frac{M_3\, h^2}{6} + \frac{\delta}{h}\] olur. Sağ taraf en küçük değerini \(h_0 = \sqrt[3]{3\delta / M_3}\) adımında alır.
İspat
İleri fark. \(\tilde f(t) = f(t) + r(t)\) yazalım; \(|r(t)| \le \delta\)’dır. Bu durumda
\[\tilde D_h^{+} f(x) = D_h^{+} f(x) + \frac{r(x + h) - r(x)}{h}\]
olur. Önerme 11.1 gereği birinci terimin \(f'(x)\)’ten farkı \(\frac{h}{2}|f''(\xi)| \le \frac{M_2 h}{2}\)’dir. İkinci terimin mutlak değeri en çok \(\frac{|r(x + h)| + |r(x)|}{h} \le \frac{2\delta}{h}\)’dir. Üçgen eşitsizliği birinci sınırı verir.
Merkezi fark. Aynı ayrıştırmayla
\[\tilde D_h f(x) - f'(x) = \big(D_h f(x) - f'(x)\big) + \frac{r(x + h) - r(x - h)}{2h}\]
olur. Parantezin mutlak değeri en çok \(\frac{M_3 h^2}{6}\), ikinci terimin mutlak değeri en çok \(\frac{2\delta}{2h} = \frac{\delta}{h}\)’dir.
En küçük değer. \(A, B > 0\) için \(g(h) = Ah + B/h\) fonksiyonunun türevi \(g'(h) = A - B/h^2\)’dir. Bu türev yalnız \(h = \sqrt{B/A}\)’da sıfır olur ve \(g''(h) = 2B/h^3 > 0\) olduğundan burası \((0, \infty)\)’daki en küçük değerdir. \(A = M_2/2\) ve \(B = 2\delta\) için \(h_{+} = \sqrt{4\delta/M_2} = 2\sqrt{\delta/M_2}\) bulunur. Benzer biçimde \(g(h) = Ah^2 + B/h\) için \(g'(h) = 2Ah - B/h^2 = 0\) denklemi \(h^3 = B/(2A)\) verir ve \(g''(h) = 2A + 2B/h^3 > 0\)’dır. \(A = M_3/6\) ve \(B = \delta\) için \(h_0^3 = 3\delta/M_3\) olur.
\(\blacksquare\)
Yani toplam hata iki kaynaktan gelir: \(h\) küçüldükçe azalan kesme hatası (truncation error) ve \(h\) küçüldükçe büyüyen yuvarlama hatası (round-off error) (bkz. Nümerik Analiz ve Nümerik Analiz). Makinede \(\delta\) kabaca \(\varepsilon\,|f(x)|\) kadardır; burada \(\varepsilon = 2^{-52} \approx 2{,}2 \cdot 10^{-16}\) makine epsilonudur. \(f\) ve türevleri \(1\) mertebesindeyse en iyi adımlar ileri fark için \(h_{+} \approx \sqrt{\varepsilon} \approx 10^{-8}\), merkezi fark için \(h_0 \approx \sqrt[3]{\varepsilon} \approx 10^{-5}\) olur. Bu adımlarda hata sınırları da yaklaşık \(\sqrt{\varepsilon} \approx 10^{-8}\) ve \(\varepsilon^{2/3} \approx 10^{-11}\)’dir. Demek ki merkezi fark yalnız daha hızlı yakınsamakla kalmaz, ulaşabildiği doğruluk da daha yüksektir.
Örnek 11.1 (Üstel Fonksiyonda En İyi Adım) \(f(x) = e^x\) için \(x = 1\)’de \(\delta = \varepsilon e\) ve \(M_2 = M_3 = e\) alın. Önerme 11.2 ile verilen en iyi adımları ve bu adımlardaki hata sınırlarını hesaplayıp adım taramasının tablosuyla karşılaştırın.
Çözüm
Formüller. \(h_{+} = 2\sqrt{\delta/M_2}\) ve \(h_0 = \sqrt[3]{3\delta/M_3}\) adımlarını ve sınırların bu adımlardaki değerlerini doğrudan koda dökeriz:
import math
eps = 2.0**-52 # makine epsilonu
x = 1.0
delta = eps * math.exp(x) # bir f değerindeki yuvarlama hatası
M = math.exp(x) # x yakınında |f''| ve |f'''| yaklaşık e
h_fwd = 2 * math.sqrt(delta / M)
bound_fwd = M * h_fwd / 2 + 2 * delta / h_fwd
h_cen = (3 * delta / M) ** (1 / 3)
bound_cen = M * h_cen**2 / 6 + delta / h_cen
print(f"ileri fark : h = {h_fwd:.2e} sınır = {bound_fwd:.2e}")
print(f"merkezi fark: h = {h_cen:.2e} sınır = {bound_cen:.2e}")Çıktı:
ileri fark : h = 2.98e-08 sınır = 8.10e-08
merkezi fark: h = 8.73e-06 sınır = 1.04e-10Karşılaştırma. İleri fark için önerilen adım \(h_{+} \approx 3 \cdot 10^{-8}\)’dir; tabloda en küçük hata (\(6{,}6 \cdot 10^{-9}\)) \(h = 10^{-8}\)’de çıkmıştı. Merkezi fark için önerilen adım \(h_0 \approx 8{,}7 \cdot 10^{-6}\)’dır; tabloda en küçük hata (\(5{,}9 \cdot 10^{-11}\)) \(h = 10^{-5}\)’te çıkmıştı. Aynı hata \(h = 10^{-7}\)’de de çıkmıştı, ama orada yuvarlama hataları tesadüfen birbirini büyük ölçüde götürmüştür. Sınır \(h = 10^{-7}\)’de \(6 \cdot 10^{-9}\)’dur, yani bu adımda hata yüz kat büyük de çıkabilirdi; nitekim aradaki \(h = 10^{-6}\) adımında hata yeniden \(1{,}6 \cdot 10^{-10}\)’dur. Güvenilir olan yalnız \(h = 10^{-5}\)’teki sonuçtur. Gerçek hatalar sınırların altında kalır, çünkü sınır en kötü durumu varsayar: iki yuvarlama hatasının hem en büyük hem de ters işaretli olduğunu. Yine de sınırlar en iyi adımın büyüklüğünü doğru tahmin eder.
\(\blacksquare\)
Aşağıdaki şekil, tablonun iki sütununu ve teoremdeki iki sınırı aynı logaritmik eksenlerde gösteriyor. Hata eğrileri V biçimindedir: sağ kolda kesme hatası, sol kolda yuvarlama hatası baskındır.
Fark bölümlerinde \(h\) küçüldükçe yuvarlama hatası \(1/h\) gibi büyür. \(10^{-12}\) gibi çok küçük bir adım, \(10^{-5}\) gibi orta bir adımdan çok daha kötü sonuç verir. Ayrıca adım, \(x\)’in ve \(f\)’nin ölçeğine göre seçilmelidir: \(x = 10^6\) civarında \(10^{-5}\) gibi bir adım \(x\)’in son basamaklarına düşer; \(x + h\) toplamında adımın yalnız birkaç anlamlı basamağı korunur.
- Merkezi farkı seç. İleri ya da geri farkı yalnız fonksiyon \(x\)’in bir yanında tanımlı değilse kullan.
- Adımı ölçekle: \(\delta \approx \varepsilon\,|f(x)|\) ve \(M_3 \approx |f'''(x)|\) alarak \(h = \sqrt[3]{3\delta/M_3}\) hesapla. Bu bilgiler yoksa \(h = \sqrt[3]{\varepsilon}\,\max(1, |x|)\) kaba ama makul bir başlangıçtır.
- Türevi \(h\) ve \(2h\) adımlarıyla hesaplayıp karşılaştır. İki sonucun farkı hatanın büyüklüğü hakkında fikir verir; fark beklenenden büyükse adımı yeniden düşün.
Örnek 11.2 (Büyük Bir Noktada Logaritmanın Türevi) \(f(x) = \ln x\) fonksiyonunun \(x = 10^6\) noktasındaki türevini merkezi farkla, adımı yukarıdaki üç adıma göre seçerek hesaplayın. Sonucu ezbere seçilmiş \(h = 10^{-5}\) adımının sonucuyla karşılaştırın. Kesin değer \(f'(10^6) = 10^{-6}\)’dır.
Çözüm
1. adım. \(\ln x\), \(10^6\)’nın iki yanında da tanımlıdır; merkezi farkı kullanırız.
2. adım. \(\delta = \varepsilon \ln(10^6) \approx 3{,}1 \cdot 10^{-15}\)’tir. \(f'''(x) = 2/x^3\) olduğundan \(M_3 = 2 \cdot 10^{-18}\) alırız. Önerilen adım \(h = \sqrt[3]{3\delta/M_3} \approx 16{,}6\) çıkar. Bu adım ilk bakışta büyük görünür, ama \(10^6\)’nın yanında küçüktür: \(\ln x\) bu ölçekte çok yavaş değişir.
3. adım. \(h\) ile \(2h\) adımlarının sonuçlarını karşılaştırırız. Kıyas için \(h = 10^{-5}\) adımını da deneriz:
import math
def central_diff(f, x, h):
return (f(x + h) - f(x - h)) / (2 * h)
eps = 2.0**-52
x = 1.0e6
exact = 1 / x # (ln x)' = 1/x
# 2. adım: delta = eps |ln x| ve M3 = |f'''(x)| = 2 / x^3
delta = eps * abs(math.log(x))
M3 = 2 / x**3
h = (3 * delta / M3) ** (1 / 3)
print(f"önerilen h = {h:.2f}")
# 3. adım: h ile 2h karşılaştırılır; kıyas için ezbere h = 1e-5
d1 = central_diff(math.log, x, h)
d2 = central_diff(math.log, x, 2 * h)
print(f"D(h) ile D(2h) arasındaki bağıl fark: {abs(d1 - d2) / d1:.1e}")
for step in [h, 2 * h, 1e-5]:
d = central_diff(math.log, x, step)
rel = abs(d - exact) / exact
print(f"h = {step:8.3g} bağıl hata = {rel:.1e}")Çıktı:
önerilen h = 16.63
D(h) ile D(2h) arasındaki bağıl fark: 2.7e-10
h = 16.6 bağıl hata = 1.0e-10
h = 33.3 bağıl hata = 3.7e-10
h = 1e-05 bağıl hata = 8.3e-08Sonuç. \(D(h)\) ile \(D(2h)\) arasındaki bağıl fark \(2{,}7 \cdot 10^{-10}\)’dur; gerçek bağıl hata da bu mertebededir (\(1{,}0 \cdot 10^{-10}\)). Teoremdeki sınır, \(h = 10^{-5}\) adımında \(\delta/h\) teriminden gelen \(3 \cdot 10^{-4}\) kadar bir bağıl hataya izin verir. Bu çalıştırmada yuvarlama hataları kısmen birbirini götürdüğü için gerçek hata daha küçük çıktı, ama yine de önerilen adımdakinin yaklaşık 800 katıdır.
\(\blacksquare\)
11.3 Izgarada Türev: np.gradient
Şimdiye kadar \(f\)’yi istediğimiz noktada hesaplayabiliyorduk. Fonksiyon bir tablo olarak verildiğinde ise adımı biz seçemeyiz; elimizdeki noktalarla yetinmek zorundayız.
Tanım 11.3 (np.gradient ile Izgarada Türev) \(x_0 < x_1 < \dots < x_n\) noktalarında \(y_k = f(x_k)\) değerleri verilsin. np.gradient(y, x) çağrısı, \(k\) indisli elemanı \(f'(x_k)\)’ye bir yaklaşım olan \(n + 1\) elemanlı bir dizi döndürür.
- İç noktalarda (\(0 < k < n\)) \(x_{k-1}\), \(x_k\), \(x_{k+1}\) noktalarından geçen ikinci dereceden polinomun \(x_k\)’deki türevini alır. Aralıklar eşitse bu, merkezi fark bölümüdür: \((y_{k+1} - y_{k-1})/(2h)\).
- Uç noktalarda varsayılan olarak tek yönlü fark bölümünü kullanır: \((y_1 - y_0)/(x_1 - x_0)\) ve \((y_n - y_{n-1})/(x_n - x_{n-1})\).
edge_order=2verilirse ilk ya da son üç noktadan geçen ikinci dereceden polinomun uçtaki türevini alır.
Noktalar eşit aralıklıysa x dizisi yerine adım da verilebilir: np.gradient(y, h).
Yani np.gradient iç noktalarda ikinci mertebeden, uçlarda ise varsayılan hâlde yalnız birinci mertebeden bir yaklaşım verir; edge_order=2 uçları da ikinci mertebeye çıkarır. Eşit aralıkta bunlar Nümerik Analiz notlarındaki üç-nokta formülleridir (bkz. Nümerik Analiz). Dokuz noktalı bir ızgarada \(\cos x\)’in türevini hesaplayıp kesin türev \(-\sin x\) ile karşılaştıralım:
import numpy as np
x = np.linspace(0, np.pi, 9) # adım pi/8
y = np.cos(x)
exact = -np.sin(x)
g1 = np.gradient(y, x) # uçlarda birinci mertebe
g2 = np.gradient(y, x, edge_order=2) # uçlarda ikinci mertebe
print(np.round(g1, 4))
print(np.round(g2, 4))
print(np.round(exact, 4) + 0.0) # + 0.0: -0. yerine 0. yazar
err1 = np.abs(g1 - exact)
err2 = np.abs(g2 - exact)
print("iç noktalarda en büyük hata:", round(err1[1:-1].max(), 4))
print("uçlarda hata, varsayılan :", np.round(err1[[0, -1]], 4))
print("uçlarda hata, edge_order=2 :", np.round(err2[[0, -1]], 4))Çıktı:
[-0.1938 -0.3729 -0.6891 -0.9003 -0.9745 -0.9003 -0.6891 -0.3729 -0.1938]
[-0.0148 -0.3729 -0.6891 -0.9003 -0.9745 -0.9003 -0.6891 -0.3729 -0.0148]
[ 0. -0.3827 -0.7071 -0.9239 -1. -0.9239 -0.7071 -0.3827 0. ]
iç noktalarda en büyük hata: 0.0255
uçlarda hata, varsayılan : [0.1938 0.1938]
uçlarda hata, edge_order=2 : [0.0148 0.0148]İç noktalardaki en büyük hata \(0{,}0255\)’tir. Bu, Önerme 11.1 ile uyumludur: \(h = \pi/8\) için \(\frac{h^2}{6} \approx 0{,}0257\) eder. Uçlarda ise varsayılan tek yönlü fark \(0{,}19\) kadar sapar; edge_order=2 bu sapmayı \(0{,}015\)’e indirir.
np.gradient kodunun dokuz noktada bulduğu türevler ve kesin türev −sin x. İç noktalardaki merkezi farklar eğriye çok yakındır. Uçlardaki tek yönlü fark ise varsayılan hâlde 0,19 kadar sapar (içi boş halkalar); edge_order=2 ile bu sapma 0,015'e iner (yeşil karolar).Örnek 11.3 (Düzensiz Aralıklı Ölçümlerden Hız) Bir cismin konumu \(t = 0;\ 0{,}5;\ 1{,}2;\ 2;\ 2{,}5;\ 3\) saniyelerinde ölçülmüş olsun ve konum \(s(t) = 2t^2\) metre kuralına uysun. Bu altı ölçümden cismin hızını np.gradient ile hesaplayın ve kesin hız \(v(t) = 4t\) ile karşılaştırın.
Çözüm
Çağrı. Ölçüm anları eşit aralıklı olmadığı için zaman dizisini ikinci argüman olarak veririz. Karşılaştırma için son satırda t dizisini bilerek vermiyoruz:
import numpy as np
t = np.array([0.0, 0.5, 1.2, 2.0, 2.5, 3.0]) # saniye
s = 2 * t**2 # metre
print(np.gradient(s, t))
print(np.round(np.gradient(s, t, edge_order=2), 6))
print(4 * t) # kesin hız
print(np.gradient(s)) # YANLIŞ: adım 1 sanılırÇıktı:
[ 1. 2. 4.8 8. 10. 11. ]
[ 0. 2. 4.8 8. 10. 12. ]
[ 0. 2. 4.8 8. 10. 12. ]
[0.5 1.44 3.75 4.81 5. 5.5 ]Sonuç. İç noktalarda sonuç tam olarak kesin hızdır, çünkü üç noktadan geçen ikinci dereceden polinom \(s\)’nin kendisidir. Uçlarda varsayılan tek yönlü fark \(1\) ve \(11\) verir (kesin değerler \(0\) ve \(12\)); edge_order=2 uçları da düzeltir. Son satır sık yapılan bir hatadır: x verilmezse np.gradient noktalar arasındaki uzaklığı \(1\) sayar ve sonuç anlamsız olur.
\(\blacksquare\)
np.gradient çok boyutlu dizilerde her eksen boyunca ayrı ayrı çalışır ve kısmi türevleri bir demet olarak döndürür. \(f(x, y) = \sin x \cos y\) fonksiyonunun kısmi türevlerini \([0, 2] \times [0, 1]\) üzerindeki bir ızgarada hesaplayalım:
import numpy as np
x = np.linspace(0, 2, 201) # adım 0.01
y = np.linspace(0, 1, 101) # adım 0.01
X, Y = np.meshgrid(x, y, indexing="ij") # X[i, j] = x[i], Y[i, j] = y[j]
F = np.sin(X) * np.cos(Y)
Fx, Fy = np.gradient(F, x, y, edge_order=2) # eksen 0 -> x, eksen 1 -> y
i, j = 100, 50
print(X[i, j], Y[i, j])
print(Fx[i, j], np.cos(1) * np.cos(0.5))
print(Fy[i, j], -np.sin(1) * np.sin(0.5))
print(np.abs(Fx - np.cos(X) * np.cos(Y)).max())
print(np.abs(Fy + np.sin(X) * np.sin(Y)).max())Çıktı:
1.0 0.5
0.4741519791538522 0.4741598817790379
-0.40341595643361927 -0.4034226801113349
3.333216667877892e-05
2.791296891380135e-05indexing="ij" seçeneği, dizinin \(0\) numaralı ekseninin \(x\)’e, \(1\) numaralı ekseninin \(y\)’ye karşılık gelmesini sağlar; böylece np.gradient(F, x, y) önce \(\partial f/\partial x\)’i, sonra \(\partial f/\partial y\)’yi verir. \((1;\ 0{,}5)\) noktasında iki kısmi türevin de ilk dört ondalık basamağı doğrudur. Bütün ızgaradaki en büyük hata \(3 \cdot 10^{-5}\) civarındadır; adım \(0{,}01\) olduğundan bu, ikinci mertebeden bir yöntemden beklenen büyüklüktür.
11.4 Yamuk ve Simpson Kuralları
İntegrale geçelim. Belirli integrali hesaplamak için \([a, b]\) aralığını \(n\) eşit parçaya böler, her parçada \(f\)’yi integrali kolay alınan bir polinomla değiştiririz.
Tanım 11.4 (Bileşik Yamuk Kuralı) \(n \ge 1\), \(h = (b - a)/n\) ve \(k = 0, 1, \dots, n\) için \(x_k = a + kh\) olsun.
\[T_n = h\left[\frac{f(x_0)}{2} + f(x_1) + \dots + f(x_{n-1}) + \frac{f(x_n)}{2}\right]\]
sayısına \(\int_a^b f(x)\,dx\) için bileşik yamuk kuralı (composite trapezoidal rule) yaklaşımı denir.
Yani her \([x_k, x_{k+1}]\) parçasında \(f\)’yi uçlardaki değerleri birleştiren doğruyla değiştirip \(n\) yamuğun alanını toplarız (tek parça için bkz. Nümerik Analiz). İç noktalar iki yamuğa birden katıldığı için ağırlıkları \(h\), uç noktaların ağırlığı \(h/2\)’dir.
Tanım 11.5 (Bileşik Simpson Kuralı) \(n\) çift, \(h = (b - a)/n\) ve \(k = 0, 1, \dots, n\) için \(x_k = a + kh\) olsun.
\[ \begin{aligned} S_n = \frac{h}{3}\Big[&f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + \dots\\[1mm] &\quad + 2f(x_{n-2}) + 4f(x_{n-1}) + f(x_n)\Big] \end{aligned} \]
sayısına \(\int_a^b f(x)\,dx\) için bileşik Simpson kuralı (composite Simpson’s rule) yaklaşımı denir.
Yani \([x_0, x_2], [x_2, x_4], \dots\) ikili parçalarının her birinde \(f\)’yi üç noktadan geçen bir parabolle değiştirip parabollerin altındaki alanları toplarız (tek parça için bkz. Nümerik Analiz). Ağırlıklar tek indisli noktalarda \(4h/3\), iç çift indisli noktalarda \(2h/3\), uçlarda \(h/3\)’tür.
İki kuralı NumPy dilimleriyle birkaç satırda yazabiliriz: y[1:-1] iç noktaları, y[1:-1:2] tek indisli noktaları, y[2:-1:2] ise iç çift indisli noktaları seçer. Örnek olarak ilkeli elementer olmayan
\[\int_0^2 e^{-x^2}\,dx = \frac{\sqrt{\pi}}{2}\,\mathrm{erf}(2)\]
integralini alalım; kesin değeri math.erf ile hesaplarız.
import math
import numpy as np
def trapezoid_rule(f, a, b, n):
"""Bileşik yamuk kuralı, n eşit alt aralık."""
x = np.linspace(a, b, n + 1)
y = f(x)
h = (b - a) / n
return h * (y[0] / 2 + y[1:-1].sum() + y[-1] / 2)
def simpson_rule(f, a, b, n):
"""Bileşik Simpson kuralı; n çift olmalı."""
if n % 2:
raise ValueError("Simpson kuralı için n çift olmalı")
x = np.linspace(a, b, n + 1)
y = f(x)
h = (b - a) / n
return h / 3 * (y[0] + 4 * y[1:-1:2].sum()
+ 2 * y[2:-1:2].sum() + y[-1])
def f(x):
return np.exp(-x**2)
exact = math.sqrt(math.pi) / 2 * math.erf(2)
print(f"kesin değer: {exact:.10f}")
print(f"T_4 = {trapezoid_rule(f, 0, 2, 4):.10f}")
print(f"S_4 = {simpson_rule(f, 0, 2, 4):.10f}")
print(" n yamuk hatası oran Simpson hatası oran")
old_t = old_s = None
for n in [4, 8, 16, 32, 64, 128, 256]:
et = abs(trapezoid_rule(f, 0, 2, n) - exact)
es = abs(simpson_rule(f, 0, 2, n) - exact)
rt = f"{old_t / et:6.2f}" if old_t else " -"
rs = f"{old_s / es:6.2f}" if old_s else " -"
print(f"{n:4d} {et:12.3e} {rt} {es:14.3e} {rs}")
old_t, old_s = et, esÇıktı:
kesin değer: 0.8820813908
T_4 = 0.8806186341
S_4 = 0.8818124253
n yamuk hatası oran Simpson hatası oran
4 1.463e-03 - 2.690e-04 -
8 3.776e-04 3.87 1.588e-05 16.94
16 9.515e-05 3.97 9.942e-07 15.97
32 2.383e-05 3.99 6.212e-08 16.01
64 5.961e-06 4.00 3.882e-09 16.00
128 1.490e-06 4.00 2.426e-10 16.00
256 3.726e-07 4.00 1.516e-11 16.00\(n = 4\) için iki yaklaşımı aşağıdaki şekilde görüyoruz. Simpson kuralı aynı beş fonksiyon değeriyle yamuk kuralından beş kat daha küçük hata veriyor.
Tablonun “oran” sütunları, \(n\) ikiye katlanınca hatanın yamuk kuralında \(4\)’e, Simpson kuralında \(16\)’ya bölündüğünü gösteriyor. Bunun nedeni hata terimleridir.
Önerme 11.3 (Bileşik Kuralların Hatası) \(h = (b - a)/n\) olsun.
- \(f \in C^2[a, b]\) ise bir \(\mu \in (a, b)\) için \[\int_a^b f(x)\,dx - T_n = -\frac{(b - a)\,h^2}{12}\, f''(\mu)\] olur.
- \(n\) çift ve \(f \in C^4[a, b]\) ise bir \(\mu \in (a, b)\) için \[\int_a^b f(x)\,dx - S_n = -\frac{(b - a)\,h^4}{180}\, f^{(4)}(\mu)\] olur.
İspat
Yamuk. Tek bir \([x_k, x_{k+1}]\) parçasında yamuk kuralının hatası, bir \(\xi_k \in (x_k, x_{k+1})\) için \(-\frac{h^3}{12} f''(\xi_k)\)’dir (bkz. Nümerik Analiz). \(n\) parçayı toplarız:
\[\int_a^b f(x)\,dx - T_n = -\frac{h^3}{12}\sum_{k=0}^{n-1} f''(\xi_k).\]
\(\frac{1}{n}\sum_{k} f''(\xi_k)\) ortalaması, \(f''(\xi_k)\) sayılarının en küçüğü ile en büyüğü arasındadır. \(f''\) sürekli olduğundan, ara değer teoremini bu iki değerin alındığı noktalar arasında uygularsak ortalamaya eşit olan bir \(f''(\mu)\) değeri buluruz; \(\mu\) bu iki nokta arasında, dolayısıyla \((a, b)\) içindedir. Böylece toplam \(n f''(\mu)\) olur ve \(n h^3 = (b - a) h^2\) eşitliği birinci formülü verir.
Simpson. Her \([x_{2j}, x_{2j+2}]\) parçasında Simpson kuralının hatası, bir \(\xi_j\) için \(-\frac{h^5}{90} f^{(4)}(\xi_j)\)’dir (bkz. Nümerik Analiz). \(n/2\) parçayı toplayıp aynı ortalama argümanını kullanırız: toplam hata \(-\frac{h^5}{90} \cdot \frac{n}{2}\, f^{(4)}(\mu)\) olur ve \(\frac{n}{2} h^5 = \frac{(b - a) h^4}{2}\) eşitliği ikinci formülü verir.
\(\blacksquare\)
Yani yamuk kuralı ikinci, Simpson kuralı dördüncü mertebedendir. \(n\) ikiye katlanınca \(h\) yarıya iner; yamuk kuralının hatası yaklaşık \(2^2 = 4\)’e, Simpson kuralınınki \(2^4 = 16\)’ya bölünür. Simpson kuralı derecesi en çok \(3\) olan polinomları tam integre eder, çünkü onların dördüncü türevi sıfırdır. Aynı tablo logaritmik eksenlerde iki doğru verir ve doğruların eğimleri mertebeleri gösterir:
11.5 Örneklenmiş Veriyle İntegral
Kendi yazdığımız fonksiyonlar \(f\)’yi bir Python fonksiyonu olarak alıyordu. Bir ölçüm tablosunda ise elimizde yalnız sayılar vardır; NumPy ve SciPy’nin hazır fonksiyonları doğrudan bu sayılarla çalışır.
Tanım 11.6 (np.trapezoid ve scipy.integrate.simpson) y dizisi, x dizisindeki noktalarda fonksiyonun değerlerini tutsun.
np.trapezoid(y, x)her \([x_k, x_{k+1}]\) parçasında bir yamuğun alanını hesaplayıp toplar: \(\sum_k (x_{k+1} - x_k)\,\frac{y_k + y_{k+1}}{2}\). Noktaların eşit aralıklı olması gerekmez. Eşit aralıktaxyerinedx=hde verilebilir.scipy.integrate.simpson(y, x=x)bileşik Simpson kuralını uygular. Alt aralık sayısı tek olduğunda da çalışır; son aralığı ayrı bir formülle hesaba katar.
Yani bu iki fonksiyon, yukarıda kendimiz yazdığımız kuralların veri dizisi alan sürümleridir. NumPy’nin eski sürümlerindeki np.trapz adı kaldırılmıştır; yeni kodda np.trapezoid kullanılır.
import numpy as np
from scipy.integrate import simpson
x = np.linspace(0, 2, 5) # n = 4 alt aralık, h = 0.5
y = np.exp(-x**2)
print(np.trapezoid(y, x)) # x noktaları verilerek
print(np.trapezoid(y, dx=0.5)) # eşit aralıkta yalnız adım
print(simpson(y, x=x))
# Düzensiz aralıklar: yamuk kuralı her parçanın kendi genişliğini kullanır
xu = np.array([0.0, 0.2, 0.5, 1.0, 1.4, 2.0])
print(np.trapezoid(np.exp(-xu**2), xu))Çıktı:
0.8806186341245394
0.8806186341245394
0.8818124252941161
0.8931873236709345İlk üç sayı, kendi kodumuzun \(n = 4\) için bulduğu \(T_4\) ve \(S_4\) değerleriyle aynıdır. Son satırda noktalar düzensiz aralıklıdır; np.trapezoid her parçanın kendi genişliğini kullanır.
- Noktaları ve değerleri kur:
x = np.linspace(a, b, n + 1)vey = f(x). Simpson kuralı için \(n\)’yi çift seç. np.trapezoid(y, x)ya dasimpson(y, x=x)ile yaklaşımı hesapla.- \(n\)’yi ikiye katlayıp yeniden hesapla. Yeni yaklaşımın hatası yamuk kuralında yaklaşık \((T_{2n} - T_n)/3\), Simpson kuralında yaklaşık \((S_{2n} - S_n)/15\)’tir.
Örnek 11.4 (İkiye Katlayarak Hata Tahmini) \(\int_1^2 \frac{dx}{x} = \ln 2\) integraline yamuk ve Simpson kurallarıyla \(n = 8, 16, 32\) alt aralık için yaklaşın. Her ikiye katlamada yeni yaklaşımın hatasını tahmin edin ve tahmini gerçek hatayla karşılaştırın.
Çözüm
Tahminin gerekçesi. \(I\) kesin değer olsun. Yamuk kuralında \(I - T_n \approx 4\,(I - T_{2n})\) olduğundan
\[T_{2n} - T_n = (I - T_n) - (I - T_{2n}) \approx 3\,(I - T_{2n})\]
olur, yani \(I - T_{2n} \approx (T_{2n} - T_n)/3\)’tür. Simpson kuralında aynı hesap \(16\) çarpanıyla \(I - S_{2n} \approx (S_{2n} - S_n)/15\) verir.
Kod. Üç adımı bir döngüde uygularız:
import math
import numpy as np
from scipy.integrate import simpson
exact = math.log(2)
old_t = old_s = None
for n in [8, 16, 32]:
# 1. adım: noktalar ve değerler
x = np.linspace(1, 2, n + 1)
y = 1 / x
# 2. adım: kurallar
t = np.trapezoid(y, x)
s = simpson(y, x=x)
# 3. adım: n'yi ikiye katlayınca hata tahmini
if old_t is not None:
est_t = (t - old_t) / 3
est_s = (s - old_s) / 15
print(f"n = {n:2d} yamuk: tahmin {est_t:.2e}"
f" gerçek {exact - t:.2e}")
print(f" Simpson: tahmin {est_s:.2e}"
f" gerçek {exact - s:.2e}")
old_t, old_s = t, sÇıktı:
n = 16 yamuk: tahmin -2.44e-04 gerçek -2.44e-04
Simpson: tahmin -4.59e-07 gerçek -4.72e-07
n = 32 yamuk: tahmin -6.10e-05 gerçek -6.10e-05
Simpson: tahmin -2.95e-08 gerçek -2.97e-08Sonuç. Yamuk kuralı için tahmin ile gerçek hata, yazılan üç basamakta aynıdır. Simpson kuralında da tahmin gerçek hatanın yüzde üçü içindedir. Hatalar negatiftir, yani iki kural da integrali fazla tahmin eder. Önerme 11.3 bunu açıklar: \(f(x) = 1/x\) için \(f''(x) = 2/x^3\) ve \(f^{(4)}(x) = 24/x^5\) pozitiftir. Yamuk kuralında bunun geometrik anlamı da vardır: \(1/x\) dışbükey olduğu için kirişler eğrinin üstünde kalır. Kesin değeri bilmediğimiz durumlarda bu tahmin, kaç alt aralığın yeteceğine karar vermek için kullanılır.
\(\blacksquare\)
Örnek 11.5 (Hız Tablosundan Yol) Bir aracın hızı ilk bir dakika boyunca beş saniyede bir okunmuş ve sırasıyla 0; 8,5; 14,6; 19,0; 22,1; 24,3; 25,9; 27,1; 27,9; 28,5; 28,9; 29,2; 29,5 m/s bulunmuş olsun. Aracın bu sürede aldığı yolu yamuk ve Simpson kurallarıyla hesaplayın. Tablo \(v(t) = 30\,(1 - e^{-t/15})\) modelinden üretilip bir ondalığa yuvarlanmıştır; sonuçları modelin verdiği \(\int_0^{60} v(t)\,dt = 1350 + 450\,e^{-4}\) değeriyle karşılaştırın.
Çözüm
Veri. 13 ölçüm 12 eşit alt aralık verir; \(n\) çift olduğundan Simpson kuralı doğrudan uygulanır.
import math
import numpy as np
from scipy.integrate import simpson
t = np.arange(0, 61, 5.0) # 0, 5, ..., 60 saniye
v = np.array([0.0, 8.5, 14.6, 19.0, 22.1, 24.3, 25.9,
27.1, 27.9, 28.5, 28.9, 29.2, 29.5]) # m/s
print("yamuk :", np.trapezoid(v, t))
print("Simpson:", simpson(v, x=t))
print("model :", 1350 + 450 * math.exp(-4))Çıktı:
yamuk : 1353.75
Simpson: 1357.8333333333335
model : 1358.2420374999303Sonuç. Simpson kuralı modelin değerinden \(0{,}4\) metre, yamuk kuralı \(4{,}5\) metre uzakta kaldı. Verideki yuvarlama tek başına sonucu \(60 \cdot 0{,}05 = 3\) metreye kadar değiştirebileceğinden Simpson sonucu, verinin izin verdiği doğruluktadır. Yamuk kuralının hatası ise çoğunlukla kuralın kendisinden gelir: \(v\) içbükeydir ve kirişler eğrinin altında kalır.
\(\blacksquare\)
Bazen integralin yalnız son değeri değil, üst sınırın fonksiyonu olarak bütün gidişatı gerekir: \(F(x) = \int_a^x f(t)\,dt\). scipy.integrate.cumulative_trapezoid(y, x, initial=0) her \(x_k\) için \(\int_{x_0}^{x_k} f(t)\,dt\) integralinin yamuk kuralıyla bulunan yaklaşımını verir; initial=0 ilk elemanı \(0\) yapar. Örnek olarak SymPy ile Sembolik Hesap bölümünde karşılaştığımız sinüs integrali fonksiyonunu,
\[\mathrm{Si}(x) = \int_0^x \frac{\sin t}{t}\,dt,\]
\([0, 20]\) aralığında hesaplayalım. \(t = 0\)’da \(\sin t / t\) ifadesi \(0/0\) olur. np.sinc(u) ise \(\sin(\pi u)/(\pi u)\) değerini hesaplar ve \(u = 0\)’da \(1\) verir; bu yüzden np.sinc(t / np.pi) sorunu ortadan kaldırır.
import numpy as np
from scipy.integrate import cumulative_trapezoid
from scipy.special import sici
t = np.linspace(0, 20, 401) # adım 0.05
y = np.sinc(t / np.pi) # sin(t)/t; t = 0'da 1
Si = cumulative_trapezoid(y, t, initial=0) # Si[k] = t[0]..t[k] integrali
exact = sici(t)[0] # sici, (Si, Ci) ikilisini verir
print(len(Si), Si[-1])
print("en büyük hata:", np.abs(Si - exact).max())
k = Si.argmax()
print("en büyük değer:", t[k], Si[k])Çıktı:
401 1.5482454765212605
en büyük hata: 9.086494840016002e-05
en büyük değer: 3.1500000000000004 1.8518598626730711Sonuçlar, SciPy’nin özel fonksiyonu scipy.special.sici ile \(10^{-4}\) duyarlıkla uyuşuyor. sici(t) iki dizi döndürür, sinüs integrali \(\mathrm{Si}\) ve kosinüs integrali \(\mathrm{Ci}\); ilkini [0] indisiyle aldık. En büyük değer \(x = \pi\) civarındadır, çünkü \(\mathrm{Si}\)’nin türevi olan \(\sin x / x\) ilk kez orada sıfır olur ve işaret değiştirir.
cumulative_trapezoid kodunun 401 noktada hesapladığı Si(x) değerleri (mavi) ve integrali alınan sin t/t fonksiyonu (gri). Si(x), sin t/t pozitifken artar, negatifken azalır; en büyük değerine x = π civarında ulaşır ve π/2 doğrusu etrafında sönen salınımlarla ona yaklaşır.Şekil, \(x \to \infty\) iken \(\mathrm{Si}(x)\)’in \(\pi/2\)’ye yaklaştığını düşündürüyor. Bu genelleştirilmiş integrali birazdan quad ile hesaplayacağız.
11.6 Uyarlamalı İntegral: quad
İntegrand bir Python fonksiyonu olarak verildiğinde en pratik araç scipy.integrate.quad fonksiyonudur. Bu fonksiyon noktaları kendisi seçer ve sonucun yanında bir hata tahmini de verir.
Tanım 11.7 (scipy.integrate.quad) quad(f, a, b) çağrısı \(\int_a^b f(x)\,dx\) için bir (değer, hata) ikilisi döndürür; buradaki hata, değerin mutlak hatası için bir tahmindir. quad her alt aralıkta integrali ve hatasını yüksek dereceli bir kuralla tahmin eder ve toplam hata tahmini istenen sınırın altına inene kadar hata tahmini en büyük olan alt aralığı ikiye böler. Sınırlar np.inf ya da -np.inf olabilir. Varsayılan hedef, mutlak ya da bağıl hatanın yaklaşık \(1{,}5 \cdot 10^{-8}\)’in altına inmesidir (epsabs ve epsrel parametreleri).
Yani quad bizden \(n\) istemez: zor bölgelerde alt aralıkları sıklaştırır, kolay bölgelerde geniş tutar. Bu tür yöntemlere uyarlamalı (adaptive) integral yöntemleri denir. quad, uzun yıllardır kullanılan QUADPACK kütüphanesinin algoritmalarına dayanır.
import math
import numpy as np
from scipy.integrate import quad
def f(x):
return np.exp(-x**2)
value, err = quad(f, 0, 2)
exact = math.sqrt(math.pi) / 2 * math.erf(2)
print(value, err)
print("gerçek hata:", abs(value - exact))
info = quad(f, 0, 2, full_output=1)[2]
print("fonksiyon çağrısı:", info["neval"])Çıktı:
0.8820813907624215 9.793070696178202e-15
gerçek hata: 0.0
fonksiyon çağrısı: 21full_output=1 verilirse quad üçüncü bir eleman olarak ayrıntılı bir bilgi sözlüğü de döndürür; kodda bu demetin [2] indisli elemanını aldık. Sözlüğün neval anahtarı, fonksiyonun kaç kez hesaplandığını verir. quad sonucu yalnız 21 fonksiyon değeriyle makine duyarlığında buldu; Simpson kuralı 257 noktayla ancak \(10^{-11}\)’e inmişti. Hata tahmini \(10^{-14}\) mertebesindedir ve gerçek hata ondan da küçüktür.
İntegrand başka parametrelere de bağlıysa bu parametrelerin değerleri args ile verilir.
Örnek 11.6 (Bessel Fonksiyonlarını İntegralle Hesaplamak) \(n\) bir tam sayı olmak üzere birinci tür Bessel fonksiyonu
\[J_n(x) = \frac{1}{\pi}\int_0^\pi \cos(n\tau - x\sin\tau)\,d\tau\]
integraliyle verilebilir. \(J_0(2{,}5)\), \(J_1(2{,}5)\) ve \(J_2(2{,}5)\) değerlerini quad ile hesaplayıp scipy.special.jv fonksiyonunun değerleriyle karşılaştırın.
Çözüm
Parametreler. İntegrand \(\tau\)’nun yanında \(n\) ve \(x\)’e de bağlıdır. quad integrali her zaman fonksiyonun ilk argümanına göre alır; geri kalan argümanların değerleri args demetiyle verilir:
import numpy as np
from scipy.integrate import quad
from scipy.special import jv
def integrand(tau, n, x):
return np.cos(n * tau - x * np.sin(tau))
x = 2.5
for n in [0, 1, 2]:
value, err = quad(integrand, 0, np.pi, args=(n, x))
print(n, value / np.pi, jv(n, x))Çıktı:
0 -0.048383776468198046 -0.048383776468197914
1 0.49709410246427405 0.4970941024642741
2 0.4460590584396172 0.44605905843961724Sonuç. Üç değer de SciPy’nin özel fonksiyonuyla son bir iki basamak dışında aynıdır. args yerine lambda tau: integrand(tau, n, x) de yazılabilirdi; args daha okunaklıdır.
\(\blacksquare\)
Gauss Kuralları
quad’ın 21 değerle bu kadar iyi sonuç vermesinin nedeni, düğümleri akıllıca seçmesidir. Yamuk ve Simpson kuralları düğümleri eşit aralıklı alır. Düğümleri de serbest bırakırsak, \(n\) düğümlü bir kural derecesi en çok \(2n - 1\) olan bütün polinomları tam integre edecek biçimde kurulabilir. Bu kurallara Gauss–Legendre kuralları denir; ispatlarını burada vermeyeceğiz. NumPy’de np.polynomial.legendre.leggauss(n) fonksiyonu \([-1, 1]\) aralığı için düğümleri ve ağırlıkları verir:
import numpy as np
t, w = np.polynomial.legendre.leggauss(3) # 3 düğüm ve ağırlık
print(t)
print(w)
for k in range(7):
approx = (w * t**k).sum()
exact = (1 - (-1) ** (k + 1)) / (k + 1) # [-1, 1] üzerinde x^k
print(f"k = {k} Gauss = {approx:.12f} kesin = {exact:.12f}")Çıktı:
[-0.77459667 0. 0.77459667]
[0.55555556 0.88888889 0.55555556]
k = 0 Gauss = 2.000000000000 kesin = 2.000000000000
k = 1 Gauss = 0.000000000000 kesin = 0.000000000000
k = 2 Gauss = 0.666666666667 kesin = 0.666666666667
k = 3 Gauss = 0.000000000000 kesin = 0.000000000000
k = 4 Gauss = 0.400000000000 kesin = 0.400000000000
k = 5 Gauss = 0.000000000000 kesin = 0.000000000000
k = 6 Gauss = 0.240000000000 kesin = 0.285714285714Üç düğümlü kural \(k = 0, 1, \dots, 5\) için kesin sonucu verir ve \(k = 6\)’da yanılır: \(2 \cdot 3 - 1 = 5\). Düğümler simetriktir ve aralığın uç noktaları düğüm değildir. \(\int_0^2 e^{-x^2}\,dx\) integralini \(x = 1 + t\) değişkeniyle \([-1, 1]\)’e taşıyıp düğüm sayısını artıralım:
import math
import numpy as np
exact = math.sqrt(math.pi) / 2 * math.erf(2)
for n in range(2, 9):
t, w = np.polynomial.legendre.leggauss(n)
# [0, 2] -> [-1, 1]: x = 1 + t, dx = dt
approx = (w * np.exp(-(1 + t) ** 2)).sum()
print(f"n = {n} hata = {abs(approx - exact):.1e}")Çıktı:
n = 2 hata = 3.7e-02
n = 3 hata = 3.2e-03
n = 4 hata = 1.5e-04
n = 5 hata = 3.4e-06
n = 6 hata = 3.5e-08
n = 7 hata = 6.8e-09
n = 8 hata = 3.4e-10Sekiz fonksiyon değeriyle hata \(3 \cdot 10^{-10}\)’a iniyor; Simpson kuralı bu doğruluk için yaklaşık 130 nokta istiyordu. quad her alt aralıkta 21 noktalı bir Gauss–Kronrod kuralı kullanır: bu kural 10 noktalı Gauss kuralının düğümlerine 11 düğüm ekler. Değer olarak 21 noktalı kuralın sonucu alınır; iki kuralın sonuçları arasındaki fark ise hata tahmininin temelidir. quad bu farkı deneyimle belirlenmiş bir formülle ölçekler ve tahminin yuvarlama düzeyinin altına inmesine izin vermez. İlk quad örneğindeki \(9{,}8 \cdot 10^{-15}\) tahmini bu yuvarlama tabanıdır: \(50\,\varepsilon \cdot 0{,}8821 \approx 9{,}8 \cdot 10^{-15}\), yani makine epsilonunun 50 katı ile integralin değerinin çarpımı.
Zor İntegrandlar
Uyarlamalı bölmenin gücü, integrand bir yerde hızla değiştiğinde ortaya çıkar. \(x = 0{,}3\)’te dar ve yüksek bir tepesi olan
\[f(x) = \frac{1}{(x - 0{,}3)^2 + 10^{-4}}\]
fonksiyonunu \([0, 1]\)’de integre edelim. Kesin değer \(100\,(\arctan 70 + \arctan 30)\)’dur. quad’ın hangi noktalara baktığını görmek için fonksiyonu, her çağrıda \(x\)’i bir listeye ekleyecek biçimde yazıyoruz:
import math
import numpy as np
from scipy.integrate import quad
calls = []
def f(x):
calls.append(x) # quad'ın baktığı noktaları kaydet
return 1 / ((x - 0.3)**2 + 1e-4)
exact = 100 * (math.atan(70) + math.atan(30))
value, err, info = quad(f, 0, 1, full_output=1)
print(value, err, abs(value - exact))
print("çağrı sayısı:", len(calls), " alt aralık:", info["last"])
left = info["alist"][:info["last"]] # alt aralıkların sol uçları
ends = np.sort(np.concatenate([left, [1.0]])) # sağ uç 1'i de ekle
print(ends)
x = np.linspace(0, 1, 1001)
print("yamuk, 1001 nokta, hata:", abs(np.trapezoid(f(x), x) - exact))Çıktı:
309.3986915124147 2.3722903395168027e-08 2.2737367544323206e-13
çağrı sayısı: 315 alt aralık: 8
[0. 0.125 0.25 0.28125 0.296875 0.3125 0.375 0.5
1. ]
yamuk, 1001 nokta, hata: 6.6448414486330876e-06Bilgi sözlüğünde info["last"] alt aralık sayısı, info["alist"] alt aralıkların sol uçlarıdır; bu sol uçlara np.concatenate ile son sağ uç \(1\)’i ekleyip sıralayınca bütün uçları elde ederiz. quad 315 fonksiyon değeriyle \(10^{-13}\) hataya ulaştı; 1001 eşit aralıklı noktayla yamuk kuralı ise \(7 \cdot 10^{-6}\) hatada kaldı. Şekil, alt aralıkların tepeye doğru nasıl daraldığını gösteriyor.
quad'ın x = 0,3'teki dar tepeyi nasıl çözdüğü: tepe fonksiyonu kodunun kaydettiği veriler. Kesikli çizgiler quad'ın böldüğü 8 alt aralığın uçları, eksenin altındaki kısa çizgiler fonksiyonun hesaplandığı 315 nokta. Alt aralıklar tepenin yanında ikiye bölüne bölüne daralır; fonksiyonun düz olduğu yerde tek geniş aralık yeter.quad, sonsuz aralıkları ve uç noktada sonsuza giden integrandları da çoğu zaman sorunsuz işler (genelleştirilmiş integraller için bkz. Analiz 3 ve Analiz 3):
import math
import numpy as np
from scipy.integrate import quad
print(quad(lambda x: np.exp(-x**2), 0, np.inf), math.sqrt(math.pi) / 2)
print(quad(lambda x: np.exp(-x**2), -np.inf, np.inf), math.sqrt(math.pi))
print(quad(lambda x: 1 / np.sqrt(x), 0, 1))
print(quad(np.log, 0, 1))Çıktı:
(0.8862269254527579, 7.101318390472462e-09) 0.8862269254527579
(1.7724538509055159, 1.4202636780944923e-08) 1.7724538509055159
(1.9999999999999984, 5.773159728050814e-15)
(-0.9999999999999999, 1.1102230246251563e-15)\(\int_0^\infty e^{-x^2}\,dx = \sqrt{\pi}/2\), \(\int_{-\infty}^\infty e^{-x^2}\,dx = \sqrt{\pi}\), \(\int_0^1 x^{-1/2}\,dx = 2\) ve \(\int_0^1 \ln x\,dx = -1\) değerleri son bir iki basamağa kadar doğru. Sonsuz aralıkta quad önce bir değişken değiştirmeyle aralığı sonlu yapar. Uçtaki tekillikleri ise integrandı uç noktada hiç hesaplamadan, uca doğru sıklaşan alt aralıklar ve bir ekstrapolasyonla aşar. Tekillik aralığın içindeyse quad’a yerini söylemek gerekir:
import numpy as np
from scipy.integrate import quad
def f(x):
return 1 / np.sqrt(abs(x))
with np.errstate(divide="ignore"): # 1/0 uyarısını sustur
print(quad(f, -1, 1))
print(quad(f, -1, 1, points=[0]))Çıktı:
(inf, inf)
(3.9999999999999813, 5.684341886080802e-14)\(1/\sqrt{|x|}\) fonksiyonu \(x = 0\)’da sonsuza gider. 21 noktalı kuralın orta düğümü \([-1, 1]\)’in tam ortası olduğundan quad fonksiyonu \(0\)’da hesaplamaya çalıştı ve inf döndürdü. points=[0] parametresi aralığı tekil noktadan böler. Tekillik böylece iki alt aralığın ucuna gelir ve sonuç doğru çıkar: \(\int_{-1}^1 |x|^{-1/2}\,dx = 4\).
Iraksak bir integralde ise quad yine bir sayı döndürür ve yalnız bir uyarı verir. \(\int_0^1 dx/x\) integrali ıraksaktır:
import warnings
from scipy.integrate import quad
with warnings.catch_warnings(record=True) as caught:
warnings.simplefilter("always") # uyarıları yakala
value, err = quad(lambda x: 1 / x, 0, 1)
print(value, err)
print(caught[0].category.__name__)
print(str(caught[0].message).splitlines()[0])Çıktı:
41.67684067538809 9.35056037314051
IntegrationWarning
The maximum number of subdivisions (50) has been achieved.quad 50 alt aralık sınırına takılıp \(41{,}68\) değerini ve \(9{,}35\) gibi büyük bir hata tahmini döndürdü; sonuç anlamsızdır. warnings.catch_warnings(record=True) bloğu, blok içinde verilen uyarıları bir listede toplar; böylece uyarıyı program içinde denetleyebiliriz.
quad başarısız olduğunda hata vermez; IntegrationWarning türünde bir uyarı basar ve bir sayı döndürmeye devam eder. Hata tahmini değere göre büyükse ya da bir uyarı çıktıysa sonuç güvenilir değildir. Tersine, küçük bir hata tahmini de kanıt değildir: integrand quad’ın baktığı noktaların arasında beklenmedik biçimde davranıyorsa tahmin yanılabilir. Önemli bir sonucu bağımsız bir yolla, örneğin değişken değiştirerek ya da integrali parçalara bölerek yeniden hesaplayın.
Yavaş sönen salınımlar da quad’ı zorlar. Şekildeki \(\mathrm{Si}(x)\)’in limiti olan \(\int_0^\infty \frac{\sin x}{x}\,dx = \frac{\pi}{2}\) integralini (bkz. Analiz 3) doğrudan hesaplamaya çalışmak başarısız olur. Çare, salınımı quad’a ağırlık olarak bildirmektir: weight="sin", wvar=1 seçenekleri \(\int g(x)\sin x\,dx\) biçimindeki integraller için özel bir yöntem kullanır.
import math
import warnings
import numpy as np
from scipy.integrate import quad
def sinc(x):
return np.sinc(x / np.pi) # sin(x)/x
with warnings.catch_warnings(record=True) as caught:
warnings.simplefilter("always")
print(quad(sinc, 0, np.inf))
print(str(caught[0].message).splitlines()[0])
# [0, 1] olağan; [1, oo) üzerinde 1/x, sin(x) ağırlığıyla
part1, _ = quad(sinc, 0, 1)
part2, _ = quad(lambda x: 1 / x, 1, np.inf, weight="sin", wvar=1)
print(part1 + part2, math.pi / 2)Çıktı:
(2.2478679634674115, 3.2903230521899287)
The integral is probably divergent, or slowly convergent.
1.5707963268467806 1.5707963267948966Doğrudan çağrı \(2{,}25\) gibi yanlış bir değer ve “integral muhtemelen ıraksak ya da yavaş yakınsak” uyarısı verdi. İntegrali \([0, 1]\) ve \([1, \infty)\) olarak ikiye bölüp ikinci parçayı \(g(x) = 1/x\) ve \(\sin x\) ağırlığıyla hesaplayınca sonuç \(\pi/2\)’ye \(5 \cdot 10^{-11}\) yakınlıkta çıktı.
- İntegrandı incele: aralık sonsuz mu, integrand bir yerde sonsuza gidiyor mu, yavaş sönen bir salınımı var mı?
- Zorlukları
quad’a bildir: tekil noktaları aralığın ucuna getir (aralığı böl ya dapointsver), sonsuz aralığı ayrı bir parça yap, salınımlar içinweightkullan. - Uyarıları ve hata tahminini oku; sonucu bilinen bir değerle ya da başka bir yolla denetle.
Örnek 11.7 (Gamma Fonksiyonunun Yarımdaki Değeri) \(\Gamma(1/2) = \int_0^\infty t^{-1/2} e^{-t}\,dt\) değerini quad ile hesaplayın ve \(\sqrt{\pi}\) ile karşılaştırın (bkz. Analiz 3).
Çözüm
1. adım. İntegrand \(t = 0\)’da sonsuza gider ve aralık sonsuzdur; iki zorluk birden vardır.
2. adım. İntegrali \(t = 1\)’den bölüp iki parçayı ayrı hesaplarız: \([0, 1]\)’de tekillik aralığın ucundadır, \([1, \infty)\)’da integrand üstel olarak söner.
import math
import numpy as np
from scipy.integrate import quad
def g(t):
return np.exp(-t) / np.sqrt(t)
left, err1 = quad(g, 0, 1) # t = 0'da tekil
right, err2 = quad(g, 1, np.inf) # sonsuz aralık
print(left, err1)
print(right, err2)
print(left + right)
print(math.sqrt(math.pi), math.gamma(0.5))Çıktı:
1.4936482656248504 4.0730441241976223e-10
0.2788055852806195 4.764513549588497e-09
1.77245385090547
1.7724538509055159 1.77245385090551593. adım. Uyarı çıkmadı ve iki parçanın hata tahminleri \(10^{-9}\) mertebesinde kaldı. Toplam, \(\sqrt{\pi}\) ile ve math.gamma(0.5) ile \(5 \cdot 10^{-14}\) yakınlıkta uyuşuyor.
\(\blacksquare\)
11.7 Katlı İntegraller: dblquad ve tplquad
quad’ı iç içe kullanarak katlı integraller de hesaplanabilir; SciPy bunu hazır fonksiyonlarla yapar.
Tanım 11.8 (dblquad ile İki Katlı İntegral) dblquad(func, a, b, gfun, hfun) çağrısı
\[\int_a^b \int_{g(x)}^{h(x)} \mathrm{func}(y, x)\,dy\,dx\]
ardışık integralini hesaplar ve bir (değer, hata) ikilisi döndürür. gfun ve hfun, \(x\)’e bağlı alt ve üst sınırları veren fonksiyonlardır; sınırlar sabitse sayı olarak da verilebilirler. func fonksiyonunun ilk argümanı iç değişken \(y\), ikincisi dış değişken \(x\)’tir.
Yani dblquad, birinci tip bir bölgede (bkz. İntegral Calculus) yazılan ardışık integrali (bkz. İntegral Calculus) iç içe iki quad çağrısıyla hesaplar: dış integral \(x\) üzerinden, iç integral \(y\) üzerinden alınır.
func(y, x) sırası, matematikteki \(f(x, y)\) alışkanlığının tersidir. lambda x, y: ... yazmak hata vermez ama sessizce başka bir integral hesaplar; integrand \(x\) ile \(y\)’de simetrik değilse sonuç yanlış çıkar. tplquad’da sıra func(z, y, x)’tir: en içteki değişken her zaman önce gelir.
Örnek 11.8 (İki Parabol Arasındaki Bölgede İntegral) \(D\), \(y = 2x^2\) ve \(y = 1 + x^2\) parabollerinin sınırladığı bölge olmak üzere \(\iint_D (x + 2y)\,dA\) integralini dblquad ile hesaplayın. Bu integral İntegral Calculus notlarında elle hesaplanmış ve \(\frac{32}{15}\) bulunmuştur (bkz. İntegral Calculus).
dblquad örneğindeki D bölgesi. Dış integral x'i a = −1'den b = 1'e götürür; her x için iç integral y'yi alttaki gfun(x) = 2x² parabolünden üstteki hfun(x) = 1 + x² parabolüne kadar tarar (ok).Çözüm
Sınırlar. Paraboller \(2x^2 = 1 + x^2\), yani \(x = \pm 1\) iken kesişir. \(-1 \le x \le 1\) için \(x^2 \le 1\) olduğundan alttaki parabol \(y = 2x^2\), üstteki \(y = 1 + x^2\)’dir. Dolayısıyla \(a = -1\), \(b = 1\), gfun fonksiyonu \(2x^2\) ve hfun fonksiyonu \(1 + x^2\)’dir.
Kod.
from fractions import Fraction
from scipy.integrate import dblquad
value, err = dblquad(lambda y, x: x + 2 * y, # önce y, sonra x
-1, 1, # x sınırları
lambda x: 2 * x**2, # alttaki parabol
lambda x: 1 + x**2) # üstteki parabol
print(value, err)
print(float(Fraction(32, 15)))Çıktı:
2.1333333333333333 2.3684757858670007e-14
2.1333333333333333Sonuç. dblquad, \(\frac{32}{15} = 2{,}1333\ldots\) değerini \(10^{-14}\) mertebesinde bir hata tahminiyle verdi; elle yapılan hesapla aynı.
\(\blacksquare\)
Tanım 11.9 (tplquad ile Üç Katlı İntegral) tplquad(func, a, b, gfun, hfun, qfun, rfun) çağrısı
\[\int_a^b \int_{g(x)}^{h(x)} \int_{q(x, y)}^{r(x, y)} \mathrm{func}(z, y, x)\,dz\,dy\,dx\]
ardışık integralini hesaplar. gfun ve hfun \(x\)’e, qfun ve rfun ise \(x\) ile \(y\)’ye bağlı sınırlardır; sabit sınırlar sayı olarak verilebilir.
Yani tplquad, tip 1 bir cisim için (bkz. İntegral Calculus) yazılan ardışık integrali hesaplar: en içte \(z\), ortada \(y\), dışta \(x\) üzerinden integral alınır.
Örnek 11.9 (Bir Yüzeyin Altındaki Cisimde Üç Katlı İntegral) \(E\), birinci oktantta \(z = 12xy\) yüzeyi ile \(y = x\), \(x = 1\) ve \(z = 0\) düzlemleri arasında kalan cisim olsun. \(\iiint_E z\,dV\) integralini tplquad ile hesaplayın ve İntegral Calculus notlarında bulunan \(4\) sonucuyla karşılaştırın (bkz. İntegral Calculus).
Çözüm
Sınırlar. Cisim \(0 \le x \le 1\), \(0 \le y \le x\) ve \(0 \le z \le 12xy\) eşitsizlikleriyle betimlenir. Sabit alt sınırları doğrudan 0 olarak yazarız:
from scipy.integrate import tplquad
value, err = tplquad(lambda z, y, x: z,
0, 1, # x sınırları
0, lambda x: x, # y sınırları
0, lambda x, y: 12 * x * y) # z sınırları
print(value, err)Çıktı:
4.0 7.890024717707689e-13Sonuç. Değer \(4\), hata tahmini \(10^{-12}\) mertebesinde; elle bulunan sonuçla aynı.
\(\blacksquare\)
Dört ve daha çok katlı integraller için scipy.integrate.nquad vardır. Ancak boyut arttıkça iç içe kuralların maliyeti üstel büyür: her boyutta 21 nokta kullanılırsa \(d\) boyutta en az \(21^d\) fonksiyon değeri gerekir. Yüksek boyutlarda bu yüzden Monte Carlo yöntemleri tercih edilir (bkz. Rastgele Sayılar ve Monte Carlo Yöntemleri).
11.8 Alıştırmalar
Alıştırma 11.1 (Küp Fonksiyonunda İleri Farkın Hatası) \(f(x) = x^3\) için \(x = 2\)’deki ileri fark bölümünün hatasının tam olarak \(6h + h^2\) olduğunu gösterin ve bunu \(h = 0{,}1;\ 0{,}01;\ 0{,}001\) için Python ile doğrulayın.
Çözüm
Hesap. \((2 + h)^3 = 8 + 12h + 6h^2 + h^3\) olduğundan
\[D_h^{+} f(2) = \frac{(2 + h)^3 - 8}{h} = 12 + 6h + h^2\]
olur. \(f'(2) = 3 \cdot 2^2 = 12\) olduğundan hata \(6h + h^2\)’dir. Bu, Önerme 11.1 ile uyumludur: \(f''(\xi) = 6\xi\) olduğundan \(\frac{h}{2} f''(\xi) = 3h\xi\)’dir ve \(\xi = 2 + h/3\) için \(3h\xi = 6h + h^2\) olur.
Kod.
def f(x):
return x**3
x = 2.0
for h in [0.1, 0.01, 0.001]:
err = (f(x + h) - f(x)) / h - 12.0 # f'(2) = 12
print(f"h = {h:5.3f} hata = {err:.10f}"
f" 6h + h^2 = {6 * h + h**2:.10f} hata/h = {err / h:.4f}")Çıktı:
h = 0.100 hata = 0.6100000000 6h + h^2 = 0.6100000000 hata/h = 6.1000
h = 0.010 hata = 0.0601000000 6h + h^2 = 0.0601000000 hata/h = 6.0100
h = 0.001 hata = 0.0060010000 6h + h^2 = 0.0060010000 hata/h = 6.0010Sonuç. Kayan noktalı hesap on ondalık basamağa kadar formülle aynı. Hata \(h\) ile orantılıdır ve hata bölü \(h\) oranı \(f''(2)/2 = 6\)’ya yaklaşır: ileri fark birinci mertebedendir.
\(\blacksquare\)
Alıştırma 11.2 (Beş Noktalı Merkezi Fark) Beş noktalı merkezi fark formülü
\[f'(x) \approx \frac{f(x - 2h) - 8f(x - h) + 8f(x + h) - f(x + 2h)}{12h}\]
ile verilir ve kesme hatası \(\frac{h^4}{30}\,|f^{(5)}(\xi)|\) kadardır. \(f(x) = \sin x\) için \(x = 1\)’de bu formülü \(h = 10^{-1}, \dots, 10^{-6}\) adımlarıyla deneyin. Hatanın \(h^4\) ile küçüldüğü adımları ve en iyi adımı bulun.
Çözüm
Toplam hata. Formüldeki katsayıların mutlak değerlerinin toplamı \(1 + 8 + 8 + 1 = 18\)’dir. Her fonksiyon değeri en çok \(\delta\) kadar yanlışsa yuvarlamadan gelen hata en çok \(\frac{18\delta}{12h} = \frac{3\delta}{2h}\) olur. Toplam hata sınırı \(E(h) = \frac{M_5 h^4}{30} + \frac{3\delta}{2h}\)’dir. \(E'(h) = \frac{4 M_5 h^3}{30} - \frac{3\delta}{2h^2} = 0\) denklemi \(h^5 = \frac{45\delta}{4M_5}\) verir. \(\delta = \varepsilon \sin 1\) ve \(M_5 = \cos 1\) alırız.
Kod.
import math
def five_point(f, x, h):
"""Beş noktalı merkezi fark formülü."""
return (f(x - 2 * h) - 8 * f(x - h)
+ 8 * f(x + h) - f(x + 2 * h)) / (12 * h)
x = 1.0
exact = math.cos(x)
old = None
for k in range(1, 7):
h = 10.0**-k
err = abs(five_point(math.sin, x, h) - exact)
ratio = f"{old / err:8.0f}" if old else " -"
print(f"h = {h:6.0e} hata = {err:.3e} oran = {ratio}")
old = err
eps = 2.0**-52
h_best = (45 * eps * math.sin(x) / (4 * math.cos(x))) ** (1 / 5)
print(f"teorik en iyi h = {h_best:.1e}")Çıktı:
h = 1e-01 hata = 1.799e-06 oran = -
h = 1e-02 hata = 1.801e-10 oran = 9988
h = 1e-03 hata = 1.497e-13 oran = 1203
h = 1e-04 hata = 5.385e-14 oran = 3
h = 1e-05 hata = 4.665e-12 oran = 0
h = 1e-06 hata = 1.847e-11 oran = 0
teorik en iyi h = 1.3e-03Sonuç. \(h = 10^{-1}\)’den \(10^{-2}\)’ye geçerken hata yaklaşık \(10^4\) kat küçülüyor: formül dördüncü mertebedendir. \(h = 10^{-1}\)’deki hata da \(\frac{h^4}{30}\cos 1 \approx 1{,}80 \cdot 10^{-6}\) tahminiyle uyuşuyor. \(h = 10^{-3}\)’te yuvarlama hatası öne geçmeye başlıyor ve oran düşüyor: bu adımda kesme hatası yalnız \(1{,}8 \cdot 10^{-14}\) kadardır, gerçek hata ise \(1{,}5 \cdot 10^{-13}\)’tür. Teorik en iyi adım \(1{,}3 \cdot 10^{-3}\)’tür. Gerçekten de en küçük hatalar \(10^{-3}\) ve \(10^{-4}\) adımlarında, \(10^{-13}\) mertebesinde çıktı. Merkezi farkla ulaşılabilen en iyi hata \(\varepsilon^{2/3} \approx 10^{-11}\) mertebesindeydi; dört fonksiyon değeri kullanan bu formül yaklaşık yüz kat daha iyi sonuç veriyor.
\(\blacksquare\)
Alıştırma 11.3 (Karmaşık Adımla Türev) \(f\), gerçel eksende gerçel değerler alan analitik bir fonksiyon olsun. Taylor açılımından \(f'(x) \approx \mathrm{Im}\, f(x + ih)/h\) yaklaşımının hatasının \(h^2\) mertebesinde olduğunu gösterin. Bu karmaşık adım (complex step) yaklaşımını \(f(x) = e^x \sin x\) için \(x = 1\)’de \(h = 10^{-2}\), \(10^{-5}\), \(10^{-10}\), \(10^{-20}\) ve \(10^{-100}\) adımlarıyla deneyin ve \(h\) çok küçükken bile sadeleşme hatası oluşmamasının nedenini açıklayın.
Çözüm
Taylor açılımı. \(i^2 = -1\) ve \(i^3 = -i\) olduğundan
\[f(x + ih) = f(x) + ih f'(x) - \frac{h^2}{2} f''(x) - i\,\frac{h^3}{6} f'''(x) + \cdots\]
olur. \(x\) gerçel ve \(f\)’nin türevleri gerçel olduğundan sanal kısım \(h f'(x) - \frac{h^3}{6} f'''(x) + \cdots\)’dır. \(h\)’ye bölünce
\[\frac{\mathrm{Im}\, f(x + ih)}{h} = f'(x) - \frac{h^2}{6} f'''(x) + \cdots\]
bulunur; hata \(h^2\) mertebesindedir.
Kod. np.exp ve np.sin karmaşık sayılarla da çalışır; math.exp ve math.sin çalışmaz. Kesin türev \(f'(x) = e^x(\sin x + \cos x)\)’tir.
import math
import numpy as np
def f(z):
return np.exp(z) * np.sin(z) # karmaşık sayılarla da çalışır
x = 1.0
exact = math.e * (math.sin(1) + math.cos(1))
for h in [1e-2, 1e-5, 1e-10, 1e-20, 1e-100]:
d = f(x + 1j * h).imag / h
print(f"h = {h:7.0e} türev = {d:.16f} hata = {abs(d - exact):.1e}")Çıktı:
h = 1e-02 türev = 3.7560765145542852 hata = 2.7e-05
h = 1e-05 türev = 3.7560492271220163 hata = 2.7e-11
h = 1e-10 türev = 3.7560492270947270 hata = 4.4e-16
h = 1e-20 türev = 3.7560492270947274 hata = 0.0e+00
h = 1e-100 türev = 3.7560492270947279 hata = 4.4e-16Sonuç. Formülde birbirine yakın iki sayının farkı yoktur: \(h f'(x)\) terimi sanal kısım olarak doğrudan hesaplanır. Bu yüzden yuvarlama hatası \(h\) küçüldükçe büyümez. \(h \le 10^{-10}\) için kesme hatası (\(h^2\) mertebesinde) makine duyarlığının altına düşer ve sonuç son basamağına kadar doğru çıkar; \(h = 10^{-100}\) bile sorun yaratmaz. Yöntemin bedeli, \(f\)’nin karmaşık sayılarla hesaplanabilmesinin gerekmesidir.
\(\blacksquare\)
Alıştırma 11.4 (Yuvarlanmış Ölçümlerden İvme) Yukarı atılan bir topun yüksekliği \(s(t) = 20t - 4{,}9t^2\) metre olsun. Yükseklik \([0, 4]\) saniye aralığında \(\Delta t\) aralıklarla ölçülüp santimetreye yuvarlanmış olsun. np.gradient’i iki kez uygulayarak hızı ve ivmeyi \(\Delta t = 0{,}1;\ 0{,}05;\ 0{,}01\) için hesaplayın ve en büyük hataları kesin değerler \(v(t) = 20 - 9{,}8t\) ve \(a = -9{,}8\) ile karşılaştırın.
Çözüm
Kesme hatası yok. \(s\) ikinci dereceden bir polinomdur; np.gradient (edge_order=2 ile) böyle bir fonksiyonun türevini tam verir. Hızın kendisi de birinci dereceden olduğundan ikinci uygulama da tamdır. Demek ki bütün hata yükseklik değerlerinin yuvarlanmasından gelir.
Kod.
import numpy as np
for dt in [0.1, 0.05, 0.01]:
t = np.arange(0, 4 + dt / 2, dt)
s = np.round(20 * t - 4.9 * t**2, 2) # santimetreye yuvarlanmış
v = np.gradient(s, t, edge_order=2)
a = np.gradient(v, t, edge_order=2)
err_v = np.abs(v - (20 - 9.8 * t)).max()
err_a = np.abs(a + 9.8).max()
print(f"dt = {dt:4.2f} hız hatası = {err_v:.3f}"
f" ivme hatası = {err_a:.2f}")Çıktı:
dt = 0.10 hız hatası = 0.040 ivme hatası = 0.45
dt = 0.05 hız hatası = 0.100 ivme hatası = 3.20
dt = 0.01 hız hatası = 0.464 ivme hatası = 40.20Sonuç. Ölçüm aralığı küçüldükçe hatalar büyüyor. Yuvarlama hatası en çok \(\delta = 0{,}005\) metredir. Hızda bu hata \(\delta/\Delta t\) ile, ivmede iki kez bölündüğü için \(\delta/\Delta t^2\) ile büyür: \(\Delta t\) on kat küçülünce ivme hatası yaklaşık yüz kat büyüdü (\(0{,}45\)’ten \(40{,}2\)’ye). Türev almak verideki hatayı büyütür. Gürültülü bir veriden türev alırken çözüm aralığı küçültmek değil, veriye önce düzgün bir eğri uydurmaktır (bkz. Polinomlar, İnterpolasyon ve Eğri Uydurma).
\(\blacksquare\)
Alıştırma 11.5 (Yamuk Kuralı İçin Yeterli Alt Aralık Sayısı) \(\int_0^1 e^x\,dx\) integralini bileşik yamuk kuralıyla \(10^{-6}\)’dan küçük hatayla hesaplamak için Önerme 11.3 sınırına göre kaç alt aralık yeter? Bu sayı için gerçek hatayı hesaplayın ve gerçekte yeten en küçük alt aralık sayısını Python ile bulun.
Çözüm
Sınır. \([0, 1]\)’de \(|f''(x)| = e^x \le e\) ve \(h = 1/n\) olduğundan hata en çok \(\frac{e}{12 n^2}\)’dir. \(\frac{e}{12 n^2} \le 10^{-6}\) eşitsizliği \(n \ge \sqrt{e/(12 \cdot 10^{-6})} \approx 475{,}9\) verir; sınır \(n = 476\)’yı önerir.
Kod. Kesin değer \(e - 1\)’dir. Gerçekte yeten en küçük \(n\)’yi \(n = 1\)’den başlayıp sayarak buluruz:
import math
import numpy as np
def trap_error(n):
x = np.linspace(0, 1, n + 1)
return abs(np.trapezoid(np.exp(x), x) - (math.e - 1))
tol = 1e-6
n_bound = math.ceil(math.sqrt(math.e / (12 * tol)))
print("sınırın verdiği n:", n_bound, " hata:", trap_error(n_bound))
n = 1
while trap_error(n) > tol:
n += 1
print("yeterli en küçük n:", n, " hata:", trap_error(n))
print("bir eksiği:", n - 1, " hata:", trap_error(n - 1))Çıktı:
sınırın verdiği n: 476 hata: 6.319740037952215e-07
yeterli en küçük n: 379 hata: 9.96861172941621e-07
bir eksiği: 378 hata: 1.0021425469464162e-06Sonuç. Sınırın önerdiği \(n = 476\)’da hata \(6{,}3 \cdot 10^{-7}\)’dir; gerçekte \(n = 379\) yeter. Sınır kötümserdir, çünkü \(f''\) yerine onun \([0, 1]\)’deki en büyük değeri \(e\)’yi kullanır. Oysa \(f''(\mu)\) değeri \(f''\)’nin bir ortalamasıdır ve \(e^x\)’in \([0, 1]\)’deki ortalaması \(e - 1 \approx 1{,}72\)’dir. Sınırdaki \(e\) yerine bu değer yazılırsa \(n \ge 378{,}4\) bulunur; bu da gerçekte yeten \(n = 379\) ile uyuşur.
\(\blacksquare\)
Alıştırma 11.6 (Simpson Kuralında Beklenmedik Hız) \(\int_0^1 \frac{4}{1 + x^2}\,dx = \pi\) integralini simpson ile \(n = 2, 4, 8, 16, 32, 64\) alt aralık için hesaplayın. Ardışık hataların oranına bakın ve beklenen \(16\) yerine neden başka bir oran çıktığını açıklayın.
Çözüm
Kod. Hataların yanında \(f'''\)’nün uç noktalardaki değerlerini SymPy ile hesaplıyoruz; nedenini birazdan göreceğiz:
import math
import numpy as np
import sympy as sp
from scipy.integrate import simpson
old = None
for n in [2, 4, 8, 16, 32, 64]:
x = np.linspace(0, 1, n + 1)
err = abs(simpson(4 / (1 + x**2), x=x) - math.pi)
ratio = f"{old / err:6.1f}" if old else " -"
print(f"n = {n:2d} hata = {err:.3e} oran = {ratio}")
old = err
t = sp.symbols("t")
d3 = sp.diff(4 / (1 + t**2), t, 3)
print(d3.subs(t, 0), d3.subs(t, 1))Çıktı:
n = 2 hata = 8.259e-03 oran = -
n = 4 hata = 2.403e-05 oran = 343.8
n = 8 hata = 1.511e-07 oran = 159.0
n = 16 hata = 2.365e-09 oran = 63.9
n = 32 hata = 3.696e-11 oran = 64.0
n = 64 hata = 5.773e-13 oran = 64.0
0 0Oran. Küçük \(n\)’lerden sonra oran \(16\) değil, \(64 = 2^6\) çıkıyor: hata \(h^6\) ile küçülüyor.
Açıklama. Önerme 11.3 hatayı \(-\frac{(b - a)h^4}{180} f^{(4)}(\mu)\) olarak verir. İspattaki toplam bir Riemann toplamı gibidir: \((b - a) f^{(4)}(\mu)\), \(\int_a^b f^{(4)}(x)\,dx = f'''(b) - f'''(a)\) değerine yaklaşır. Daha ayrıntılı bir analiz (Euler–Maclaurin formülü) hatanın
\[\int_a^b f(x)\,dx - S_n = -\frac{h^4}{180}\big[f'''(b) - f'''(a)\big] + O(h^6)\]
biçiminde açıldığını gösterir. Burada \(f'''(0) = f'''(1) = 0\) olduğu için \(h^4\) terimi yok olur ve geriye \(h^6\) mertebesindeki terim kalır.
\(\blacksquare\)
Alıştırma 11.7 (Periyodik Fonksiyonda Yamuk Kuralı) \(\int_0^{2\pi} e^{\cos x}\,dx = 2\pi I_0(1)\)’dir; burada \(I_0\), scipy.special.i0 ile hesaplanan değiştirilmiş Bessel fonksiyonudur. Bu integrali yamuk kuralıyla \(n = 2, 4, 8, 16, 32\) alt aralık için hesaplayın ve hatanın neden \(h^2\)’den çok daha hızlı küçüldüğünü açıklayın.
Çözüm
Kod.
import numpy as np
from scipy.special import i0
exact = 2 * np.pi * i0(1)
print(exact)
for n in [2, 4, 8, 16, 32]:
x = np.linspace(0, 2 * np.pi, n + 1)
err = abs(np.trapezoid(np.exp(np.cos(x)), x) - exact)
print(f"n = {n:2d} hata = {err:.1e}")Çıktı:
7.954926521012844
n = 2 hata = 1.7e+00
n = 4 hata = 3.4e-02
n = 8 hata = 1.3e-06
n = 16 hata = 0.0e+00
n = 32 hata = 0.0e+00Açıklama. Hata \(n = 4\)’ten \(n = 8\)’e geçerken yirmi binden fazla kat küçüldü; \(n = 16\)’da makine duyarlığına ulaştı. Simpson kuralındaki akıl yürütmenin yamuk kuralı için karşılığı, hatanın
\[\int_a^b f(x)\,dx - T_n = -\frac{h^2}{12}\big[f'(b) - f'(a)\big] + O(h^4)\]
biçiminde açılmasıdır; sonraki terimler de \(f'''\), \(f^{(5)}, \dots\) türevlerinin uç noktalardaki değerlerinin farklarını içerir. \(e^{\cos x}\) fonksiyonu \(2\pi\) periyotlu ve sonsuz kez türevlenebilir olduğundan bütün türevleri \(0\) ile \(2\pi\)’de aynı değeri alır. Bu yüzden açılımın bütün terimleri yok olur ve hata \(h\)’nin her kuvvetinden hızlı küçülür. Periyodik fonksiyonların tam periyot üzerindeki integralinde en iyi kural, en basit kural olan yamuk kuralıdır.
\(\blacksquare\)
Alıştırma 11.8 (Sinüs Eğrisinin Uzunluğu) \(y = \sin x\) eğrisinin \(0 \le x \le \pi\) arasındaki yayının uzunluğunu
\[L = \int_0^\pi \sqrt{1 + \cos^2 x}\,dx\]
formülüyle quad kullanarak hesaplayın. Sonucu, eğri üzerindeki \(n + 1\) eşit aralıklı noktayı birleştiren kırık çizgilerin \(n = 4, 16, 64, 256\) için uzunluklarıyla karşılaştırın.
Çözüm
Kırık çizgi. Ardışık noktalar arasındaki kirişin boyu \(\sqrt{(\Delta x)^2 + (\Delta y)^2}\)’dir; np.diff farkları, np.hypot bu karekökü verir.
import numpy as np
from scipy.integrate import quad
L, err = quad(lambda x: np.sqrt(1 + np.cos(x)**2), 0, np.pi)
print(L, err)
for n in [4, 16, 64, 256]:
x = np.linspace(0, np.pi, n + 1)
chords = np.hypot(np.diff(x), np.diff(np.sin(x))) # kiriş boyları
print(f"n = {n:3d} kırık çizgi = {chords.sum():.10f}"
f" fark = {L - chords.sum():.1e}")Çıktı:
3.8201977890277115 1.3015768476030976e-13
n = 4 kırık çizgi = 3.7900913085 fark = 3.0e-02
n = 16 kırık çizgi = 3.8182749650 fark = 1.9e-03
n = 64 kırık çizgi = 3.8200775044 fark = 1.2e-04
n = 256 kırık çizgi = 3.8201902708 fark = 7.5e-06Sonuç. İntegral elementer fonksiyonlarla hesaplanamaz (bir eliptik integraldir); quad onu \(10^{-13}\) hata tahminiyle \(L \approx 3{,}8201977890\) olarak verdi. Kırık çizgilerin uzunlukları bu değere alttan yaklaşır, çünkü her kiriş, gerdiği yaydan kısadır. \(n\) dört katına çıkınca fark yaklaşık 16 kat küçülüyor; kırık çizginin hatası \(h^2\) ile orantılıdır.
\(\blacksquare\)
Alıştırma 11.9 (Normal Dağılımda Bir, İki ve Üç Standart Sapma) Standart normal dağılımın yoğunluğu \(\varphi(x) = e^{-x^2/2}/\sqrt{2\pi}\)’dir (bkz. Olasılık Teorisi). \(Z\) standart normal dağılımlı olmak üzere \(k = 1, 2, 3\) için \(P(|Z| \le k) = \int_{-k}^{k} \varphi(x)\,dx\) olasılıklarını quad ile hesaplayın ve \(\mathrm{erf}(k/\sqrt{2})\) değerleriyle karşılaştırın.
Çözüm
Denetim. Önce \(\varphi\)’nin bütün doğru üzerindeki integralinin \(1\) olduğunu, yani \(\varphi\)’nin gerçekten bir yoğunluk olduğunu doğrularız. Sonra her \(k\) için integrali alırız. \(x = \sqrt{2}\,u\) değişken değiştirmesi \(\int_{-k}^{k} \varphi(x)\,dx = \mathrm{erf}(k/\sqrt{2})\) verir.
import math
import numpy as np
from scipy.integrate import quad
def phi(x):
return np.exp(-x**2 / 2) / np.sqrt(2 * np.pi)
print(quad(phi, -np.inf, np.inf))
for k in [1, 2, 3]:
p, err = quad(phi, -k, k)
print(f"k = {k} olasılık = {p:.10f}"
f" erf = {math.erf(k / math.sqrt(2)):.10f}")Çıktı:
(0.9999999999999998, 1.0178191320905743e-08)
k = 1 olasılık = 0.6826894921 erf = 0.6826894921
k = 2 olasılık = 0.9544997361 erf = 0.9544997361
k = 3 olasılık = 0.9973002039 erf = 0.9973002039Sonuç. Olasılıklar yaklaşık \(0{,}6827\), \(0{,}9545\) ve \(0{,}9973\)’tür ve on ondalık basamağa kadar math.erf ile aynıdır. İstatistikteki “68–95–99,7 kuralı” bu üç sayıdır.
\(\blacksquare\)
Alıştırma 11.10 (Paraboloidin Altındaki Hacim) \(z = 1 - x^2 - y^2\) paraboloidi ile \(z = 0\) düzlemi arasındaki cismin hacmini dblquad ile önce dik koordinatlarda, sonra kutupsal koordinatlarda hesaplayın ve sonuçları \(\pi/2\) ile karşılaştırın (bkz. İntegral Calculus).
Çözüm
Sınırlar. Taban birim disktir. Dik koordinatlarda \(-1 \le x \le 1\) ve \(-\sqrt{1 - x^2} \le y \le \sqrt{1 - x^2}\)’dir. Kutupsal koordinatlarda \(0 \le \theta \le 2\pi\), \(0 \le r \le 1\)’dir ve \(dA = r\,dr\,d\theta\) olduğundan integrand \((1 - r^2)\,r\) olur. dblquad’da iç değişken önce gelir: dik koordinatlarda func(y, x), kutupsal koordinatlarda func(r, theta).
import numpy as np
from scipy.integrate import dblquad
# Dik koordinatlar: -1 <= x <= 1, -sqrt(1 - x^2) <= y <= sqrt(1 - x^2)
v1, e1 = dblquad(lambda y, x: 1 - x**2 - y**2, -1, 1,
lambda x: -np.sqrt(1 - x**2),
lambda x: np.sqrt(1 - x**2))
# Kutupsal koordinatlar: 0 <= theta <= 2 pi, 0 <= r <= 1, dA = r dr dtheta
v2, e2 = dblquad(lambda r, theta: (1 - r**2) * r, 0, 2 * np.pi, 0, 1)
print(v1, e1)
print(v2, e2)
print(np.pi / 2)Çıktı:
1.5707963267740037 2.2114190221070963e-08
1.5707963267948966 1.743934249004316e-14
1.5707963267948966Sonuç. Kutupsal koordinatlardaki sonuç son basamağına kadar \(\pi/2\)’dir; dik koordinatlardaki sonuç ise \(2 \cdot 10^{-11}\) kadar sapıyor. Fark integrandın düzgünlüğünden gelir. Kutupsal koordinatlarda iki integral de polinom integralidir. Dik koordinatlarda iç integral \(\frac{4}{3}(1 - x^2)^{3/2}\) fonksiyonunu verir; bu fonksiyonun \(x = \pm 1\)’de ikinci türevi sınırsızdır ve dış integrali zorlaştırır. İki sonuç da varsayılan \(1{,}5 \cdot 10^{-8}\) hedefinin çok altındadır.
\(\blacksquare\)
Alıştırma 11.11 (Bir Dörtyüzlünün Hacmi) \(x + 2y + z = 2\), \(x = 2y\), \(x = 0\) ve \(z = 0\) düzlemleriyle sınırlı dörtyüzlünün hacmini tplquad ile hesaplayın (bkz. İntegral Calculus).
Çözüm
Sınırlar. Dörtyüzlü alttan \(z = 0\), üstten \(z = 2 - x - 2y\) düzlemiyle sınırlıdır. Üst düzlem \(z = 0\) düzlemini \(y = 1 - x/2\) doğrusu boyunca keser. Bu doğru \(y = x/2\) doğrusuyla \(x = 1\)’de buluşur. Dolayısıyla cisim
\[0 \le x \le 1, \qquad \frac{x}{2} \le y \le 1 - \frac{x}{2}, \qquad 0 \le z \le 2 - x - 2y\]
eşitsizlikleriyle betimlenir. Hacim, \(1\) fonksiyonunun bu cisim üzerindeki integralidir.
from scipy.integrate import tplquad
value, err = tplquad(lambda z, y, x: 1.0,
0, 1, # x
lambda x: x / 2, lambda x: 1 - x / 2, # y
0, lambda x, y: 2 - x - 2 * y) # z
print(value, err, 1 / 3)Çıktı:
0.3333333333333333 2.210813483580885e-14 0.3333333333333333Sonuç. Hacim \(\frac{1}{3}\)’tür; İntegral Calculus notlarındaki elle hesapla aynı.
\(\blacksquare\)
Bu bölümde türev ve integrali fonksiyonun sonlu sayıda noktadaki değerlerinden hesapladık. Fark bölümlerinde kesme hatası ile yuvarlama hatasının dengesini, yamuk ve Simpson kurallarında mertebenin anlamını, quad ile de uyarlamalı yöntemlerin gücünü ve sınırlarını gördük. Bu formüllerin çoğu aslında bir interpolasyon polinomunun türevi ya da integralidir. Sonraki bölümde, Polinomlar, İnterpolasyon ve Eğri Uydurma, verilen noktalardan geçen ya da onlara en iyi uyan fonksiyonları kuracağız.