12 Polinomlar, İnterpolasyon ve Eğri Uydurma
Polinomlar bilgisayarın en rahat hesapladığı fonksiyonlardır: değerleri yalnız toplama ve çarpmayla bulunur, türevleri ve integralleri yine polinomdur. Bu yüzden sayısal hesabın büyük bir kısmı, karmaşık bir fonksiyonu ya da yalnız birkaç noktada bilinen bir büyüklüğü bir polinomla veya polinom parçalarıyla değiştirme fikrine dayanır. Sayısal Türev ve İntegral bölümündeki yamuk ve Simpson kuralları da aslında bir interpolasyon polinomunun integralidir.
Elimizdeki veriyle iki farklı soru sorabiliriz. Değerler kesinse (bir tablodan okunmuş ya da pahalı bir hesaptan gelmişse) verilen noktaların hepsinden geçen bir fonksiyon isteriz; buna interpolasyon denir. Değerler ölçümden geliyor ve gürültü içeriyorsa bütün noktalardan geçmeye çalışmak gürültüyü de kopyalamak olur. O zaman noktaların gidişine en iyi uyan basit bir eğri ararız; buna da eğri uydurma denir (bkz. Nümerik Analiz).
Bu bölümde önce NumPy’ın Polynomial sınıfıyla polinomlar üzerinde hesap yapıyor ve köklerin katsayılara ne kadar duyarlı olabildiğini görüyoruz. Ardından interpolasyon polinomunu üç yoldan kuruyor, eşit aralıklı noktalarda derece artınca ortaya çıkan Runge olayını ve onu gideren Çebişev düğümlerini inceliyoruz. Sonra SciPy’ın kübik spline’larına geçiyoruz. En sonda polyfit ve curve_fit ile ölçüm verisine model uyduruyor ve artıklara bakarak modelin yeterli olup olmadığına karar vermeyi öğreniyoruz. Yöntemlerin teorisi Nümerik Analiz notlarındadır; burada onların Python’daki karşılığını kuruyoruz.
12.1 NumPy’da Polinomlar
Polinomlarla çalışmak için NumPy’ın numpy.polynomial paketini kullanacağız.
Tanım 12.1 (Polynomial Nesnesi) numpy.polynomial.Polynomial sınıfının bir örneğine Polynomial nesnesi denir. Polynomial([a0, a1, ..., an]) çağrısı
\[ p(x) = a_0 + a_1 x + a_2 x^2 + \dots + a_n x^n \]
polinomunu temsil eden bir nesne kurar. Katsayılar artan derece sırasıyla verilir ve nesnenin coef niteliğinde saklanır.
Yani bir Polynomial nesnesi, katsayı listesini taşıyan ve bir polinom gibi davranan bir Python nesnesidir: fonksiyon gibi çağrılır, başka polinomlarla toplanıp çarpılır, türevi ve integrali alınır. Sınıflar, Hata Yakalama ve Dosyalar bölümünün alıştırmalarında yazdığımız Poly sınıfı bu fikrin küçük bir örneğiydi. \(p(x) = (x-1)(x-2)(x-3)\) polinomuyla başlayalım. Açılmış hâli \(-6 + 11x - 6x^2 + x^3\) olduğundan katsayı listesi [-6, 11, -6, 1]’dir:
import numpy as np
from numpy.polynomial import Polynomial
p = Polynomial([-6, 11, -6, 1]) # -6 + 11x - 6x^2 + x^3
print(p)
print(p.coef, p.degree())
print(p(4))
print(p(np.array([0.0, 1.0, 2.5])))Çıktı:
-6.0 + 11.0 x - 6.0 x**2 + 1.0 x**3
[-6. 11. -6. 1.] 3
6.0
[-6. 0. -0.375]degree metodu polinomun derecesini verir. Nesne bir sayıyla çağrılınca bir sayı, bir NumPy dizisiyle çağrılınca dizinin her elemanı için ayrı bir değer döndürür; örneğin \(p(2{,}5) = 1{,}5 \cdot 0{,}5 \cdot (-0{,}5) = -0{,}375\)’tir.
Polinomun yazdırılış biçimi işletim sistemine bağlıdır. Bu notlardaki çıktılar Windows’ta alındı ve kuvvetleri x**2 biçiminde gösteriyor. Linux, macOS ve Colab’da aynı polinom -6.0 + 11.0·x - 6.0·x² + 1.0·x³ biçiminde yazılır. İsterseniz programın başında np.polynomial.set_default_printstyle("ascii") çağırarak biçimi her yerde aynı yapabilirsiniz.
Polinomlar arasındaki işlemler alışılmış operatörlerle yazılır. divmod(p, q) bölümü ve kalanı birlikte verir; deriv türevi, integ ise sabit terimi sıfır olan ilkel fonksiyonu döndürür:
from numpy.polynomial import Polynomial
p = Polynomial([-6, 11, -6, 1])
q = Polynomial([1, 1]) # 1 + x
print(p + q)
print(p * q)
print(q**3)
quot, rem = divmod(p, q)
print(quot, "|", rem)
print(p.deriv())
P = p.integ() # P(0) = 0 olan ilkel fonksiyon
print(P)
print(P(1) - P(0))Çıktı:
-5.0 + 12.0 x - 6.0 x**2 + 1.0 x**3
-6.0 + 5.0 x + 5.0 x**2 - 5.0 x**3 + 1.0 x**4
1.0 + 3.0 x + 3.0 x**2 + 1.0 x**3
18.0 - 7.0 x + 1.0 x**2 | -24.0
11.0 - 12.0 x + 3.0 x**2
0.0 - 6.0 x + 5.5 x**2 - 2.0 x**3 + 0.25 x**4
-2.25Bölme satırı \(p(x) = (x + 1)(x^2 - 7x + 18) - 24\) eşitliğini söylüyor. Kalan teoremine göre kalan \(p(-1) = (-2)(-3)(-4) = -24\) olmalıdır ve öyle. Son satır, ilkel fonksiyonun uç noktalardaki farkı olarak
\[ \int_0^1 (x-1)(x-2)(x-3)\,dx = -\frac{9}{4} \]
integralini hesapladı.
Bir polinomun değerini hesaplamanın en ekonomik yolu, onu iç içe çarpımlar biçiminde yazmaktır.
Tanım 12.2 (Horner Yöntemi) \(p(x) = a_0 + a_1 x + \dots + a_n x^n\) polinomunun bir \(x\) noktasındaki değerini
\[ b_n = a_n, \qquad b_k = a_k + x\,b_{k+1} \quad (k = n-1, \dots, 1, 0) \]
yinelemesiyle hesaplayıp \(p(x) = b_0\) almaya Horner yöntemi denir.
Yani polinom
\[ p(x) = a_0 + x\Big(a_1 + x\big(a_2 + \dots + x\,(a_{n-1} + x\,a_n)\cdots\big)\Big) \]
biçiminde içten dışa doğru hesaplanır. Her adımda bir çarpma ve bir toplama yapıldığından toplam maliyet \(n\) çarpma ve \(n\) toplamadır; her \(a_k x^k\) terimini ayrı ayrı hesaplamak ise yaklaşık \(n^2/2\) çarpma ister. Yöntem birkaç satırda yazılır. NumPy işlemleri dizilere eleman eleman uygulandığından aynı kod bütün bir diziyle de çalışır:
import numpy as np
from numpy.polynomial import Polynomial
def horner(coef, x):
"""Katsayıları artan sırada verilen polinomun x'teki değeri."""
result = 0.0
for a in reversed(coef):
result = result * x + a
return result
p = Polynomial([-6, 11, -6, 1])
x = np.array([0.0, 1.0, 2.5, 4.0])
print(horner(p.coef, x))
print(p(x))Çıktı:
[-6. 0. -0.375 6. ]
[-6. 0. -0.375 6. ]Polynomial nesnesinin değer hesabı da aynı sonucu veriyor.
NumPy’da polinomlar için iki arayüz vardır. Eski arayüzün fonksiyonları (np.polyval, np.polyfit, np.roots, np.poly1d) katsayıları en yüksek dereceden başlayarak ister; Polynomial sınıfı ise artan sırada. Aynı listeyi iki arayüze vermek hiçbir hata mesajı üretmez, yalnız yanlış sonuç verir:
import numpy as np
from numpy.polynomial import Polynomial
c = [1, -6, 11, -6] # x^3 - 6x^2 + 11x - 6, büyük dereceden başlar
print(np.polyval(c, 4))
print(np.roots(c))
print(Polynomial(c)(4)) # yanlış: Polynomial c'yi artan sırada okur
print(Polynomial(c[::-1])(4))Çıktı:
6
[3. 2. 1.]
-231.0
6.0Polynomial(c) listeyi \(1 - 6x + 11x^2 - 6x^3\) olarak okudu ve \(x = 4\)’te \(-231\) verdi. Yeni kodda Polynomial sınıfını kullanın; eski fonksiyonlarla çalışan bir kodla karşılaşınca listeyi c[::-1] ile ters çevirin.
12.2 Kökler
Bir Polynomial nesnesinin köklerini roots metodu verir; Polynomial.fromroots ise tersini yaparak verilen köklerden baş katsayısı \(1\) olan polinomu kurar.
from numpy.polynomial import Polynomial
p = Polynomial([-6, 11, -6, 1])
print(p.roots())
q = Polynomial([5, 0, 1]) # 5 + x^2
print(q.roots())
r = Polynomial.fromroots([1, 1, -2]) # (x - 1)^2 (x + 2)
print(r)
print(r.roots())Çıktı:
[1. 2. 3.]
[0.-2.23606798j 0.+2.23606798j]
2.0 - 3.0 x + 0.0 x**2 + 1.0 x**3
[-2. 0.99999998 1.00000002]\(x^2 + 5\) polinomunun kökleri \(\pm i\sqrt{5} \approx \pm 2{,}23607\,i\)’dir; kökler karmaşık olunca sonuç da karmaşık sayılardan oluşan bir dizidir ve Python sanal birimi j ile yazar. Son satıra dikkat: \((x-1)^2(x+2)\) polinomunun çift katlı kökü \(1\), \(0{,}99999998\) ve \(1{,}00000002\) olarak bulundu. Hata \(10^{-8}\) mertebesindedir, yani makine hassasiyeti \(2{,}2 \cdot 10^{-16}\)’nın kareköküne yakındır. Katlı köklerin bu duyarlılığını alıştırmalarda yeniden göreceğiz (Alıştırma 12.3).
Beşinci ve daha yüksek dereceden polinomlar için genel bir kök formülü olmadığından kökleri bulan yöntem sayısal olmak zorundadır. NumPy’ın yolu, problemi bir özdeğer problemine çevirmektir.
Tanım 12.3 (Eşlik Matrisi) Baş katsayısı \(1\) olan
\[ p(x) = a_0 + a_1x + \dots + a_{n-1}x^{n-1} + x^n \]
polinomu verilsin. Alt köşegeni birlerden, son sütunu \(-a_0, -a_1, \dots, -a_{n-1}\) sayılarından oluşan ve öteki elemanları sıfır olan
\[ C = \begin{pmatrix} 0 & 0 & \cdots & 0 & -a_0 \\ 1 & 0 & \cdots & 0 & -a_1 \\ 0 & 1 & \cdots & 0 & -a_2 \\ \vdots & & \ddots & & \vdots \\ 0 & 0 & \cdots & 1 & -a_{n-1} \end{pmatrix} \]
\(n \times n\) matrisine \(p\)’nin eşlik matrisi (companion matrix) denir.
Yani eşlik matrisi, polinomun katsayılarını bir matrise yerleştirmenin standart bir yoludur. Önemi, özdeğerlerinin (bkz. Lineer Cebir) polinomun kökleri olmasından gelir.
Önerme 12.1 (Eşlik Matrisinin Özdeğerleri) \(r\), baş katsayısı \(1\) olan \(p\) polinomunun bir kökü ise \(r\), \(p\)’nin eşlik matrisi \(C\)’nin bir özdeğeridir. Özel olarak \(p\)’nin \(n\) farklı kökü varsa \(C\)’nin özdeğerleri tam olarak bu köklerdir.
İspat
\(v = (1, r, r^2, \dots, r^{n-1})\) satır vektörünü alıp \(vC\) çarpımını sütun sütun hesaplayalım. \(j < n\) için \(C\)’nin \(j\)’inci sütununda yalnız \((j+1)\)’inci satırda bir \(1\) vardır; bu yüzden \(vC\)’nin \(j\)’inci bileşeni \(v\)’nin \((j+1)\)’inci bileşenidir, yani \(r^{j} = r \cdot r^{j-1}\)’dir. Son sütunda \(-a_0, \dots, -a_{n-1}\) sayıları bulunur ve \(p(r) = 0\) olduğundan
\[ -\left(a_0 + a_1 r + \dots + a_{n-1} r^{n-1}\right) = r^n = r \cdot r^{n-1} \]
olur. Demek ki \(vC = r\,v\)’dir, yani \(C^T v^T = r\,v^T\). İlk bileşeni \(1\) olduğundan \(v \ne 0\)’dır; dolayısıyla \(r\), \(C^T\)’nin bir özdeğeridir. \(\det(C^T - \lambda I) = \det(C - \lambda I)\) olduğundan \(C\) ile \(C^T\)’nin özdeğerleri aynıdır ve \(r\), \(C\)’nin de bir özdeğeridir. \(n \times n\) bir matrisin en çok \(n\) özdeğeri olduğundan \(p\)’nin \(n\) farklı kökü varsa \(C\)’nin özdeğerleri tam olarak bunlardır. \(\blacksquare\)
Aslında daha fazlası doğrudur: \(C\)’nin karakteristik polinomu (bkz. Lineer Cebir) tam olarak \(p\)’dir; böylece katlı kökler de katlılıklarıyla birlikte özdeğer olur. Bunu \((x-1)(x-2)(x-3)\) için hem sayısal hem sembolik olarak görelim. polycompanion fonksiyonu, NumPy’ın kullandığı eşlik matrisini verir:
import numpy as np
import sympy as sp
from numpy.polynomial import polynomial as P
c = np.array([-6.0, 11.0, -6.0, 1.0]) # (x - 1)(x - 2)(x - 3)
C = P.polycompanion(c)
print(C)
print(np.linalg.eigvals(C))
x = sp.symbols("x")
print(sp.Matrix(C.astype(int)).charpoly(x).as_expr())Çıktı:
[[ 0. 0. 6.]
[ 1. 0. -11.]
[ 0. 1. 6.]]
[1.+0.j 2.+0.j 3.+0.j]
x**3 - 6*x**2 + 11*x - 6Matris tanımdaki biçimdedir: son sütunda \(6\), \(-11\), \(6\) sayıları \(-a_0\), \(-a_1\), \(-a_2\)’dir. np.linalg.eigvals özdeğerleri karmaşık sayı olarak döndürdü, ama sanal kısımları sıfırdır. SymPy’ın charpoly metodu (SymPy ile Sembolik Hesap) karakteristik polinomun \(p\)’nin kendisi olduğunu gösteriyor. roots metodu da bu matrisin özdeğerlerini np.linalg.eigvals ile hesaplar.
12.3 Köklerin Duyarlılığı
Özdeğer hesabı çok güvenilir bir algoritmadır, ama bu, katsayılardan hesaplanan köklerin her zaman doğru olduğu anlamına gelmez. Kökleri \(1, 2, \dots, 20\) olan
\[ w(x) = (x - 1)(x - 2) \cdots (x - 20) \]
Wilkinson polinomuna bakalım. Kökleri tam sayıdır ve birbirinden iyi ayrılmıştır; hesaplaması en kolay görünen polinomlardan biri olmalı.
import numpy as np
from numpy.polynomial import Polynomial
w = Polynomial.fromroots(np.arange(1, 21)) # (x - 1)(x - 2)...(x - 20)
print(w.coef[19], w.coef[0])
r = w.roots()
err = np.abs(r - np.arange(1, 21))
print(np.round(r[10:16], 4))
print(f"en büyük hata {err.max():.3f}, kök {np.argmax(err) + 1}")Çıktı:
-210.0 2.43290200817664e+18
[11.0485 11.9087 13.1831 13.774 15.2444 15.792 ]
en büyük hata 0.244, kök 15\(x^{19}\)’un katsayısı \(-210\), sabit terim \(20! \approx 2{,}43 \cdot 10^{18}\)’dir. Yazdırılan kökler \(11\)’den \(16\)’ya kadar olanlardır: \(15\) yerine \(15{,}2444\) bulundu ve en büyük hata \(0{,}244\). Kodda bir yanlışlık yok. Bulunan kökler, katsayıları \(w\)’nin katsayılarından yuvarlama hatası düzeyinde farklı olan bir polinomun neredeyse tam kökleridir; sorun \(w\)’nin köklerinin katsayılara aşırı duyarlı olmasıdır. Bunu James Wilkinson’ın ünlü deneyiyle görelim: yalnız \(x^{19}\)’un katsayısını \(2^{-23} \approx 1{,}2 \cdot 10^{-7}\) kadar azaltalım.
import numpy as np
from numpy.polynomial import Polynomial
w = Polynomial.fromroots(np.arange(1, 21))
c = w.coef.copy()
c[19] -= 2.0**-23 # x^19 katsayısı: -210 - 2^(-23)
r = Polynomial(c).roots()
r = r[np.argsort(r.real)]
print("karmaşık kök sayısı:", np.sum(r.imag != 0))
for z in r[r.imag >= 0][6:]: # eşlenik çiftlerden yalnız biri
print(f"{z.real:8.4f} {z.imag:+8.4f}i")Çıktı:
karmaşık kök sayısı: 10
6.9997 +0.0000i
8.0075 +0.0000i
8.9153 +0.0000i
10.0943 +0.6478i
11.7945 +1.6542i
13.9930 +2.5192i
16.7310 +2.8127i
19.5025 +1.9403i
20.8469 +0.0000iYirmi kökün onu karmaşık oldu; eşlenik çiftlerden yalnız üst yarı düzlemdekileri yazdırdık (ilk altı kök \(1, \dots, 6\)’ya çok yakın olduğu için atlandı). Bulunan köklerin hepsi karmaşık düzlemde şöyle görünüyor:
Katsayıdaki bağıl değişiklik \(2^{-23}/210 \approx 5{,}7 \cdot 10^{-10}\) iken kökler \(2{,}8\) birime kadar yer değiştirdi. Köklerin katsayılardan hesaplanması bu yüzden kötü koşullu bir problem olabilir: girdideki çok küçük bir değişiklik çıktıyı çok değiştirir. Aynı kavramı lineer sistemler için NumPy ile Lineer Cebir bölümünde koşul sayısıyla ölçmüştük. Bir polinom katsayılarıyla değil de kökleriyle ya da değerleriyle biliniyorsa, onu katsayılara çevirmeden çalışmak çoğu zaman daha güvenlidir. Bölümün geri kalanında interpolasyonda da aynı ilkeyi izleyeceğiz.
12.4 İnterpolasyon Polinomu
Şimdi bölümün asıl sorusuna geçelim: verilen noktalardan geçen polinomu bulmak.
Teorem 12.1 (İnterpolasyon Polinomunun Varlığı ve Tekliği) \(x_0, x_1, \dots, x_n\) birbirinden farklı sayılar ve \(y_0, y_1, \dots, y_n\) keyfi sayılar olsun. Derecesi en fazla \(n\) olan ve \(k = 0, 1, \dots, n\) için \(p(x_k) = y_k\) koşulunu sağlayan bir ve yalnız bir \(p\) polinomu vardır.
İspatı için bkz. Nümerik Analiz. Bu polinoma verinin interpolasyon polinomu, \(x_k\) noktalarına da düğümler denir.
Polinomu bulmanın en doğrudan yolu, bilinmeyen katsayıları bir lineer sistemin çözümü olarak yazmaktır. \(p(x) = a_0 + a_1x + \dots + a_nx^n\) için \(p(x_k) = y_k\) koşulları
\[ \begin{pmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^n \\ 1 & x_1 & x_1^2 & \cdots & x_1^n \\ \vdots & \vdots & \vdots & & \vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^n \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ \vdots \\ a_n \end{pmatrix} = \begin{pmatrix} y_0 \\ y_1 \\ \vdots \\ y_n \end{pmatrix} \]
sistemini verir. Katsayı matrisi, NumPy Dizileri bölümünde np.vander ile kurduğumuz Vandermonde matrisidir ve teorem, düğümler farklıyken bu sistemin tek bir çözümü olduğunu söyler. \((0, 1)\), \((1, 3)\), \((2, 2)\) ve \((3, 5)\) noktalarından geçen kübik polinomu bulalım:
import numpy as np
from numpy.polynomial import Polynomial
x = np.array([0.0, 1.0, 2.0, 3.0])
y = np.array([1.0, 3.0, 2.0, 5.0])
V = np.vander(x, increasing=True) # sütunlar: 1, x, x^2, x^3
a = np.linalg.solve(V, y)
p = Polynomial(a)
print(p)
print(p(x))
print(p(1.5))Çıktı:
1.0 + 5.83333333 x - 5.0 x**2 + 1.16666667 x**3
[1. 3. 2. 5.]
2.437500000000001Polinom dört düğümde de verilen değerleri alıyor. Kesin katsayılar \(a_1 = \frac{35}{6}\) ve \(a_3 = \frac{7}{6}\)’dır, yani
\[ p(x) = 1 + \frac{35}{6}\,x - 5x^2 + \frac{7}{6}\,x^3 . \]
\(p(1{,}5) = \frac{39}{16} = 2{,}4375\) değeri son basamakta küçük bir yuvarlama hatasıyla çıktı.
Aynı polinomu Polynomial.fit(x, y, deg) de verir. Bu yöntem asıl olarak en küçük kareler uydurması için yazılmıştır ve onu bölümün sonunda kullanacağız. Nokta sayısı deg + 1 olduğunda interpolasyon polinomu artıkları sıfır yaptığından en küçük kareler çözümü tam olarak interpolasyon polinomudur.
import numpy as np
from numpy.polynomial import Polynomial
x = np.array([0.0, 1.0, 2.0, 3.0])
y = np.array([1.0, 3.0, 2.0, 5.0])
p = Polynomial.fit(x, y, deg=3)
print(p)
print(p.domain, p.window)
print(p.convert())
print(p(1.5))Çıktı:
2.4375 - 1.9375 (-1.0 + 0.66666667x) + 0.5625 (-1.0 + 0.66666667x)**2 +
3.9375 (-1.0 + 0.66666667x)**3
[0. 3.] [-1. 1.]
1.0 + 5.83333333 x - 5.0 x**2 + 1.16666667 x**3
2.4374999999999996Çıktının ilk iki satırı şaşırtıcı görünebilir. fit, verinin bulunduğu \([0, 3]\) aralığını (domain) \([-1, 1]\) aralığına (window) taşıyan \(u = -1 + \frac{2}{3}x\) dönüşümünü yapar ve katsayıları \(u\)’nun kuvvetlerine göre saklar. Yazdırılan polinom
\[ 2{,}4375 - 1{,}9375u + 0{,}5625u^2 + 3{,}9375u^3 \]
polinomudur. \(u = 0\) noktası \(x = 1{,}5\)’e karşılık geldiğinden sabit terim \(p(1{,}5)\)’tir. Değer hesapları doğru çıkar, çünkü nesne dönüşümü kendisi yapar; ama p.coef alışılmış katsayılar değildir. Alışılmış katsayılar için convert() çağrılmalıdır. Ölçeklemenin nedeni sayısaldır: \([-1, 1]\) aralığında \(u\)’nun hiçbir kuvveti \(1\)’i aşmaz ve Vandermonde matrisi çok daha iyi koşullu olur (bkz. NumPy ile Lineer Cebir bölümündeki Vandermonde alıştırması).
Lineer sistem çözmeden yazılan bir biçim de vardır. Lagrange katsayı polinomları
\[ L_k(x) = \prod_{j \ne k} \frac{x - x_j}{x_k - x_j}, \qquad k = 0, 1, \dots, n \]
ile interpolasyon polinomu \(p(x) = \sum_{k=0}^{n} y_k L_k(x)\) biçiminde yazılır. \(L_k\) polinomu \(x_k\)’de \(1\), öteki düğümlerde \(0\) değerini alır; bu yüzden toplam her düğümde doğru değeri verir. \(L_0, \dots, L_n\) polinomları, derecesi en fazla \(n\) olan polinomların uzayı için bir taban oluşturduğundan bunlara Lagrange taban polinomları da denir. Formülü broadcasting ile vektörel olarak yazabiliriz (Vektörizasyon ve Broadcasting): np.delete(xn, k), xn’den k indisli elemanı çıkarılmış yeni bir dizi döndürür; x[:, None] - others farkı, her değerlendirme noktası için bir satır, her öteki düğüm için bir sütun içeren bir matristir ve np.prod(..., axis=1) her satırın çarpımını alır.
import numpy as np
def lagrange_basis(xn, k, x):
"""k'inci Lagrange taban polinomunun x noktalarındaki değerleri."""
others = np.delete(xn, k) # x_j, j != k
return np.prod((x[:, None] - others) / (xn[k] - others), axis=1)
def lagrange_eval(xn, yn, x):
"""Lagrange biçimindeki interpolasyon polinomunun değerleri."""
return sum(yn[k] * lagrange_basis(xn, k, x) for k in range(len(xn)))
xn = np.array([0.0, 1.0, 2.0, 3.0])
yn = np.array([1.0, 3.0, 2.0, 5.0])
print(lagrange_basis(xn, 1, xn)) # düğümlerde 0, 1, 0, 0
print(lagrange_eval(xn, yn, np.array([1.5, 2.5])))Çıktı:
[ 0. 1. -0. 0.]
[2.4375 2.5625]İlk satır \(L_1\)’in dört düğümdeki değerleridir. Oradaki -0. eksi işaretli sıfırdır: kayan noktalı sayılar sıfırın işaretini de saklar, ama değeri yine \(0\)’dır. İkinci satırdaki \(p(1{,}5) = 2{,}4375\) Vandermonde yoluyla bulduğumuz değerdir; \(p(2{,}5) = \frac{41}{16} = 2{,}5625\)’tir. Taban polinomlarını ve onların ağırlıklı toplamını yan yana çizelim:
Lagrange biçimi kavramsal olarak açıktır, ama her yeni noktada bütün \(L_k\)’leri baştan hesaplar ve nokta başına maliyeti \(n^2\) mertebesindedir. Aynı formülün küçük bir cebirle elde edilen bir yazılışı hem daha hızlı hem de sayısal olarak daha güvenilirdir.
Önerme 12.2 (Barisentrik İnterpolasyon Formülü) \(x_0, \dots, x_n\) birbirinden farklı düğümler ve
\[ w_k = \frac{1}{\prod_{j \ne k} (x_k - x_j)}, \qquad k = 0, 1, \dots, n \]
olsun. Düğümlerden farklı her \(x\) için interpolasyon polinomu
\[ p(x) = \frac{\displaystyle\sum_{k=0}^{n} \frac{w_k\, y_k}{x - x_k}}{\displaystyle\sum_{k=0}^{n} \frac{w_k}{x - x_k}} \]
eşitliğini sağlar.
İspat
\(\ell(x) = \prod_{j=0}^{n} (x - x_j)\) yazalım. \(x\) bir düğüm değilse \(L_k(x)\)’in payı \(\ell(x)/(x - x_k)\)’ye, paydası da \(1/w_k\)’ye eşittir; yani
\[ L_k(x) = \ell(x)\,\frac{w_k}{x - x_k} \]
olur ve Lagrange biçimi
\[ p(x) = \ell(x) \sum_{k=0}^{n} \frac{w_k\, y_k}{x - x_k} \]
hâlini alır. Aynı eşitliği bütün \(y_k\) değerleri \(1\) olan veriye uygulayalım. Bu verinin interpolasyon polinomu sabit \(1\) polinomudur, çünkü bu polinom bütün düğümlerden geçer ve teklik (Teorem 12.1) başka bir polinoma yer bırakmaz. Dolayısıyla
\[ 1 = \ell(x) \sum_{k=0}^{n} \frac{w_k}{x - x_k} \]
olur. Önceki eşitliği bu eşitliğe bölmek \(\ell(x)\) çarpanını yok eder ve iddia edilen formülü verir. \(\blacksquare\)
Yani \(w_k\) ağırlıkları düğümlerden bir kez, \(n^2\) mertebesinde işlemle hesaplanır; sonra her yeni noktada polinomun değeri \(n\) mertebesinde işlemle bulunur. Formülün bir güzelliği de bütün \(w_k\)’ler aynı sabitle çarpıldığında sonucun değişmemesidir, çünkü pay ve payda aynı sabitle çarpılır. SciPy’ın BarycentricInterpolator sınıfı bu formülü kullanır ve ağırlıkları, çok büyük ya da çok küçük sayılar oluşmasın diye ortak bir sabitle ölçekler:
import numpy as np
from scipy.interpolate import BarycentricInterpolator
xn = np.array([0.0, 1.0, 2.0, 3.0])
yn = np.array([1.0, 3.0, 2.0, 5.0])
diff = xn[:, None] - xn[None, :] # diff[k, j] = x_k - x_j
np.fill_diagonal(diff, 1.0) # j = k çarpanını dışarıda bırak
w = 1 / diff.prod(axis=1)
print(w)
b = BarycentricInterpolator(xn, yn)
print(b(np.array([1.5, 2.5])))
print(b.wi / w) # SciPy'ın ağırlıkları w'nin katıÇıktı:
[-0.16666667 0.5 -0.5 0.16666667]
[2.4375 2.5625]
[0.421875 0.421875 0.421875 0.421875]Ağırlıklarımız \(-\frac16\), \(\frac12\), \(-\frac12\), \(\frac16\)’dır; SciPy’ın ağırlıkları bunların \(0{,}421875 = \frac{27}{64}\) katıdır ve değerler yine \(2{,}4375\) ile \(2{,}5625\) çıktı. SciPy’daki eski scipy.interpolate.lagrange fonksiyonu sonucu tek terimlilerin katsayılarıyla verir, çok sayıda düğümde güvenilmezdir ve kullanımdan kaldırılmıştır; onun yerine BarycentricInterpolator kullanılmalıdır.
İnterpolasyon polinomu bir \(f\) fonksiyonunun değerlerinden kurulduğunda, polinomun \(f\)’den ne kadar saptığını şu teorem söyler.
Teorem 12.2 (İnterpolasyon Hatası) \(x_0, \dots, x_n\) noktaları \([a, b]\) aralığında birbirinden farklı, \(f \in C^{n+1}[a, b]\) ve \(p\), \(f\)’nin bu düğümlerdeki interpolasyon polinomu olsun. Her \(x \in [a, b]\) için
\[ f(x) - p(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!}\,\omega(x), \qquad \omega(x) = \prod_{k=0}^{n} (x - x_k) \]
eşitliğini sağlayan ve \(x\)’e bağlı olan bir \(\xi \in (a, b)\) sayısı vardır.
İspatı için bkz. Nümerik Analiz. Hata iki çarpandan oluşur: \(f^{(n+1)}(\xi)/(n+1)!\) çarpanı fonksiyona, \(\omega(x)\) çarpanı ise yalnız düğümlerin yerine bağlıdır. Fonksiyonu değiştiremeyiz, ama düğümleri seçebiliyorsak ikinci çarpanı küçük tutmaya çalışmalıyız. Bunun neden önemli olduğunu Runge olayı gösterir.
12.5 Runge Olayı
Düğüm sayısını artırınca interpolasyon polinomunun fonksiyona yaklaşmasını bekleriz. Weierstrass Yaklaşım Teoremi de kapalı bir aralıkta sürekli her fonksiyona polinomlarla istenildiği kadar yaklaşılabildiğini söyler. Ama bu, eşit aralıklı noktalardaki interpolasyonun fonksiyona yaklaştığı anlamına gelmez. Carl Runge’un 1901’de verdiği örneği inceleyelim:
\[ f(x) = \frac{1}{1 + 25x^2}, \qquad x \in [-1, 1]. \]
\(f\)’yi \(n + 1\) eşit aralıklı düğümde interpolasyonla yaklaşıyor ve en büyük hatayı \(2001\) noktalı sık bir ızgarada ölçüyoruz:
import numpy as np
from scipy.interpolate import BarycentricInterpolator
def f(x):
return 1 / (1 + 25 * x**2)
xx = np.linspace(-1, 1, 2001) # hatayı ölçtüğümüz sık ızgara
for n in range(2, 21, 2):
xn = np.linspace(-1, 1, n + 1) # n + 1 eşit aralıklı düğüm
p = BarycentricInterpolator(xn, f(xn))
err = np.max(np.abs(p(xx) - f(xx)))
print(f"n = {n:2d} en büyük hata = {err:9.3e}")Çıktı:
n = 2 en büyük hata = 6.462e-01
n = 4 en büyük hata = 4.384e-01
n = 6 en büyük hata = 6.169e-01
n = 8 en büyük hata = 1.045e+00
n = 10 en büyük hata = 1.916e+00
n = 12 en büyük hata = 3.663e+00
n = 14 en büyük hata = 7.195e+00
n = 16 en büyük hata = 1.439e+01
n = 18 en büyük hata = 2.919e+01
n = 20 en büyük hata = 5.982e+01Hata \(n = 4\)’ten sonra azalmıyor, büyüyor ve \(n = 20\)’de \(60\)’a yaklaşıyor. Oysa \(f\) son derece düzgün bir fonksiyondur, her mertebeden türevi vardır.
Tanım 12.4 (Runge Olayı) Eşit aralıklı düğümlerde kurulan interpolasyon polinomlarının, derece arttıkça aralığın uçlarına yakın giderek büyüyen salınımlar yapmasına ve en büyük hatanın sıfıra gitmek yerine büyümesine Runge olayı (Runge’s phenomenon) denir.
Yani daha çok veri her zaman daha iyi bir polinom demek değildir. Hata teoremindeki iki çarpana bakalım. Runge fonksiyonunun yüksek mertebeden türevleri \(n\) ile çok hızlı büyür. \(\omega(x)\) ise eşit aralıklı düğümlerde aralığın ortasında küçük, uçlarına yakın büyüktür; salınımların uçlarda başlaması bundandır. Çare, düğümleri \(\omega\)’yı aralık boyunca küçük tutacak biçimde seçmektir.
12.6 Çebişev Düğümleri
Düğümleri aralığın uçlarına doğru sıklaştırırsak \(\omega\)’yı aralığın her yerinde aynı ölçüde küçük tutabiliriz.
Tanım 12.5 (Çebişev Düğümleri) \(n \ge 0\) için \([-1, 1]\) aralığındaki
\[ x_k = \cos\frac{(2k+1)\pi}{2n+2}, \qquad k = 0, 1, \dots, n \]
sayılarına \(n + 1\) Çebişev düğümü (Chebyshev nodes) denir. Bir \([a, b]\) aralığındaki Çebişev düğümleri, bu sayıların \(t \mapsto \frac{a+b}{2} + \frac{b-a}{2}\,t\) dönüşümüyle taşınmasıyla elde edilir.
Yani birim çemberin üst yarısını \(n + 1\) eşit yaya bölüp her yayın orta noktasını \(x\) eksenine dik olarak izdüşürüyoruz. Çember üzerinde eşit aralıklı olan noktalar, eksen üzerinde uçlara doğru sıklaşır. Bu seçimin değeri şu önermededir.
Önerme 12.3 (Çebişev Düğümlerinde Düğüm Polinomu) \(x_0, \dots, x_n\), \([-1, 1]\) aralığındaki \(n + 1\) Çebişev düğümü ve \(\omega(x) = (x - x_0) \cdots (x - x_n)\) olsun. Her \(x \in [-1, 1]\) için \(|\omega(x)| \le 2^{-n}\)’dir ve bu sınıra \(x = \pm 1\)’de ulaşılır.
İspat
\(T_0(x) = 1\), \(T_1(x) = x\) ve \(k \ge 1\) için \(T_{k+1}(x) = 2x\,T_k(x) - T_{k-1}(x)\) yinelemesiyle tanımlanan polinomlara bakalım. Tümevarımla, \(k \ge 1\) için \(T_k\)’nin derecesi \(k\) ve baş katsayısı \(2^{k-1}\)’dir: \(T_1\) için doğrudur ve yinelemede en yüksek dereceli terim yalnız \(2x\,T_k\)’den gelir.
Her \(\theta\) için \(T_k(\cos\theta) = \cos k\theta\) olduğunu yine tümevarımla gösterelim. \(k = 0\) ve \(k = 1\) için açıktır. Kosinüsün toplam formülünden gelen
\[ \cos(k+1)\theta + \cos(k-1)\theta = 2\cos\theta\cos k\theta \]
özdeşliğinden ve tümevarım hipotezinden
\[ T_{k+1}(\cos\theta) = 2\cos\theta\,\cos k\theta - \cos(k-1)\theta = \cos(k+1)\theta \]
bulunur.
Şimdi \(T_{n+1}\)’e bakalım. \(\theta_k = \frac{(2k+1)\pi}{2n+2}\) için \((n+1)\theta_k = \frac{(2k+1)\pi}{2}\) olduğundan \(T_{n+1}(x_k) = \cos(n+1)\theta_k = 0\)’dır. \(\theta_0 < \theta_1 < \dots < \theta_n\) sayıları \((0, \pi)\) aralığında olduğundan \(x_k = \cos\theta_k\) düğümleri birbirinden farklıdır. Böylece derecesi \(n + 1\) olan \(T_{n+1}\)’in \(n + 1\) kökü tam olarak düğümlerdir. Baş katsayısı \(2^n\) olduğundan \(T_{n+1}(x) = 2^n\,\omega(x)\) olur. \(x \in [-1, 1]\) ise bir \(\theta \in [0, \pi]\) için \(x = \cos\theta\) yazılabilir ve
\[ |\omega(x)| = 2^{-n}\,|\cos(n+1)\theta| \le 2^{-n} \]
bulunur. \(x = 1\) (\(\theta = 0\)) ve \(x = -1\) (\(\theta = \pi\)) noktalarında eşitlik vardır. \(\blacksquare\)
İspatta kullanılan \(T_k\) polinomlarına Çebişev polinomları denir; NumPy’da numpy.polynomial.Chebyshev sınıfıyla temsil edilirler. Önermeyi eşit aralıklı düğümlerle karşılaştıralım. \(\omega\)’nın en büyük mutlak değerini yine sık bir ızgarada ölçüyoruz:
import numpy as np
def omega(xn, x):
"""omega(x) = (x - x_0)(x - x_1)...(x - x_n)"""
return np.prod(x[:, None] - xn, axis=1)
def cheb_nodes(n):
"""[-1, 1] aralığında n + 1 Çebişev düğümü."""
k = np.arange(n + 1)
return np.cos((2 * k + 1) * np.pi / (2 * n + 2))
xx = np.linspace(-1, 1, 2001)
for n in (5, 10, 15, 20):
eq = np.max(np.abs(omega(np.linspace(-1, 1, n + 1), xx)))
ch = np.max(np.abs(omega(cheb_nodes(n), xx)))
print(f"n = {n:2d} eşit: {eq:9.3e} Çebişev: {ch:9.3e}"
f" 2^(-n): {2.0**-n:9.3e}")Çıktı:
n = 5 eşit: 6.923e-02 Çebişev: 3.125e-02 2^(-n): 3.125e-02
n = 10 eşit: 8.532e-03 Çebişev: 9.766e-04 2^(-n): 9.766e-04
n = 15 eşit: 1.345e-03 Çebişev: 3.052e-05 2^(-n): 3.052e-05
n = 20 eşit: 2.336e-04 Çebişev: 9.537e-07 2^(-n): 9.537e-07Çebişev düğümlerindeki en büyük değer önermedeki \(2^{-n}\) sınırına eşit; eşit aralıklı düğümlerdeki değer ise \(n = 20\)’de bunun yaklaşık \(245\) katı. \(n = 10\) için iki \(\omega\) grafiği farkı açıkça gösteriyor:
Şimdi Runge fonksiyonunu iki tür düğümle, \(n = 10\) için çizelim. Kod iki alt grafikli bir şekil kurar (Matplotlib ile Grafik Çizimi); sharey=True iki grafiğe aynı \(y\) eksenini verir.
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import BarycentricInterpolator
def f(x):
return 1 / (1 + 25 * x**2)
n = 10
k = np.arange(n + 1)
nodes = {"eşit aralıklı düğümler": np.linspace(-1, 1, n + 1),
"Çebişev düğümleri": np.cos((2 * k + 1) * np.pi / (2 * n + 2))}
xx = np.linspace(-1, 1, 401)
fig, axes = plt.subplots(1, 2, figsize=(9, 3.5), sharey=True)
for ax, (title, xn) in zip(axes, nodes.items()):
p = BarycentricInterpolator(xn, f(xn))
ax.plot(xx, f(xx), color="black", label="f")
ax.plot(xx, p(xx), label="P10")
ax.plot(xn, f(xn), "o")
ax.set_title(title)
ax.set_ylim(-0.6, 2.1)
ax.legend()
err = np.max(np.abs(p(xx) - f(xx)))
print(f"{title}: en büyük hata {err:.4f}")
fig.savefig("runge.png", dpi=150)Çıktı:
eşit aralıklı düğümler: en büyük hata 1.9156
Çebişev düğümleri: en büyük hata 0.1092Kod en büyük hataları yazdırır ve grafiği runge.png dosyasına kaydeder. Aşağıdaki şekil aynı grafiğin aynı verilerle çizilmiş hâlidir:
Derece arttıkça iki düğüm türünün nasıl davrandığını bir tabloda görelim. NumPy’ın Chebyshev.interpolate sınıf metodu, bir fonksiyonu Çebişev düğümlerinde interpolasyonla yaklaşan polinomu tek satırda kurar:
import numpy as np
from numpy.polynomial import Chebyshev
from scipy.interpolate import BarycentricInterpolator
def f(x):
return 1 / (1 + 25 * x**2)
xx = np.linspace(-1, 1, 2001)
for n in range(4, 41, 4):
k = np.arange(n + 1)
xe = np.linspace(-1, 1, n + 1)
xc = np.cos((2 * k + 1) * np.pi / (2 * n + 2))
ee = np.max(np.abs(BarycentricInterpolator(xe, f(xe))(xx) - f(xx)))
ec = np.max(np.abs(BarycentricInterpolator(xc, f(xc))(xx) - f(xx)))
print(f"n = {n:2d} eşit: {ee:9.2e} Çebişev: {ec:9.2e}")
c = Chebyshev.interpolate(f, 10) # Çebişev düğümlerinde interpolasyon
print(f"Chebyshev.interpolate, n = 10: {np.max(np.abs(c(xx) - f(xx))):.2e}")Çıktı:
n = 4 eşit: 4.38e-01 Çebişev: 4.02e-01
n = 8 eşit: 1.05e+00 Çebişev: 1.71e-01
n = 12 eşit: 3.66e+00 Çebişev: 6.92e-02
n = 16 eşit: 1.44e+01 Çebişev: 3.26e-02
n = 20 eşit: 5.98e+01 Çebişev: 1.53e-02
n = 24 eşit: 2.57e+02 Çebişev: 6.95e-03
n = 28 eşit: 1.13e+03 Çebişev: 3.08e-03
n = 32 eşit: 5.06e+03 Çebişev: 1.40e-03
n = 36 eşit: 2.29e+04 Çebişev: 6.41e-04
n = 40 eşit: 1.05e+05 Çebişev: 2.89e-04
Chebyshev.interpolate, n = 10: 1.09e-01Eşit aralıklı düğümlerde hata \(n = 40\)’ta \(10^5\)’i geçti. Çebişev düğümlerinde ise düzenli olarak azalıyor: \(n\) her \(4\) arttığında kabaca \(2{,}2\) kat küçülüyor. Chebyshev.interpolate aynı \(0{,}109\) hatasını verdi, çünkü o da aynı düğümleri kullanır. Türevi sınırlı her \(f\) için Çebişev düğümlerindeki interpolasyon polinomlarının \(f\)’ye düzgün yakınsadığı bilinir; eşit aralıklı düğümlerde böyle bir güvence yoktur. Düğümleri seçebildiğimiz her durumda Çebişev düğümlerini tercih etmeliyiz.
12.7 Parçalı Polinomlar ve Spline
Düğümleri seçemiyorsak, örneğin veri eşit aralıklı bir tablodan geliyorsa, başka bir yol izleriz: dereceyi artırmak yerine aralığı küçük parçalara bölüp her parçada düşük dereceli bir polinom kullanırız. En basit seçim, ardışık düğümleri doğru parçalarıyla birleştirmektir ve NumPy’da bunu np.interp yapar:
import numpy as np
def f(x):
return 1 / (1 + 25 * x**2)
xn = np.linspace(-1, 1, 11) # 10 eşit alt aralık
xx = np.linspace(-1, 1, 2001)
L = np.interp(xx, xn, f(xn)) # parçalı doğrusal
print(f"en büyük hata: {np.max(np.abs(L - f(xx))):.4f}")
print(np.interp([0.05, 0.5], xn, f(xn)))Çıktı:
en büyük hata: 0.0674
[0.875 0.15 ]np.interp(x, xn, yn), artan sırada verilmiş xn düğümleri arasında doğrusal interpolasyon yapar: \(0{,}05\) noktası \(0\) ile \(0{,}2\) düğümleri arasında kaldığından değer \(1\) ile \(0{,}5\)’in ağırlıklı ortalaması olan \(0{,}875\)’tir. Aynı \(11\) düğümle en büyük hata \(0{,}067\); aynı düğümlerdeki \(P_{10}\) polinomunun hatası ise \(1{,}92\)’ydi. Bu yöntemin hatasına hata teoremi doğrudan bir sınır verir.
Önerme 12.4 (Parçalı Doğrusal İnterpolasyonun Hatası) \(f \in C^2[a, b]\), her \(t \in [a, b]\) için \(|f''(t)| \le M\) ve \(a = x_0 < x_1 < \dots < x_n = b\) ardışık düğümleri arasındaki en büyük uzaklık \(h\) olsun. Düğümlerde \(f\) ile çakışan parçalı doğrusal \(L\) fonksiyonu için her \(x \in [a, b]\) noktasında
\[ |f(x) - L(x)| \le \frac{M h^2}{8} \]
olur.
İspat
\(x\), bir \([x_i, x_{i+1}]\) aralığında olsun. Bu aralıkta \(L\), \(f\)’nin \(x_i\) ve \(x_{i+1}\) düğümlerindeki birinci dereceden interpolasyon polinomudur. Hata teoremi (Teorem 12.2) \(n = 1\) ile bir \(\xi\) için
\[ |f(x) - L(x)| = \frac{|f''(\xi)|}{2}\,(x - x_i)(x_{i+1} - x) \]
verir. Toplamı sabit olan iki pozitif sayının çarpımı, sayılar eşitken en büyüktür; bu yüzden \((x - x_i)(x_{i+1} - x)\) çarpımı aralığın orta noktasında en büyük değeri olan \(\frac{(x_{i+1} - x_i)^2}{4}\)’ü alır. \(|f''(\xi)| \le M\) ve \(x_{i+1} - x_i \le h\) olduğundan
\[ |f(x) - L(x)| \le \frac{M}{2} \cdot \frac{h^2}{4} = \frac{M h^2}{8} \]
bulunur. \(\blacksquare\)
Yani düğüm sayısını iki katına çıkarmak hatayı yaklaşık dörtte birine indirir. Daha iyisini yapmak için parçaları daha yüksek dereceden seçer ve birleşme noktalarında türevlerin de uyuşmasını isteriz.
Tanım 12.6 (Kübik Spline) \(a = x_0 < x_1 < \dots < x_n = b\) düğümleri ve \(y_0, y_1, \dots, y_n\) değerleri verilsin. Aşağıdaki üç koşulu sağlayan \(S : [a, b] \to \mathbb{R}\) fonksiyonuna bu verinin bir kübik spline’ı (cubic spline) denir:
- Her \([x_i, x_{i+1}]\) alt aralığında \(S\), derecesi en fazla \(3\) olan bir polinomdur.
- Her \(i = 0, 1, \dots, n\) için \(S(x_i) = y_i\)’dir.
- \(S\), \(S'\) ve \(S''\) fonksiyonları \([a, b]\) aralığında süreklidir.
Yani spline, düğümlerde birbirine türevleriyle birlikte pürüzsüzce eklenmiş kübik parçalardan oluşur. Adını, teknik çizimde eğri çizmek için kullanılan esnek cetvelden alır. Bilinmeyenleri sayalım. \(n\) parçanın her biri \(4\) katsayı taşır; toplam \(4n\) bilinmeyen vardır. Her parçanın iki ucunda veriyle çakışması \(2n\) koşul, iç düğümlerde \(S'\) ile \(S''\)’nün sürekliliği \(2(n-1)\) koşul verir. Toplam \(4n - 2\) koşul eder; spline’ı tek kılmak için iki koşul daha gerekir. Bu iki koşul genellikle uç noktalara konur ve her seçim ayrı bir spline türü verir.
Tanım 12.7 (Doğal Spline) \(S''(a) = S''(b) = 0\) koşullarını sağlayan kübik spline’a doğal spline (natural spline) denir.
Yani doğal spline’ın uçlarda eğriliği sıfırdır: uç düğümlerin ötesinde bir doğru gibi devam etmeye hazırdır. Fonksiyon hakkında ek bilgi istemez, ama gerçek \(f''\) uçlarda sıfır değilse orada hatası büyür (Alıştırma 12.8).
Tanım 12.8 (Sabitlenmiş Spline) Değerleri spline’a verilen fonksiyon \(f\) olsun. \(S'(a) = f'(a)\) ve \(S'(b) = f'(b)\) koşullarını sağlayan kübik spline’a sabitlenmiş spline (clamped spline) denir.
Yani uçlardaki eğimi biz veririz. Ek bilgi kullandığı için genellikle en doğru seçenektir; ama uçlardaki türevler çoğu zaman bilinmez.
Tanım 12.9 (Düğüm Olmama Koşulu) İlk iki parçanın ve son iki parçanın aynı kübik polinom olmasını, yani \(S'''\)’nün \(x_1\) ve \(x_{n-1}\) düğümlerinde de sürekli olmasını isteyen koşula düğüm olmama koşulu (not-a-knot condition) denir.
Yani \(x_1\) ve \(x_{n-1}\) noktalarında parça değişmez; bu iki nokta yalnız spline’ın geçtiği noktalar olarak kalır. Değeri, türevleri ve üçüncü türevi bir noktada çakışan iki kübik polinom özdeş olduğundan iki tanım aynı şeyi söyler. Bu koşul ek bilgi istemez ve uçlarda doğal spline’dan daha doğrudur. SciPy’ın varsayılan seçimi budur.
- Düğümleri artan sırada bir
xdizisine, değerleri birydizisine koyun. S = CubicSpline(x, y, bc_type=...)ile spline’ı kurun. Seçenekler"not-a-knot"(varsayılan),"natural"ve uç türevlerini veren((1, da), (1, db))biçimidir.S(t)ile değer,S(t, 1)veS(t, 2)ile türev,S.integrate(a, b)ile integral hesaplayın.
Örnek 12.1 (Runge Fonksiyonunun Spline’ı) \(f(x) = \frac{1}{1 + 25x^2}\) fonksiyonunun \([-1, 1]\) aralığındaki \(11\) eşit aralıklı düğümdeki değerlerinden kübik spline kurunuz. En büyük hatayı bulunuz; \(x = 0{,}3\)’teki değeri ve türevi ile \([-1, 1]\) üzerindeki integrali gerçek değerlerle karşılaştırınız.
Çözüm
Tarifin üç adımını uygulayalım. Gerçek türev \(f'(x) = -\frac{50x}{(1 + 25x^2)^2}\), \(x = 0{,}3\)’te \(-\frac{15}{(1 + 2{,}25)^2}\)’dir. Gerçek integral, \(\int \frac{dx}{1 + 25x^2} = \frac{1}{5}\arctan 5x\) olduğundan \(\frac{2}{5}\arctan 5\)’tir.
import numpy as np
from scipy.interpolate import CubicSpline
def f(x):
return 1 / (1 + 25 * x**2)
xn = np.linspace(-1, 1, 11) # 10 eşit alt aralık
xx = np.linspace(-1, 1, 2001)
S = CubicSpline(xn, f(xn)) # varsayılan: not-a-knot
print(f"en büyük hata: {np.max(np.abs(S(xx) - f(xx))):.4f}")
print(S.c.shape) # 4 katsayı, 10 parça
print(S(0.3), f(0.3))
print(S(0.3, 1), -15 / (1 + 25 * 0.3**2)**2) # birinci türev
print(S.integrate(-1, 1), 0.4 * np.arctan(5)) # integralÇıktı:
en büyük hata: 0.0220
(4, 10)
0.29733288239958944 0.3076923076923077
-1.3660027056024633 -1.4201183431952662
0.5519677815614747 0.5493603067780064Spline’ın en büyük hatası \(0{,}022\)’dir; bu, aynı düğümlerdeki \(P_{10}\) polinomunun hatasının (\(1{,}92\)) yaklaşık \(87\)’de biridir. S.c dizisinin şekli \((4, 10)\)’dur: \(i\)’inci sütun, \([x_i, x_{i+1}]\) parçasında
\[ S(x) = c_0 (x - x_i)^3 + c_1 (x - x_i)^2 + c_2 (x - x_i) + c_3 \]
yazılışının dört katsayısıdır. \(x = 0{,}3\)’te değer \(0{,}2973\) (gerçeği \(0{,}3077\)), türev \(-1{,}366\) (gerçeği \(-1{,}420\)), integral \(0{,}55197\) (gerçeği \(0{,}54936\)) çıktı. \(\blacksquare\)
Spline’ın parçalarını farklı renklerle çizersek birleşme yerleri görünür:
CubicSpline ile kurulan kübik spline. Renkler on kübik parçayı birbirinden ayırır; parçalar düğümlerde birinci ve ikinci türevleriyle birlikte birleşir. En büyük hata 0,022'dir; aynı düğümlerdeki P10 polinomunun hatası 1,92'ydi.Düğüm sayısı artınca spline’ın hatası, parçalı doğrusal interpolasyonunkinden çok daha hızlı azalır. Bu hızı ölçelim:
import numpy as np
from scipy.interpolate import CubicSpline
def f(x):
return 1 / (1 + 25 * x**2)
xx = np.linspace(-1, 1, 20001)
for n in (10, 20, 40, 80, 160, 320): # alt aralık sayısı
xn = np.linspace(-1, 1, n + 1)
e_lin = np.max(np.abs(np.interp(xx, xn, f(xn)) - f(xx)))
e_spl = np.max(np.abs(CubicSpline(xn, f(xn))(xx) - f(xx)))
print(f"n = {n:3d} h = {2 / n:.5f} doğrusal: {e_lin:8.2e}"
f" spline: {e_spl:8.2e}")Çıktı:
n = 10 h = 0.20000 doğrusal: 6.74e-02 spline: 2.20e-02
n = 20 h = 0.10000 doğrusal: 4.18e-02 spline: 3.18e-03
n = 40 h = 0.05000 doğrusal: 1.40e-02 spline: 2.78e-04
n = 80 h = 0.02500 doğrusal: 3.80e-03 spline: 1.61e-05
n = 160 h = 0.01250 doğrusal: 9.70e-04 spline: 9.68e-07
n = 320 h = 0.00625 doğrusal: 2.44e-04 spline: 5.98e-08\(n = 40\)’tan sonra alt aralık sayısı iki katına çıkınca parçalı doğrusal interpolasyonun hatası yaklaşık \(4\) kat, spline’ınki ise yaklaşık \(16\) kat küçülüyor; örneğin \(9{,}70 \cdot 10^{-4} / 2{,}44 \cdot 10^{-4} \approx 4{,}0\) ve \(9{,}68 \cdot 10^{-7} / 5{,}98 \cdot 10^{-8} \approx 16{,}2\). Daha küçük \(n\)’de oranlar daha küçüktür (doğrusal için \(1{,}6\) ve \(3{,}0\), spline için \(6{,}9\) ve \(11{,}4\)), çünkü \(h\), fonksiyonun \(0\) çevresindeki dar tepesine göre henüz büyüktür. İki eksen de logaritmik olunca, \(h^p\) ile orantılı bir hata eğimi \(-p\) olan bir doğru çizer:
np.interp) ve kübik spline'ın (CubicSpline) en büyük hatası, iki ekseni de logaritmik ölçekte; tablodaki sayıların aynısı. n = 40'tan sonra alt aralık sayısı iki katına çıkınca hatalar yaklaşık 4 ve 16 kat küçülür; daha küçük n'de iki eğri de daha yavaş iner. Kesikli çizgiler bu hızlara karşılık gelen −2 ve −4 eğimlerini gösterir.Spline’daki \(h^4\) hızı genel bir teoremdir.
Teorem 12.3 (Sabitlenmiş Spline’ın Hatası) \(f \in C^4[a, b]\), her \(t \in [a, b]\) için \(|f^{(4)}(t)| \le M\) ve \(S\), \(f\)’nin \(a = x_0 < \dots < x_n = b\) düğümlerindeki sabitlenmiş spline’ı olsun. Ardışık düğümler arasındaki en büyük uzaklık \(h\) ise her \(x \in [a, b]\) için
\[ |f(x) - S(x)| \le \frac{5}{384}\,M h^4 \]
olur.
İspatı bu bölümün kapsamını aşıyor; sonucu sayısal olarak doğrulayacağız. Düğüm olmama koşulu için de hata \(h^4\) ile orantılı olarak azalır. Doğal spline’da ise \(f''\) uçlarda sıfır değilse hata uçlara yakın yalnız \(h^2\) hızıyla azalır.
Örnek 12.2 (Sinüsün Sabitlenmiş Spline’ı) \([0, \pi]\) aralığında \(\sin x\) fonksiyonunu \(n = 4, 8, 16, 32, 64\) eşit alt aralıkla sabitlenmiş spline ile yaklaşınız. En büyük hatayı teoremdeki sınırla (Teorem 12.3) karşılaştırınız.
Çözüm
Uçlardaki türevler \(\cos 0 = 1\) ve \(\cos\pi = -1\)’dir. bc_type=((1, 1.0), (1, -1.0)) seçeneği her uçta birinci türevi (1) verilen değere eşitler. \(|\sin^{(4)} x| = |\sin x| \le 1\) olduğundan \(M = 1\) alırız; \(h = \pi/n\)’dir.
import numpy as np
from scipy.interpolate import CubicSpline
xx = np.linspace(0, np.pi, 20001)
prev = None
for n in (4, 8, 16, 32, 64):
xn = np.linspace(0, np.pi, n + 1)
h = np.pi / n
# sabitlenmiş spline: S'(0) = cos 0 = 1, S'(pi) = cos pi = -1
S = CubicSpline(xn, np.sin(xn), bc_type=((1, 1.0), (1, -1.0)))
err = np.max(np.abs(S(xx) - np.sin(xx)))
bound = 5 / 384 * h**4 # max |sin''''| = 1
ratio = "" if prev is None else f" oran: {prev / err:5.2f}"
print(f"n = {n:2d} hata: {err:8.2e} sınır: {bound:8.2e}{ratio}")
prev = errÇıktı:
n = 4 hata: 1.12e-03 sınır: 4.95e-03
n = 8 hata: 6.32e-05 sınır: 3.10e-04 oran: 17.77
n = 16 hata: 3.89e-06 sınır: 1.94e-05 oran: 16.26
n = 32 hata: 2.42e-07 sınır: 1.21e-06 oran: 16.06
n = 64 hata: 1.51e-08 sınır: 7.56e-08 oran: 16.01Hata her satırda sınırın altında kalıyor ve ardışık iki hatanın oranı \(16\)’ya yaklaşıyor: \(h\) yarıya inince hata \(2^4 = 16\) kat küçülüyor. Sınır gerçek hatanın yaklaşık \(5\) katıdır; teorem en kötü durumu kapsar. \(\blacksquare\)
12.8 En Küçük Kareler ile Polinom Uydurma
Ölçüm verisinde her değer bir miktar gürültü taşır ve verinin bütün noktalarından geçmeye çalışmak gürültüyü de modele katmak olur. Bunun yerine düşük dereceli bir \(p\) polinomu seçip \(m\) ölçüm \((x_i, y_i)\) için
\[ \sum_{i=1}^{m} \big(y_i - p(x_i)\big)^2 \]
toplamını en küçük yapan katsayıları ararız. Bu, sütunları \(1, x, \dots, x^d\) olan bir Vandermonde matrisiyle kurulan ve denklem sayısı bilinmeyen sayısından fazla olan bir sistemin en küçük kareler çözümüdür (bkz. NumPy ile Lineer Cebir ve Nümerik Analiz). NumPy’da bu işi iki fonksiyon yapar: Polynomial.fit(x, y, deg) ve eski arayüzden np.polyfit(x, y, deg). Uydurmadan sonra \(r_i = y_i - p(x_i)\) artıklarına bakmak, modelin yeterli olup olmadığını anlamanın en iyi yoludur.
Tanım 12.10 (Artık Grafiği) \((x_i, y_i)\) verisine bir \(\hat y\) modeli uydurulmuş olsun. \(r_i = y_i - \hat y(x_i)\) artıklarının \(x_i\)’ye karşı çizildiği grafiğe artık grafiği (residual plot) denir.
Yani artık grafiği, modelin açıklayamadığı kısmı gösterir. Model verinin yapısını yakalamışsa geriye yalnız gürültü kalır ve artıklar sıfırın iki yanına düzensizce dağılır. Artıklarda bir desen görülüyorsa (bir kemer, bir dalga ya da uçlara doğru açılan bir yelpaze) model bir şeyi kaçırıyordur.
- Bir derece seçin; mümkünse dereceyi problemin kendisinden ya da verinin grafiğinden belirleyin.
p = Polynomial.fit(x, y, deg)ile uydurun; alışılmış katsayılar içinp.convert().coefkullanın.r = y - p(x)artıklarını hesaplayıp çizin. Artıklarda bir desen varsa dereceyi artırın ya da modeli değiştirin.
Örnek 12.3 (Yukarı Atılan Bir Cismin Yüksekliği) \(25\) m yükseklikten \(3\) m/s hızla yukarı atılan bir cismin yüksekliği \(t \in [0, 2]\) s aralığında \(0{,}1\) s arayla \(21\) kez ölçülüyor. Ölçümleri üretmek için \(h(t) = 25 + 3t - 4{,}905t^2\) değerlerine standart sapması \(0{,}05\) m olan normal dağılımlı gürültü ekleniyor. Veriye birinci ve ikinci dereceden polinom uydurup artık kareleri toplamlarını karşılaştırınız.
Çözüm
rng.normal(0, 0.05, t.size) ortalaması \(0\), standart sapması \(0{,}05\) olan normal dağılımdan t.size tane sayı üretir (Rastgele Sayılar ve Monte Carlo Yöntemleri); tohum sabit olduğundan her çalıştırmada aynı veri çıkar. Tarifin ilk iki adımını iki derece için uygulayalım:
import numpy as np
from numpy.polynomial import Polynomial
rng = np.random.default_rng(16)
t = np.linspace(0, 2, 21) # s
h = 25 + 3 * t - 4.905 * t**2 + rng.normal(0, 0.05, t.size) # m
for deg in (1, 2):
p = Polynomial.fit(t, h, deg).convert()
r = h - p(t) # artıklar
print(f"derece {deg}: katsayılar {np.round(p.coef, 4)}")
print(f" artık kareleri toplamı {np.sum(r**2):.4f}")Çıktı:
derece 1: katsayılar [28.1307 -6.8326]
artık kareleri toplamı 54.0890
derece 2: katsayılar [25.022 2.9842 -4.9084]
artık kareleri toplamı 0.0435Doğrunun artık kareleri toplamı \(54{,}1\), parabolünkü \(0{,}0435\)’tir; parabol yaklaşık \(1200\) kat daha iyi uyuyor. Parabolün katsayıları \(25{,}022\); \(2{,}984\); \(-4{,}908\) ve verinin üretildiği \(25\); \(3\); \(-4{,}905\) değerlerine çok yakın. \(\blacksquare\)
Tarifin üçüncü adımı artıklara bakmaktır. Aşağıdaki kod üç alt grafik kurar: ilkinde veri ve iki polinom, ötekilerde iki modelin artıkları vardır.
import matplotlib.pyplot as plt
import numpy as np
from numpy.polynomial import Polynomial
rng = np.random.default_rng(16)
t = np.linspace(0, 2, 21)
h = 25 + 3 * t - 4.905 * t**2 + rng.normal(0, 0.05, t.size)
tt = np.linspace(0, 2, 201)
fig, axes = plt.subplots(1, 3, figsize=(12, 3.5))
axes[0].plot(t, h, "o", color="black", label="ölçüm")
for deg in (1, 2):
p = Polynomial.fit(t, h, deg)
r = h - p(t) # artıklar
axes[0].plot(tt, p(tt), label=f"derece {deg}")
axes[deg].plot(t, r, "o")
axes[deg].axhline(0, color="gray", lw=0.8)
axes[deg].set_title(f"derece {deg}: artıklar")
print(f"derece {deg}: en büyük |artık| {np.max(np.abs(r)):.3f} m")
axes[0].set_xlabel("t (s)")
axes[0].set_ylabel("h (m)")
axes[0].legend()
fig.savefig("artiklar.png", dpi=150)Çıktı:
derece 1: en büyük |artık| 3.161 m
derece 2: en büyük |artık| 0.078 mKod en büyük artıkları yazdırır ve grafiği artiklar.png dosyasına kaydeder. Aşağıdaki şekil aynı grafiğin aynı verilerle çizilmiş hâlidir:
Doğrunun artıkları düzgün bir kemer çiziyor: ortada pozitif, uçlarda negatif. Bu, modelin verideki eğriliği yakalayamadığını söyler. Parabolün artıkları ise sıfırın iki yanına düzensizce dağılıyor ve büyüklükleri, gürültünün \(0{,}05\) m olan standart sapmasıyla uyumlu.
Modelin katsayılarının fiziksel bir anlamı da var. \(h(t) = h_0 + v_0 t - \frac{g}{2}t^2\) olduğundan \(t^2\)’nin katsayısının \(-2\) katı yerçekimi ivmesinin bir tahminidir. np.polyfit, cov=True seçeneğiyle katsayıların kovaryans matrisini de verir; bu matrisin köşegenindeki sayıların karekökleri katsayıların standart hatalarıdır. polyfit’in katsayıları en yüksek dereceden başlattığına dikkat edin:
import numpy as np
rng = np.random.default_rng(16)
t = np.linspace(0, 2, 21)
h = 25 + 3 * t - 4.905 * t**2 + rng.normal(0, 0.05, t.size)
c, C = np.polyfit(t, h, 2, cov=True) # c = [c2, c1, c0]
se = np.sqrt(np.diag(C)) # katsayıların standart hataları
g, se_g = -2 * c[0], 2 * se[0] # h = h0 + v0 t - (g/2) t^2
print(np.round(c, 4))
print(np.round(se, 4))
print(f"g ≈ {g:.3f} ± {se_g:.3f} m/s^2")Çıktı:
[-4.9084 2.9842 25.022 ]
[0.0328 0.068 0.0293]
g ≈ 9.817 ± 0.066 m/s^2Tahmin \(g \approx 9{,}817 \pm 0{,}066\) m/s² oldu; gerçek \(9{,}81\) değeri tahminin bir standart hata yakınında. NumPy bu standart hataları artıkların büyüklüğünden kestirir: gürültünün varyansını, artık kareleri toplamını serbestlik derecesine, yani \(21 - 3 = 18\)’e bölerek tahmin eder.
Derecesi en fazla \(d\) olan polinomlar, derecesi en fazla \(d + 1\) olanların arasındadır. Bu yüzden en küçük kareler uydurmasında artık kareleri toplamı dereceyle hiçbir zaman artmaz:
import numpy as np
from numpy.polynomial import Polynomial
rng = np.random.default_rng(16)
t = np.linspace(0, 2, 21)
h = 25 + 3 * t - 4.905 * t**2 + rng.normal(0, 0.05, t.size)
for deg in range(0, 9):
p = Polynomial.fit(t, h, deg)
rss = np.sum((h - p(t))**2)
print(f"derece {deg}: artık kareleri toplamı {rss:.5f}")Çıktı:
derece 0: artık kareleri toplamı 413.55893
derece 1: artık kareleri toplamı 54.08900
derece 2: artık kareleri toplamı 0.04348
derece 3: artık kareleri toplamı 0.04281
derece 4: artık kareleri toplamı 0.03167
derece 5: artık kareleri toplamı 0.02490
derece 6: artık kareleri toplamı 0.02473
derece 7: artık kareleri toplamı 0.02216
derece 8: artık kareleri toplamı 0.01933İkinci dereceden sonraki düşüşler çok küçüktür. Bunlar modelin iyileştiğini değil, polinomun gürültüye uymaya başladığını gösterir; buna aşırı uyum (overfitting) denir. Dereceyi artık kareleri toplamının küçüklüğüne bakarak değil, artık grafiğine ve modelin anlamına bakarak seçin. Uydurmada kullanılmayan test noktaları da iyi bir ölçüttür (Alıştırma 12.10).
12.9 Doğrusal Olmayan Modeller
Her model katsayılarına göre doğrusal değildir ve böyle modeller için polyfit yetmez.
Tanım 12.11 (Parametrelerine Göre Doğrusal Model) \(\varphi_1, \dots, \varphi_m\) bilinen fonksiyonlar ve \(c_1, \dots, c_m\) parametreler olmak üzere
\[ y = c_1\varphi_1(x) + c_2\varphi_2(x) + \dots + c_m\varphi_m(x) \]
biçimindeki modele parametrelerine göre doğrusal model denir. Bu biçimde yazılamayan modellere doğrusal olmayan model denir.
Yani “doğrusal” sözü \(x\)’e göre değil, parametrelere göredir. \(y = c_0 + c_1x + c_2x^2\) parabolü ya da \(y = c_1 \sin x + c_2 e^x\) doğrusal modellerdir ve bunların en küçük kareler problemi bir lineer sistemin en küçük kareler çözümüne indirgenir. \(y = a e^{-kx}\) ya da \(y = A\sin \omega t\) ise \(k\) ve \(\omega\)’ya göre doğrusal değildir. Doğrusal olmayan bir modelde artık kareleri toplamı parametrelerin genel bir fonksiyonudur; en küçük değeri, bir başlangıç tahmininden yola çıkan yinelemeli bir yöntemle aranır. Bu bir optimizasyon problemidir (Kök Bulma ve Optimizasyon) ve scipy.optimize.curve_fit fonksiyonu onu varsayılan olarak Levenberg–Marquardt yöntemiyle çözer.
- Modeli, ilk argümanı bağımsız değişken, ötekileri parametreler olan bir Python fonksiyonu olarak yazın:
def model(t, a, b, c). popt, pcov = curve_fit(model, t, y, p0=...)ile uydurun.p0, parametreler için verinin grafiğinden okunan makul bir başlangıç tahminidir.np.sqrt(np.diag(pcov))ile parametrelerin standart hatalarını hesaplayın ve artıklara bakın.
Örnek 12.4 (Newton Soğuma Yasası) Ortam sıcaklığı \(T_a\) olan bir odada soğuyan bir sıvının sıcaklığı
\[ T(t) = T_a + (T_0 - T_a)\,e^{-kt} \]
yasasına uyar. \(t = 0, 2, \dots, 30\) dakikalarında ölçülen sıcaklıklar, \(T_a = 22\), \(T_0 = 90\), \(k = 0{,}12\) değerlerine standart sapması \(0{,}5\) °C olan gürültü eklenerek üretiliyor. \(T_a\), \(T_0\) ve \(k\) parametrelerini curve_fit ile tahmin ediniz.
Çözüm
Model \(k\)’ye göre doğrusal değildir, bu yüzden tarifi uygularız. Başlangıç tahminini verinin kendisinden okuruz. Kod önce ilk ve son ölçümü yazdırıyor: son ölçüm \(23\) °C civarında olduğundan ortam sıcaklığı bunun biraz altında olmalı ve \(T_a\) için \(20\) deneriz. İlk ölçüm \(T_0\)’ın \(90\) civarında olduğunu gösteriyor; kaba bir tahmin olarak \(80\) yazmak yeter. \(k\) için \(0{,}1\) deneriz.
import numpy as np
from scipy.optimize import curve_fit
def model(t, Ta, T0, k):
"""Newton soğuma yasası: ortam Ta, başlangıç T0, hız sabiti k."""
return Ta + (T0 - Ta) * np.exp(-k * t)
rng = np.random.default_rng(5)
t = np.arange(0, 31, 2.0) # dakika
T = model(t, 22, 90, 0.12) + rng.normal(0, 0.5, t.size) # °C
print("ilk ve son ölçüm:", np.round(T[[0, -1]], 2))
popt, pcov = curve_fit(model, t, T, p0=(20, 80, 0.1))
perr = np.sqrt(np.diag(pcov))
for name, v, e in zip(("Ta", "T0", "k"), popt, perr):
print(f"{name} = {v:.4f} ± {e:.4f}")
r = T - model(t, *popt)
print(f"artık kareleri toplamı: {np.sum(r**2):.3f}")Çıktı:
ilk ve son ölçüm: [89.6 22.99]
Ta = 21.6914 ± 0.3530
T0 = 89.3752 ± 0.4294
k = 0.1173 ± 0.0022
artık kareleri toplamı: 3.439Tahminler \(T_a \approx 21{,}69 \pm 0{,}35\), \(T_0 \approx 89{,}38 \pm 0{,}43\) ve \(k \approx 0{,}1173 \pm 0{,}0022\) oldu. Verinin üretildiği \(22\), \(90\) ve \(0{,}12\) değerlerinin hepsi tahminlerden iki standart hatadan daha az uzakta. Artık kareleri toplamı \(3{,}44\)’tür. Bunu serbestlik derecesine, yani \(16 - 3 = 13\)’e bölünce \(0{,}26\) çıkar; bu, gürültünün \(0{,}25\) olan varyansına çok yakındır. \(\blacksquare\)
Uydurulan eğri ölçümlerin arasından geçiyor ve zamanla tahmin edilen ortam sıcaklığına yaklaşıyor:
curve_fit ile uydurulan Newton soğuma eğrisi. Kesikli çizgi, tahmin edilen ortam sıcaklığıdır; eğri zamanla ona yaklaşır.Yinelemeli yöntem, başlangıç tahmininin yakınındaki bir yerel minimuma gider. \(y = A\sin \omega t\) modelinde frekans için kötü bir başlangıç değeri yanlış bir sonuç verir. Aşağıdaki veri \(A = 2\), \(\omega = 3\) ile üretildi. Kodda \(\omega\), w adıyla yazılıyor; artık kareleri toplamı (residual sum of squares) da kısaca RSS diye yazdırılıyor:
import numpy as np
from scipy.optimize import curve_fit
def model(t, A, w):
return A * np.sin(w * t)
rng = np.random.default_rng(4)
t = np.linspace(0, 6, 40)
y = model(t, 2.0, 3.0) + rng.normal(0, 0.2, t.size)
for w0 in (1.0, 2.5, 4.0):
(A, w), _ = curve_fit(model, t, y, p0=(1.0, w0))
rss = np.sum((y - model(t, A, w))**2)
print(f"w0 = {w0}: A = {A:.4f} w = {w:.4f} RSS = {rss:.2f}")Çıktı:
w0 = 1.0: A = 0.2805 w = 1.7477 RSS = 78.95
w0 = 2.5: A = 1.9694 w = 2.9947 RSS = 1.75
w0 = 4.0: A = 0.2980 w = 4.2162 RSS = 78.73Başlangıç frekansı \(1\) ve \(4\) iken yöntem yanlış yerel minimumlarda durdu ve artık kareleri toplamı \(79\) civarında kaldı. Başlangıç frekansı \(2{,}5\) iken gerçek değerlere ulaştı ve toplam \(1{,}75\)’e indi. Başlangıç tahminini verinin grafiğinden okuyun (burada ardışık iki tepe arasındaki uzaklık \(2\pi/\omega\)’dır) ve sonucu her zaman artıklarla denetleyin.
12.10 Alıştırmalar
Aşağıdaki alıştırmaları önce kendiniz çözmeye çalışın; her çözümde kod ve çıktısı verilmiştir.
Alıştırma 12.1 (Bir Polinomun Yerel Ekstremumları) \(p(x) = 3x^4 - 4x^3 - 12x^2 + 5\) polinomunun yerel ekstremumlarını Polynomial sınıfının deriv ve roots metodlarıyla bulunuz ve ikinci türev testiyle sınıflandırınız.
Çözüm
Adım 1. Bir polinom her yerde türevlenebilir olduğundan yerel ekstremumlar ancak \(p'(x) = 0\) denkleminin köklerinde, yani kritik noktalarda olabilir. Bir kritik noktada \(p'' > 0\) ise yerel minimum, \(p'' < 0\) ise yerel maksimum vardır (bkz. Analiz 2).
Adım 2. deriv() birinci, deriv(2) ikinci türevi verir; kritik noktalar birinci türevin kökleridir:
from numpy.polynomial import Polynomial
p = Polynomial([5, 0, -12, -4, 3]) # 3x^4 - 4x^3 - 12x^2 + 5
dp, d2p = p.deriv(), p.deriv(2)
print(dp)
for c in dp.roots(): # kritik noktalar
kind = "yerel minimum" if d2p(c) > 0 else "yerel maksimum"
print(f"x = {c:5.2f} p(x) = {p(c):6.2f} p''(x) = {d2p(c):5.1f}"
f" {kind}")Çıktı:
0.0 - 24.0 x - 12.0 x**2 + 12.0 x**3
x = -1.00 p(x) = 0.00 p''(x) = 36.0 yerel minimum
x = 0.00 p(x) = 5.00 p''(x) = -24.0 yerel maksimum
x = 2.00 p(x) = -27.00 p''(x) = 72.0 yerel minimumAdım 3. Türevin çarpanlarına ayrılmış hâli
\[ p'(x) = 12x^3 - 12x^2 - 24x = 12x(x + 1)(x - 2) \]
olduğundan kritik noktalar \(-1\), \(0\) ve \(2\)’dir. \(x = -1\) ve \(x = 2\)’de yerel minimum (\(p(-1) = 0\), \(p(2) = -27\)), \(x = 0\)’da yerel maksimum (\(p(0) = 5\)) vardır. \(\blacksquare\)
Alıştırma 12.2 (Horner Yöntemiyle Sentetik Bölme) Horner yöntemindeki ara değerlerin, \(p(x)\) polinomu \(x - r\)’ye bölündüğünde bölümün katsayıları olduğunu gösteriniz. Bunu kullanarak bölümü ve kalanı veren bir synthetic_division(coef, r) fonksiyonu yazınız ve \((x-1)(x-2)(x-3)\) polinomunu \(x - 4\)’e bölüp sonucu divmod ile karşılaştırınız.
Çözüm
Adım 1. Horner yinelemesini \(x = r\) noktasında yapalım: \(b_n = a_n\) ve \(b_k = a_k + r\,b_{k+1}\). \(q(x) = b_1 + b_2 x + \dots + b_n x^{n-1}\) diyelim. Çarpımı açarsak
\[ (x - r)\,q(x) = b_n x^n + \sum_{k=1}^{n-1} (b_k - r\,b_{k+1})\, x^k - r\,b_1 \]
bulunur. Yineleme gereği \(b_k - r\,b_{k+1} = a_k\) ve \(b_n = a_n\)’dir. Ayrıca \(b_0 - r\,b_1 = a_0\) olduğundan \(-r\,b_1 = a_0 - b_0\)’dır. Dolayısıyla \((x - r)\,q(x) = p(x) - b_0\), yani \(p(x) = (x - r)\,q(x) + b_0\) olur. Bölüm \(q\), kalan da \(b_0 = p(r)\)’dir.
Adım 2. Fonksiyon ara değerleri bir listede biriktirir; son değer kalandır, öncekiler bölümün katsayılarıdır:
from numpy.polynomial import Polynomial
def synthetic_division(coef, r):
"""p(x) = (x - r) q(x) + p(r) ayrışması; katsayılar artan sırada."""
b = [coef[-1]]
for a in reversed(coef[:-1]):
b.append(a + r * b[-1]) # Horner adımı
rem = b.pop() # son değer p(r)
return b[::-1], rem
coef = [-6, 11, -6, 1] # (x - 1)(x - 2)(x - 3)
q, rem = synthetic_division(coef, 4)
print(q, rem)
quot, rem2 = divmod(Polynomial(coef), Polynomial([-4, 1]))
print(quot, "|", rem2)Çıktı:
[3, -2, 1] 6
3.0 - 2.0 x + 1.0 x**2 | 6.0Adım 3. Bölüm \(3 - 2x + x^2\), kalan \(6 = p(4) = 3 \cdot 2 \cdot 1\)’dir; divmod aynı sonucu verdi. Fonksiyon bölümün katsayılarını artan sırada döndürür, bu yüzden liste [3, -2, 1] biçimindedir. \(\blacksquare\)
Alıştırma 12.3 (Beş Katlı Bir Kökün Hesabı) \((x-1)^5\) polinomunun köklerini roots metoduyla hesaplayınız. Bulunan köklerin \(1\)’e uzaklıklarını açıklayınız.
Çözüm
Adım 1. Polinomu kökünden kurup köklerini ve köklerin \(1\)’e uzaklıklarını yazdıralım; karşılaştırma için makine hassasiyeti \(\varepsilon\)’un beşinci kökünü de hesaplayalım:
import numpy as np
from numpy.polynomial import Polynomial
p = Polynomial.fromroots([1, 1, 1, 1, 1]) # (x - 1)^5
print(p)
r = p.roots()
print(np.round(r, 5))
print(np.abs(r - 1))
print(np.finfo(float).eps ** (1 / 5))Çıktı:
-1.0 + 5.0 x - 10.0 x**2 + 10.0 x**3 - 5.0 x**4 + 1.0 x**5
[0.99877+0.j 0.99962-0.00117j 0.99962+0.00117j 1.00099-0.00072j
1.00099+0.00072j]
[0.0012277 0.00122816 0.00122816 0.0012289 0.0012289 ]
0.000740095979741405Adım 2. Kökler \(1\) çevresinde, yarıçapı yaklaşık \(1{,}23 \cdot 10^{-3}\) olan bir çemberin üzerinde, düzgün bir beşgenin köşelerinde duruyor. Hesaplanan kökler, katsayıları yuvarlama hatası düzeyinde değişmiş bir polinomun kökleridir. \((x - 1)^5 = \delta\) denkleminin kökleri \(1\)’in çevresinde, yarıçapı \(|\delta|^{1/5}\) olan bir çemberin üzerinde eşit aralıklarla dizilir. \(|\delta|^{1/5} \approx 1{,}23 \cdot 10^{-3}\) ise \(|\delta| \approx 2{,}8 \cdot 10^{-15}\)’tir; bu, makine hassasiyetinin yaklaşık \(13\) katıdır.
Adım 3. Demek ki \(10^{-15}\) mertebesindeki bir değişiklik, beş katlı kökü \(10^{-3}\) mertebesinde oynatır. Genel olarak \(m\) katlı bir kök, katsayılardaki \(\varepsilon\) büyüklüğündeki bir değişiklikle \(\varepsilon^{1/m}\) mertebesinde yer değiştirir; son satırdaki \(\varepsilon^{1/5} \approx 7{,}4 \cdot 10^{-4}\) bu ölçeği gösteriyor. \(\blacksquare\)
Alıştırma 12.4 (Bölünmüş Farklarla Newton Biçimi) \((0, 1)\), \((1, 3)\), \((2, 2)\), \((3, 5)\) verisinin bölünmüş farklarını hesaplayan bir fonksiyon yazınız. İnterpolasyon polinomunun Newton biçimini Polynomial nesneleriyle kurup Vandermonde yoluyla bulduğumuz polinomla karşılaştırınız.
Çözüm
Adım 1. Newton biçimi
\[ p(x) = f[x_0] + f[x_0, x_1](x - x_0) + \dots + f[x_0, \dots, x_n](x - x_0) \cdots (x - x_{n-1}) \]
yazılışıdır (bkz. Nümerik Analiz). Bölünmüş farklar \(f[x_i] = y_i\) ve
\[ f[x_i, \dots, x_{i+j}] = \frac{f[x_{i+1}, \dots, x_{i+j}] - f[x_i, \dots, x_{i+j-1}]}{x_{i+j} - x_i} \]
yinelemesiyle hesaplanır.
Adım 2. Fark tablosunu tek bir dizi üzerinde sütun sütun güncelleriz: \(j\)’inci adımdan sonra c[i] (\(i \ge j\)) elemanı \(f[x_{i-j}, \dots, x_i]\)’yi tutar. Atamanın sağ tarafı eski değerlerle hesaplandığından güncelleme vektörel yapılabilir. Newton biçimini kurarken \((x - x_0) \cdots (x - x_{k-1})\) çarpımını basis adlı bir Polynomial nesnesinde biriktiririz:
import numpy as np
from numpy.polynomial import Polynomial
def divided_differences(x, y):
"""Newton biçiminin katsayıları: f[x0], f[x0,x1], ..., f[x0,...,xn]."""
c = np.array(y, dtype=float)
for j in range(1, len(x)):
# c[i] <- (c[i] - c[i-1]) / (x[i] - x[i-j]), i = j, ..., n
c[j:] = (c[j:] - c[j - 1:-1]) / (x[j:] - x[:-j])
return c
x = np.array([0.0, 1.0, 2.0, 3.0])
y = np.array([1.0, 3.0, 2.0, 5.0])
c = divided_differences(x, y)
print(c)
p = Polynomial([0.0])
basis = Polynomial([1.0]) # (x - x0)...(x - x_{k-1})
for ck, xk in zip(c, x):
p = p + ck * basis
basis = basis * Polynomial([-xk, 1.0])
print(p)Çıktı:
[ 1. 2. -1.5 1.16666667]
1.0 + 5.83333333 x - 5.0 x**2 + 1.16666667 x**3Adım 3. Bölünmüş farklar \(1\), \(2\), \(-1{,}5\) ve \(\frac{7}{6}\)’dır, yani
\[ p(x) = 1 + 2x - \frac{3}{2}\,x(x - 1) + \frac{7}{6}\,x(x - 1)(x - 2) . \]
Açılmış hâli Vandermonde yoluyla bulduğumuz \(1 + \frac{35}{6}x - 5x^2 + \frac{7}{6}x^3\) polinomunun aynısıdır. Son bölünmüş fark, beklendiği gibi baş katsayıya eşittir (bkz. Nümerik Analiz). \(\blacksquare\)
Alıştırma 12.5 (Sinüs İnterpolasyonunun Hata Sınırı) \([0, \pi]\) aralığındaki \(5\) eşit aralıklı düğümde \(\sin x\) fonksiyonunun interpolasyon polinomunu kurunuz. Gerçek en büyük hatayı interpolasyon hatası teoreminin verdiği sınırla karşılaştırınız.
Çözüm
Adım 1. \(5\) düğüm \(n = 4\) demektir. \(|\sin^{(5)} x| = |\cos x| \le 1\) olduğundan teorem (Teorem 12.2) her \(x \in [0, \pi]\) için
\[ |\sin x - p(x)| \le \frac{\max |\omega|}{5!} \]
sınırını verir.
Adım 2. Beş noktaya dördüncü dereceden Polynomial.fit interpolasyon polinomunu verir. \(\max|\omega|\)’yı sık bir ızgarada ölçeriz:
import math
import numpy as np
from numpy.polynomial import Polynomial
xn = np.linspace(0, np.pi, 5) # n = 4, beş düğüm
p = Polynomial.fit(xn, np.sin(xn), 4)
xx = np.linspace(0, np.pi, 4001)
err = np.max(np.abs(p(xx) - np.sin(xx)))
omega = np.max(np.abs(np.prod(xx[:, None] - xn, axis=1)))
bound = omega / math.factorial(5) # max |sin^(5)| = 1
print(f"gerçek hata: {err:.2e}")
print(f"max |omega|: {omega:.4f} sınır: {bound:.2e}")Çıktı:
gerçek hata: 1.81e-03
max |omega|: 1.0852 sınır: 9.04e-03Adım 3. Gerçek hata \(1{,}81 \cdot 10^{-3}\), sınır \(9{,}04 \cdot 10^{-3}\)’tür; sınır sağlanıyor ve gerçek hatanın yaklaşık \(5\) katı. Sınırın kötümser olmasının nedeni, hata formülündeki \(|\sin^{(5)} \xi| = |\cos \xi|\) çarpanının yerine en kötü değeri olan \(1\)’i koymasıdır. Teoremdeki \(\xi\) noktası \(x\)’e bağlıdır ama \(x\)’in yakınında olmak zorunda değildir; bu yüzden \(|\cos x|\)’in \(1\)’e yakın olduğu uçlarda bile \(|\cos \xi|\) küçük olabilir. Gerçek hatadan geri hesaplanan \(|\cos \xi(x)|\) değeri, yani \(5!\,|\sin x - p(x)| / |\omega(x)|\) oranı, düğümler dışındaki her noktada \(0{,}24\)’ün altında kalıyor. \(\blacksquare\)
Alıştırma 12.6 (Runge’un Özgün Örneği) Runge örneğini \([-5, 5]\) aralığında \(g(x) = \frac{1}{1 + x^2}\) fonksiyonuyla vermişti. \(11\) eşit aralıklı düğümle ve \([-5, 5]\) aralığına taşınmış \(11\) Çebişev düğümüyle en büyük interpolasyon hatalarını hesaplayınız. Sonuçların bölümdeki \(1{,}9156\) ve \(0{,}1092\) değerleriyle aynı çıkmasını açıklayınız.
Çözüm
Adım 1. Çebişev düğümlerini \([-1, 1]\) aralığında hesaplayıp \(t \mapsto \frac{a+b}{2} + \frac{b-a}{2}\,t\) dönüşümüyle \([a, b] = [-5, 5]\) aralığına taşırız (Tanım 12.5). Chebyshev.interpolate aynı işi domain seçeneğiyle yapar.
Adım 2. Hataları sık bir ızgarada ölçelim:
import numpy as np
from numpy.polynomial import Chebyshev
from scipy.interpolate import BarycentricInterpolator
def f(x):
return 1 / (1 + x**2)
a, b, n = -5.0, 5.0, 10
xx = np.linspace(a, b, 4001)
xe = np.linspace(a, b, n + 1)
k = np.arange(n + 1)
t = np.cos((2 * k + 1) * np.pi / (2 * n + 2)) # [-1, 1] düğümleri
xc = (a + b) / 2 + (b - a) / 2 * t # [a, b] aralığına taşı
for name, xn in (("eşit aralıklı", xe), ("Çebişev", xc)):
p = BarycentricInterpolator(xn, f(xn))
print(f"{name:13s}: {np.max(np.abs(p(xx) - f(xx))):.4f}")
c = Chebyshev.interpolate(f, n, domain=[a, b])
print(f"{'interpolate':13s}: {np.max(np.abs(c(xx) - f(xx))):.4f}")Çıktı:
eşit aralıklı: 1.9156
Çebişev : 0.1092
interpolate : 0.1092Adım 3. Hatalar bölümdeki \(f(x) = \frac{1}{1 + 25x^2}\) için bulunanlarla aynı. Nedeni, \(x = 5u\) dönüşümüyle \(g(5u) = \frac{1}{1 + 25u^2} = f(u)\) olmasıdır. Bu dönüşüm \([-1, 1]\)’deki eşit aralıklı düğümleri \([-5, 5]\)’teki eşit aralıklı düğümlere, Çebişev düğümlerini de Çebişev düğümlerine taşır. \(g\)’nin interpolasyon polinomu \(P\) ise \(P(5u)\), \(u\)’nun derecesi en fazla \(10\) olan ve \(f\)’yi \(u\)-düğümlerinde interpolasyonla yaklaşan bir polinomdur; teklik gereği \(f\)’nin interpolasyon polinomudur. Dolayısıyla iki problemin hataları nokta nokta aynıdır. \(\blacksquare\)
Alıştırma 12.7 (Hız Ölçümlerinden Yol ve İvme) Bir aracın hızı \(t = 0, 1, \dots, 6\) saniyelerinde ölçülüyor; ölçümler \(v(t) = 10\,(1 - e^{-t/2})\) m/s fonksiyonundan alınıyor. Kübik spline ile \([0, 6]\) aralığında alınan yolu ve \(t = 3\)’teki ivmeyi tahmin edip gerçek değerlerle karşılaştırınız.
Çözüm
Adım 1. Yol hızın integrali, ivme hızın türevidir. Gerçek değerler
\[ \int_0^6 v(t)\,dt = 10\Big[t + 2e^{-t/2}\Big]_0^6 = 40 + 20e^{-3}, \qquad v'(3) = 5e^{-3/2} \]
olur.
Adım 2. Spline’ı kurup integrate ve birinci türevle hesaplayalım:
import numpy as np
from scipy.interpolate import CubicSpline
t = np.arange(0.0, 7.0) # 0, 1, ..., 6 saniye
v = 10 * (1 - np.exp(-t / 2)) # m/s, ölçülen hızlar
S = CubicSpline(t, v)
print(np.round(v, 4))
print(f"yol: {S.integrate(0, 6):.4f} gerçek: {40 + 20 * np.exp(-3):.4f}")
print(f"ivme: {S(3, 1):.4f} gerçek: {5 * np.exp(-1.5):.4f}")Çıktı:
[0. 3.9347 6.3212 7.7687 8.6466 9.1792 9.5021]
yol: 40.9908 gerçek: 40.9957
ivme: 1.1164 gerçek: 1.1157Adım 3. Yalnız \(7\) ölçümle yol \(40{,}9908\) m bulundu; gerçek değer \(40{,}9957\) m, bağıl hata yaklaşık \(10^{-4}\). İvme \(1{,}1164\) m/s², gerçeği \(1{,}1157\) m/s²’dir. \(\blacksquare\)
Alıştırma 12.8 (Doğal Spline’ın Uçlardaki Hatası) \([0, 1]\) aralığında \(e^x\) fonksiyonunu \(n = 10, 20, 40, 80\) eşit alt aralıkla doğal spline ve düğüm olmama koşullu spline ile yaklaşınız. Hataların azalma hızlarını karşılaştırıp farkı açıklayınız.
Çözüm
Adım 1. İki spline’ı her \(n\) için kurup en büyük hatalarını ölçelim:
import numpy as np
from scipy.interpolate import CubicSpline
xx = np.linspace(0, 1, 20001)
for n in (10, 20, 40, 80):
xn = np.linspace(0, 1, n + 1)
errs = []
for bc in ("natural", "not-a-knot"):
S = CubicSpline(xn, np.exp(xn), bc_type=bc)
errs.append(np.max(np.abs(S(xx) - np.exp(xx))))
print(f"n = {n:2d} doğal: {errs[0]:8.2e} not-a-knot: {errs[1]:8.2e}")Çıktı:
n = 10 doğal: 1.33e-03 not-a-knot: 6.93e-06
n = 20 doğal: 3.34e-04 not-a-knot: 4.56e-07
n = 40 doğal: 8.34e-05 not-a-knot: 2.92e-08
n = 80 doğal: 2.09e-05 not-a-knot: 1.85e-09Adım 2. \(n\) iki katına çıkınca doğal spline’ın hatası \(4\) kat (\(1{,}33 \cdot 10^{-3} / 3{,}34 \cdot 10^{-4} \approx 4{,}0\)), düğüm olmama koşullu spline’ınki yaklaşık \(16\) kat (\(4{,}56 \cdot 10^{-7} / 2{,}92 \cdot 10^{-8} \approx 15{,}6\)) küçülüyor. Doğal spline yalnız \(h^2\), öteki \(h^4\) hızıyla yakınsıyor.
Adım 3. Nedeni uç koşuludur. \((e^x)'' = e^x\), uçlarda \(1\) ve \(e\) değerlerini alır; doğal spline ise uçlarda \(S'' = 0\)’ı zorlar. Bu yanlış bilgi uçlara yakın parçalarda \(h^2\) mertebesinde bir hata doğurur. Gerçek \(f''\) uçlarda sıfır olmadıkça doğal spline’ın bu zayıflığı vardır. \(\blacksquare\)
Alıştırma 12.9 (Doğal Spline’ın Koşullarını Denetlemek) \((0;\ 1)\), \((1;\ 2)\), \((2{,}5;\ 0)\), \((3;\ 1)\) ve \((4;\ 3)\) noktalarından geçen doğal spline’ı kurunuz. S.c katsayılarını kullanarak \(S''\)’nün iç düğümlerde sürekli ve uçlarda sıfır olduğunu doğrulayınız.
Çözüm
Adım 1. \(i\)’inci parçada \(u = x - x_i\) ile \(S(x) = c_0u^3 + c_1u^2 + c_2u + c_3\)’tür, burada \(c_0, c_1\) sırasıyla S.c[0, i] ve S.c[1, i]’dir. İkinci türev \(S'' = 6c_0u + 2c_1\) olur. Parçanın uzunluğu \(h_i\) ise \(S''\) parçanın sağ ucunda \(6c_0h_i + 2c_1\), sol ucunda \(2c_1\) değerini alır. Bir iç düğümde, soldaki parçanın sağ uç değeri ile sağdaki parçanın sol uç değeri eşit olmalıdır.
Adım 2. Bu değerleri bütün parçalar için birden hesaplayalım; kodda \(c_0\) ve \(c_1\), c3 ve c2 adlarıyla tutuluyor:
import numpy as np
from scipy.interpolate import CubicSpline
x = np.array([0.0, 1.0, 2.5, 3.0, 4.0])
y = np.array([1.0, 2.0, 0.0, 1.0, 3.0])
S = CubicSpline(x, y, bc_type="natural")
h = np.diff(x)
c3, c2 = S.c[0], S.c[1] # i'inci parça: c3 u^3 + c2 u^2 + ...
left = 6 * c3 * h + 2 * c2 # S'' parçanın sağ ucunda (soldan)
right = 2 * c2 # S'' parçanın sol ucunda (sağdan)
print(np.round(left[:-1], 10)) # iç düğümlerde soldan
print(np.round(right[1:], 10)) # iç düğümlerde sağdan
print(S(x[0], 2), S(x[-1], 2)) # uç noktalarda S''Çıktı:
[-4.89423077 6.98076923 -1.16346154]
[-4.89423077 6.98076923 -1.16346154]
0.0 0.0Adım 3. İç düğümlerde (\(x = 1;\ 2{,}5;\ 3\)) soldan ve sağdan hesaplanan \(S''\) değerleri on basamağa kadar aynı ve uç noktalarda \(S'' = 0\). Spline, tanımdaki süreklilik koşulunu ve doğal spline’ın uç koşulunu sağlıyor. \(\blacksquare\)
Alıştırma 12.10 (Test Noktalarıyla Derece Seçimi) \([0, 1]\) aralığından rastgele seçilmiş \(40\) noktada \(\sin 2\pi x\) değerlerine standart sapması \(0{,}2\) olan gürültü eklenerek bir veri üretiliyor. Her dördüncü noktayı test için ayırıp kalan \(30\) noktaya \(1\)’den \(10\)’a kadar her dereceden polinom uydurunuz. Uydurmada kullanılan noktalardaki artık kareleri toplamını ve test noktalarındaki hatanın karesel ortalamasını karşılaştırarak bir derece seçiniz.
Çözüm
Adım 1. test adlı mantıksal dizi, test noktalarını işaretler; ~test onun değilidir (Vektörizasyon ve Broadcasting). Uydurma yalnız ~test noktalarıyla yapılır. Test hatası olarak
\[ \text{RMS} = \sqrt{\frac{1}{10} \sum_{\text{test}} \big(y_i - p(x_i)\big)^2} \]
karesel ortalamasını kullanırız.
Adım 2. Bütün dereceleri bir döngüde deneyelim:
import numpy as np
from numpy.polynomial import Polynomial
rng = np.random.default_rng(3)
x = np.sort(rng.uniform(0, 1, 40))
y = np.sin(2 * np.pi * x) + rng.normal(0, 0.2, x.size)
test = np.zeros(x.size, dtype=bool)
test[::4] = True # her dördüncü nokta test için
for deg in range(1, 11):
p = Polynomial.fit(x[~test], y[~test], deg)
rss = np.sum((y[~test] - p(x[~test]))**2)
rms = np.sqrt(np.mean((y[test] - p(x[test]))**2))
print(f"derece {deg:2d} eğitim RSS: {rss:6.3f} test RMS: {rms:.3f}")Çıktı:
derece 1 eğitim RSS: 8.859 test RMS: 0.565
derece 2 eğitim RSS: 8.511 test RMS: 0.652
derece 3 eğitim RSS: 1.294 test RMS: 0.184
derece 4 eğitim RSS: 1.212 test RMS: 0.204
derece 5 eğitim RSS: 0.927 test RMS: 0.284
derece 6 eğitim RSS: 0.927 test RMS: 0.285
derece 7 eğitim RSS: 0.902 test RMS: 0.325
derece 8 eğitim RSS: 0.899 test RMS: 0.303
derece 9 eğitim RSS: 0.895 test RMS: 0.281
derece 10 eğitim RSS: 0.889 test RMS: 0.251Adım 3. Eğitim noktalarındaki artık kareleri toplamı dereceyle hep azalıyor, uyarıda gördüğümüz gibi. Test hatası ise \(3\)’üncü derecede en küçük değerini (\(0{,}184\)) alıyor ve daha yüksek derecelerde hep bundan büyük kalıyor; polinom gürültüye uymaya başlıyor. En küçük test hatasının gürültünün \(0{,}2\) olan standart sapmasına yakın olması, modelin açıklanabilen yapıyı yakaladığını düşündürür. Seçilecek derece \(3\)’tür. \(\blacksquare\)
Alıştırma 12.11 (Üstel Bir Modeli İki Yoldan Uydurmak) \(x \in [0, 4]\) aralığında eşit aralıklı \(17\) noktada \(y = 2e^{0{,}8x}\) değerlerine standart sapması \(0{,}5\) olan gürültü ekleniyor. \(y = ae^{bx}\) modelini önce \(\ln y = \ln a + bx\) doğrusuna polyfit uydurarak, sonra curve_fit ile uydurunuz. İki sonucu artık kareleri toplamıyla karşılaştırınız.
Çözüm
Adım 1. Logaritma almak modeli \(\ln a\) ve \(b\) parametrelerine göre doğrusal yapar; ama bunun için bütün ölçümlerin pozitif olması gerekir. curve_fit için başlangıç tahminini logaritma yönteminin sonucundan alırız.
Adım 2. İki uydurmayı yapıp artık kareleri toplamlarını özgün ölçekte hesaplayalım:
import numpy as np
from scipy.optimize import curve_fit
def model(x, a, b):
return a * np.exp(b * x)
rng = np.random.default_rng(5)
x = np.linspace(0, 4, 17)
y = model(x, 2.0, 0.8) + rng.normal(0, 0.5, x.size)
print(f"en küçük ölçüm: {y.min():.3f}") # log için y > 0 gerekir
b_log, ln_a = np.polyfit(x, np.log(y), 1) # ln y = ln a + b x
a_log = np.exp(ln_a)
(a_cf, b_cf), _ = curve_fit(model, x, y, p0=(a_log, b_log))
for name, a, b in (("log + polyfit", a_log, b_log),
("curve_fit", a_cf, b_cf)):
rss = np.sum((y - model(x, a, b))**2)
print(f"{name:13s} a = {a:.4f} b = {b:.4f} RSS = {rss:.2f}")Çıktı:
en küçük ölçüm: 1.599
log + polyfit a = 1.8214 b = 0.8343 RSS = 14.61
curve_fit a = 2.0259 b = 0.7954 RSS = 3.93Adım 3. Logaritma yöntemi \(a = 1{,}82\), \(b = 0{,}834\) ve \(14{,}61\) artık kareleri toplamı verdi; curve_fit ise \(a = 2{,}026\), \(b = 0{,}795\) ve \(3{,}93\) verdi. Fark, iki yöntemin farklı şeyleri en küçük yapmasından gelir: logaritma yöntemi \(\sum (\ln y_i - \ln \hat y_i)^2\)’yi, curve_fit ise \(\sum (y_i - \hat y_i)^2\)’yi en küçük yapar. Gürültü her noktada aynı büyüklükte olduğundan küçük \(y\) değerlerinde göreli olarak büyüktür (\(x = 0\) civarında \(2\)’nin yanında \(0{,}5\)). Logaritma bu noktaların sapmalarını büyütür ve uydurmayı onlara doğru çeker. Asıl ölçekteki hatayı en küçük yapmak istiyorsak curve_fit kullanılmalıdır; logaritma yöntemi yine de iyi bir başlangıç tahmini verir. \(\blacksquare\)
Bu bölümde polinomlarla NumPy’da hesap yapmayı, köklerin katsayılara duyarlılığını, interpolasyonu ve Runge olayı gibi tuzaklarını, spline’larla parça parça yaklaşımı ve ölçüm verisine model uydurmayı gördük. Bu fikirler bir sonraki bölümde de iş görecek: adi diferansiyel denklemler için kullanılan sayısal yöntemler, çözümü küçük adımlarda polinomlarla yaklaşır ve SciPy’ın çözücüleri adımlar arasındaki değerleri de interpolasyonla verir. Diferansiyel Denklemlerin Sayısal Çözümü bölümünde Euler ve Runge–Kutta yöntemlerini kendimiz yazıp solve_ivp ile tanışacağız.