13 Diferansiyel Denklemlerin Sayısal Çözümü
Bir nüfusun büyümesi, bir sarkacın salınımı, iki türün birbirini beslemesi birer diferansiyel denklemle anlatılır. Diferansiyel Denklemler notlarında ayrılabilen, homojen ve tam denklemleri kapalı biçimde çözdük; SymPy ile Sembolik Hesap bölümünde de dsolve aynı işi bilgisayara yaptırdı. Ne var ki çözümü bir formülle yazılabilen denklemler azınlıktadır. \(y' = t - y^2\) gibi masum görünen bir denklemin çözümü elementer fonksiyonlarla yazılamaz; sarkacın \(\theta'' = -\sin\theta\) denkleminin çözümü ise eliptik fonksiyonlar adı verilen özel fonksiyonlar gerektirir.
Böyle durumlarda çözümün formülünü değil, sonlu sayıda noktadaki değerlerini hesaplarız. Bu bölümde önce bu fikrin en basit iki uygulamasını, Euler ve dördüncü mertebe Runge–Kutta yöntemlerini, birkaç satırlık Python fonksiyonları olarak kendimiz yazacak ve hatalarının adım uzunluğuyla nasıl küçüldüğünü ölçeceğiz. Ardından aynı işi çok daha akıllıca yapan SciPy fonksiyonu solve_ivp’yi, katı denklemleri, denklem sistemlerini ve faz düzlemini göreceğiz. Bölüm, Lotka–Volterra av-avcı modeli ve sarkaçla bitiyor.
Bağımsız değişkeni SciPy’nin alışkanlığına uyarak \(t\) ile göstereceğiz ve onu zaman olarak düşüneceğiz; Diferansiyel Denklemler notlarındaki \(x\)’in yerini tutar. Kodlar NumPy, SciPy ve bir yerde SymPy kullanır; grafik kodları Matplotlib ile Grafik Çizimi bölümündeki nesne yönelimli yazımı izler.
13.1 Başlangıç Değer Problemi ve Yön Alanı
Sayısal yöntemlerin çözdüğü problem, bir diferansiyel denklemle bir başlangıç koşulunun ikilisidir. Birinci mertebeden bir başlangıç değer problemi (bkz. Diferansiyel Denklemler)
\[y' = f(t, y), \qquad y(t_0) = y_0\]
biçimindedir. \(f\) ve \(\partial f/\partial y\), \((t_0, y_0)\) noktasını içeren bir bölgede sürekliyse problemin \(t_0\)’ın bir komşuluğunda tek bir çözümü vardır (bkz. Diferansiyel Denklemler). Sayısal yöntemler bu çözümü bulmaz, ona yaklaşır; varlık ve teklik teoremi de yaklaştığımız şeyin gerçekten var ve tek olduğunu güvenceye alır. Bu bölümdeki denklemlerin sağ tarafları polinom, üstel ve trigonometrik fonksiyonlardan kurulu olduğundan bu koşullar her yerde sağlanır.
Denklemi çözmeden önce bile \(f\)’nin bize bir şey söylediğini fark edebiliriz: her noktada çözümün eğimini.
Tanım 13.1 (Yön Alanı) \(y' = f(t, y)\) denklemi verilsin. \(ty\)-düzleminin her \((t, y)\) noktasına, o noktadan geçen ve eğimi \(f(t, y)\) olan kısa bir doğru parçası iliştirelim. Bu doğru parçalarının oluşturduğu resme denklemin yön alanı (direction field) denir.
Yani \((t, y)\) noktasından geçen çözüm eğrisi o noktada tam olarak \(f(t, y)\) eğimiyle ilerler; çözüm eğrileri bu yüzden yön alanındaki doğru parçalarına her noktada teğettir. Başlangıç değer problemini çözmek, \((t_0, y_0)\) noktasından çıkıp bu doğru parçalarının gösterdiği yönü izlemek demektir.
Yön alanını çizmek için \(f\)’yi bir ızgaranın her noktasında hesaplamak yeter. Vektörizasyon ve Broadcasting bölümündeki np.meshgrid bunu tek satırda yapar. \(y' = t - y^2\) denkleminin eğimlerini küçük bir ızgarada hesaplayalım:
import numpy as np
def f(t, y):
return t - y**2
t = np.linspace(0, 2, 5) # 0; 0.5; 1; 1.5; 2
y = np.linspace(1, -1, 5) # yukarıdan aşağıya: 1; 0.5; ...; -1
T, Y = np.meshgrid(t, y)
S = f(T, Y) # her ızgara noktasındaki eğim
print("t :", t)
for yk, row in zip(y, S):
print(f"y = {yk:4.1f}:", row)Çıktı:
t : [0. 0.5 1. 1.5 2. ]
y = 1.0: [-1. -0.5 0. 0.5 1. ]
y = 0.5: [-0.25 0.25 0.75 1.25 1.75]
y = 0.0: [0. 0.5 1. 1.5 2. ]
y = -0.5: [-0.25 0.25 0.75 1.25 1.75]
y = -1.0: [-1. -0.5 0. 0.5 1. ]Her satır sabit bir \(y\) değerine, her sütun sabit bir \(t\) değerine karşılık gelir; satırları, şekildeki gibi yukarıdan aşağıya dizdik. Eğim \(y\)’ye yalnız \(y^2\) üzerinden bağlı olduğundan \(y = 1\) ile \(y = -1\) satırları ve \(y = 0{,}5\) ile \(y = -0{,}5\) satırları aynıdır. \(t = y^2\) olan noktalarda, örneğin \((1, 1)\) ve \((1, -1)\)’de, eğim sıfırdır. Bu parabolün iç tarafında (\(t > y^2\)) eğim pozitiftir ve çözümler artar; dış tarafında çözümler azalır. Aşağıdaki şekil aynı hesabı daha sık bir ızgarada yapıp her noktaya eğimi kadar eğik bir çizgi koyar:
Şekil, denklemi çözmeden çok şey söylüyor. \(y(0) \ge -0{,}6\) olan çözümler \(t\) büyüdükçe birbirine yaklaşıp parabolün üst kolunun hemen altında ilerliyor. \(y(0) \le -0{,}8\) olanlar ise parabolün alt kolunun dışına düşüyor ve orada eğim gittikçe dikleşen negatif değerler aldığından hızla aşağı iniyor. Kalın eğri \(y(0) = 0\) koşulunu sağlayan çözümdür; onun \(t = 2\)’deki değerini bu bölümde birkaç yöntemle hesaplayacağız.
13.2 Euler Yöntemi
Yön alanını izlemenin en doğrudan yolu, bulunduğumuz noktadaki eğimle kısa bir adım atmak ve vardığımız yerde eğimi yeniden okumaktır.
Tanım 13.2 (Euler Yöntemi) \(y' = f(t, y)\), \(y(t_0) = y_0\) problemi ve bir \(h > 0\) adım uzunluğu (step size) verilsin. \(t_i = t_0 + ih\) noktalarında çözüme yaklaşan \(y_i\) sayıları
\[y_{i+1} = y_i + h\, f(t_i, y_i), \qquad i = 0, 1, 2, \dots\]
kuralıyla hesaplanır. Bu kurala Euler yöntemi denir.
Yani her adımda çözümü bulunduğumuz noktadaki teğet doğrusuyla değiştiririz: \((t_i, y_i)\) noktasından eğimi \(f(t_i, y_i)\) olan doğru boyunca \(h\) kadar ilerleyip \((t_{i+1}, y_{i+1})\) noktasına varırız. Formül Taylor açılımından da okunabilir:
\[y(t + h) = y(t) + h\,y'(t) + \frac{h^2}{2}\,y''(\xi)\]
eşitliğinde son terimi atıp \(y'(t) = f\bigl(t, y(t)\bigr)\) yazarsak Euler adımı çıkar (bkz. Nümerik Analiz). Atılan \(\frac{h^2}{2}\,y''(\xi)\) terimi, bir adımda yapılan hatadır.
Yöntemi bir Python fonksiyonu olarak yazalım. Fonksiyon \(f\)’yi, başlangıç noktasını, adım uzunluğunu ve adım sayısını alır; \(t_i\) ve \(y_i\) değerlerini iki NumPy dizisi olarak döndürür. İlk denemeyi, çözümü \(y = e^t\) olan \(y' = y\), \(y(0) = 1\) problemiyle \(h = 0{,}25\) alarak yapalım:
import numpy as np
def euler(f, t0, y0, h, n):
"""y' = f(t, y), y(t0) = y0 problemini n Euler adımıyla çözer."""
t = t0 + h * np.arange(n + 1) # t_0, t_1, ..., t_n
y = np.zeros(n + 1)
y[0] = y0
for i in range(n):
y[i + 1] = y[i] + h * f(t[i], y[i])
return t, y
t, y = euler(lambda t, y: y, 0.0, 1.0, 0.25, 4)
print(" t y_n e^t hata")
for tk, yk in zip(t, y):
print(f"{tk:.2f} {yk:.6f} {np.exp(tk):.6f} {np.exp(tk) - yk:.6f}")Çıktı:
t y_n e^t hata
0.00 1.000000 1.000000 0.000000
0.25 1.250000 1.284025 0.034025
0.50 1.562500 1.648721 0.086221
0.75 1.953125 2.117000 0.163875
1.00 2.441406 2.718282 0.276876Hesap elle de kolayca izlenir. \(f(t, y) = y\) olduğundan her adım \(y_{i+1} = y_i + 0{,}25\,y_i = 1{,}25\,y_i\) olur, yani \(y_i = 1{,}25^i\)’dir. \(t = 1\)’de \(1{,}25^4 = 2{,}441406\ldots\) değeri \(e = 2{,}718281\ldots\)’den \(0{,}277\) kadar küçük kalır ve hata her adımda büyür. Nedenini şekil gösteriyor:
\(e^t\) konveks olduğundan teğet doğruları eğrinin altında kalır. Her Euler adımı bu yüzden bir alttaki çözüm eğrisine iner ve sonraki adım artık o eğrinin teğetini izler; küçük hatalar böylece birikir.
- Denklemin sağ tarafını
f(t, y)biçiminde bir Python fonksiyonu olarak yazın. - Başlangıç noktası \((t_0, y_0)\)’ı, varmak istediğiniz \(T\) noktasını ve adım sayısı \(n\)’yi seçin; adım uzunluğu \(h = (T - t_0)/n\) olur.
euler(f, t0, y0, h, n)çağrısıyla \(y_1, \dots, y_n\) değerlerini hesaplayın.- Sonucu denetleyin: tam çözüm biliniyorsa onunla karşılaştırın; bilinmiyorsa \(n\)’yi ikiye katlayıp sonucun ne kadar değiştiğine bakın.
Örnek 13.1 (Tam Çözümü Bilinmeyen Bir Problemde Euler Yöntemi) \(y' = t - y^2\), \(y(0) = 0\) probleminde \(y(2)\)’yi Euler yöntemiyle \(n = 25, 50, 100, 200, 400\) adımla hesaplayınız. Sonucun kaç basamağına güvenilebileceğini ardışık sonuçlara bakarak kestiriniz.
Çözüm
Adım 1 ve 2. \(f(t, y) = t - y^2\), \(t_0 = 0\), \(y_0 = 0\), \(T = 2\) ve \(h = 2/n\)’dir.
Adım 3. Yukarıda tanımladığımız euler fonksiyonunu kullanıyoruz; kod bir önceki kodun devamıdır:
def f(t, y):
return t - y**2
ends = []
for n in [25, 50, 100, 200, 400]:
t, y = euler(f, 0.0, 0.0, 2 / n, n)
ends.append(y[-1])
print(f"n = {n:3d} h = {2 / n:.3f} y(2) ≈ {y[-1]:.6f}")
d = np.diff(ends)
print("ardışık farklar:", np.round(d, 6))
print("farkların oranı:", np.round(d[:-1] / d[1:], 3))Çıktı:
n = 25 h = 0.080 y(2) ≈ 1.194340
n = 50 h = 0.040 y(2) ≈ 1.193908
n = 100 h = 0.020 y(2) ≈ 1.193731
n = 200 h = 0.010 y(2) ≈ 1.193651
n = 400 h = 0.005 y(2) ≈ 1.193613
ardışık farklar: [-4.31e-04 -1.78e-04 -8.00e-05 -3.80e-05]
farkların oranı: [2.429 2.214 2.107]Adım 4. Tam çözüm bilinmediğinden ardışık sonuçların farklarına bakıyoruz. Adım yarıya indikçe fark da kabaca yarıya iniyor; farkların oranı \(2\)’ye yaklaşıyor. Hata adımla orantılıysa, yani \(y_h \approx y(2) + Ch\) ise
\[y_h - y_{h/2} \approx Ch - C\,\frac{h}{2} = C\,\frac{h}{2} \approx y_{h/2} - y(2)\]
olur: son farkın büyüklüğü, son sonucun hatası için iyi bir kestirimdir. \(h = 0{,}005\) ile bulunan \(1{,}193613\) değerinin hatası bu yüzden yaklaşık \(3{,}8 \cdot 10^{-5}\)’tir ve güvenebileceğimiz kısım \(1{,}1936\)’dır. Hatta kestirilen hatayı sonuçtan çıkarırsak çok daha iyi bir değer elde ederiz: \(1{,}193613 - 0{,}000038 = 1{,}193575\). Bu sayıları aşağıda solve_ivp ile çok daha doğru bir hesapla karşılaştıracağız.
\(\blacksquare\)
Euler yönteminin hatasının adımla orantılı olduğunu tam çözümü bilinen bir problemde kesin olarak görebiliriz.
Örnek 13.2 (Euler Yöntemi ve e Sayısı) \(y' = y\), \(y(0) = 1\) problemine \([0, 1]\) aralığında \(n\) adımlık Euler yöntemi uygulandığında \(y_n = \left(1 + \frac{1}{n}\right)^n\) olduğunu gösteriniz ve bunu \(n = 10, 100, 1000, 10000\) için sayısal olarak doğrulayınız.
Çözüm
Formül. \(h = 1/n\) ve \(f(t, y) = y\) için Euler adımı
\[y_{i+1} = y_i + \frac{1}{n}\,y_i = \left(1 + \frac{1}{n}\right) y_i\]
olur. \(y_0 = 1\) olduğundan tümevarımla \(y_i = \left(1 + \frac{1}{n}\right)^i\), özel olarak \(y_n = \left(1 + \frac{1}{n}\right)^n\) bulunur. Bu dizinin limiti \(e\)’dir (bkz. Analiz 1); yani Euler yöntemi \(h \to 0\) iken doğru değere yakınsar.
Sayısal doğrulama. Yine euler fonksiyonunu kullanıyoruz; kod önceki kodların devamıdır:
print(" n Euler y(1) (1 + 1/n)^n n·(e - y(1))")
for n in [10, 100, 1000, 10000]:
t, y = euler(lambda t, y: y, 0.0, 1.0, 1 / n, n)
print(f"{n:6d} {y[-1]:.10f} {(1 + 1 / n)**n:.10f} "
f"{n * (np.e - y[-1]):.6f}")
print("e / 2 =", np.e / 2)Çıktı:
n Euler y(1) (1 + 1/n)^n n·(e - y(1))
10 2.5937424601 2.5937424601 1.245394
100 2.7048138294 2.7048138294 1.346800
1000 2.7169239322 2.7169239322 1.357896
10000 2.7181459268 2.7181459268 1.359016
e / 2 = 1.3591409142295225İlk iki sütun aynıdır. Son sütundaki \(n\,(e - y_n)\) çarpımı \(e/2 \approx 1{,}3591\)’e yaklaşıyor; yani hata yaklaşık \(\frac{e}{2n} = \frac{e}{2}\,h\)’dir ve adım uzunluğuyla orantılıdır. Aynı yaklaşıklığı Python ile İlk Adımlar bölümündeki örnekte de görmüştük.
\(\blacksquare\)
13.3 Dördüncü Mertebe Runge–Kutta Yöntemi
Euler yönteminin zayıflığı, eğimi yalnız adımın başında okumasıdır; adım boyunca eğim değiştikçe hata birikir. Runge–Kutta yöntemleri eğimi adımın içinde birkaç noktada ölçer ve bu ölçümlerin ağırlıklı ortalamasıyla ilerler.
Tanım 13.3 (Dördüncü Mertebe Runge–Kutta Yöntemi) \(y' = f(t, y)\), \(y(t_0) = y_0\) problemi ve \(h > 0\) verilsin. \(t_i = t_0 + ih\) olmak üzere her adımda
\[ \begin{aligned} k_1 &= f(t_i,\ y_i),\\[1mm] k_2 &= f\left(t_i + \tfrac{h}{2},\ y_i + \tfrac{h}{2}\,k_1\right),\\[1mm] k_3 &= f\left(t_i + \tfrac{h}{2},\ y_i + \tfrac{h}{2}\,k_2\right),\\[1mm] k_4 &= f\left(t_i + h,\ y_i + h\,k_3\right) \end{aligned} \]
eğimleri hesaplanır ve
\[y_{i+1} = y_i + \frac{h}{6}\,\bigl(k_1 + 2k_2 + 2k_3 + k_4\bigr)\]
alınır. Bu kurala dördüncü mertebe Runge–Kutta yöntemi ya da kısaca RK4 denir.
Yani RK4 eğimi dört kez ölçer: adımın başında (\(k_1\)), ortasında iki kez (\(k_2\) ve \(k_3\)) ve sonunda (\(k_4\)). Ortadaki ve sondaki ölçümler için gereken \(y\) değerleri, bir önceki eğimle atılan kısa Euler adımlarıdır. Sonra bu eğimlerin \(1, 2, 2, 1\) ağırlıklı ortalamasıyla tek bir adım atılır. Ağırlıkların kaynağı Simpson kuralıdır: \(f\) yalnız \(t\)’ye bağlıysa \(k_2 = k_3\) olur ve adım
\[y_{i+1} - y_i = \frac{h}{6}\left(f(t_i) + 4f\left(t_i + \tfrac{h}{2}\right) + f(t_i + h)\right)\]
biçimini alır. Bu, \(\int_{t_i}^{t_i + h} f(t)\,dt\) integrali için Simpson kuralıdır (bkz. Nümerik Analiz); kuralın Python’daki kullanımı Sayısal Türev ve İntegral bölümündedir.
Tek bir adımı elle izleyelim. \(y' = y\), \(y(0) = 1\) için \(h = 1\) alırsak \(f(t, y) = y\) olduğundan
\[ \begin{aligned} k_1 &= f(0;\ 1) = 1,\\[1mm] k_2 &= f(0{,}5;\ 1 + 0{,}5 \cdot 1) = 1{,}5,\\[1mm] k_3 &= f(0{,}5;\ 1 + 0{,}5 \cdot 1{,}5) = 1{,}75,\\[1mm] k_4 &= f(1;\ 1 + 1 \cdot 1{,}75) = 2{,}75 \end{aligned} \]
ve
\[y_1 = 1 + \frac{1}{6}\,(1 + 3 + 3{,}5 + 2{,}75) = 1 + \frac{10{,}25}{6} = 2{,}708\overline{3}\]
bulunur. Tek adımda Euler yöntemi \(2\)’ye varırken RK4, \(e = 2{,}71828\ldots\) sayısına \(0{,}01\)’den daha yakın bir değer verir:
RK4’ü de bir fonksiyon olarak yazalım. Gövde Euler’inkiyle aynıdır; yalnız döngünün içi dört eğim hesaplar. Aynı problemi \(h = 0{,}25\) ile çözelim:
import numpy as np
def rk4(f, t0, y0, h, n):
"""y' = f(t, y), y(t0) = y0 problemini n RK4 adımıyla çözer."""
t = t0 + h * np.arange(n + 1)
y = np.zeros(n + 1)
y[0] = y0
for i in range(n):
k1 = f(t[i], y[i])
k2 = f(t[i] + h / 2, y[i] + h / 2 * k1)
k3 = f(t[i] + h / 2, y[i] + h / 2 * k2)
k4 = f(t[i] + h, y[i] + h * k3)
y[i + 1] = y[i] + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
return t, y
t, y = rk4(lambda t, y: y, 0.0, 1.0, 0.25, 4)
print(" t y_n e^t hata")
for tk, yk in zip(t, y):
print(f"{tk:.2f} {yk:.8f} {np.exp(tk):.8f} {np.exp(tk) - yk:.2e}")Çıktı:
t y_n e^t hata
0.00 1.00000000 1.00000000 0.00e+00
0.25 1.28401693 1.28402542 8.49e-06
0.50 1.64869947 1.64872127 2.18e-05
0.75 2.11695803 2.11700002 4.20e-05
1.00 2.71820994 2.71828183 7.19e-05Aynı dört adımda Euler’in \(0{,}277\) olan hatası RK4’te \(7{,}19 \cdot 10^{-5}\)’e iner; yaklaşık \(3850\) kat küçüktür. RK4 adım başına \(f\)’yi dört kez hesapladığından adil bir karşılaştırma için Euler’e dört kat, yani \(16\) adım vermek gerekir. O zaman bile Euler’in hatası \(e - \left(\frac{17}{16}\right)^{16} \approx 0{,}080\) olur.
Bu farkın nedenini görmek için tek bir RK4 adımını SymPy ile Sembolik Hesap bölümündeki gibi sembolik olarak açalım:
import sympy as sp
h, y = sp.symbols("h y")
def f(t, y):
return y
k1 = f(0, y)
k2 = f(h / 2, y + h / 2 * k1)
k3 = f(h / 2, y + h / 2 * k2)
k4 = f(h, y + h * k3)
step = sp.expand(y + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4))
print(sp.collect(step, y))
print(sp.series(sp.exp(h), h, 0, 6))Çıktı:
y*(h**4/24 + h**3/6 + h**2/2 + h + 1)
1 + h + h**2/2 + h**3/6 + h**4/24 + h**5/120 + O(h**6)\(y' = y\) denkleminde bir RK4 adımı \(y_i\)’yi \(1 + h + \frac{h^2}{2} + \frac{h^3}{6} + \frac{h^4}{24}\) ile çarpar. Tam çözüm ise bir adımda \(e^h\) ile çarpılır; RK4’ün çarpanı tam olarak \(e^h\)’nin dördüncü dereceden Taylor polinomudur. Adım başına hata bu yüzden \(\frac{h^5}{120}\) mertebesindedir. Euler’in çarpanı \(1 + h\) ise yalnız birinci dereceden Taylor polinomudur ve adım başına \(\frac{h^2}{2}\) mertebesinde hata yapar.
13.4 Adım Uzunluğu ve Hata
İki yöntemi tek bir adım uzunluğunda karşılaştırmak yetmez; adım küçüldükçe hatanın nasıl küçüldüğünü ölçmek gerekir. Bundan sonra iki fonksiyonu sık kullanacağımız için önce onları bir modülde toplayalım.
Sınıflar, Hata Yakalama ve Dosyalar bölümünde gördüğümüz gibi aşağıdaki kodu odetools.py adıyla kaydedelim. Fonksiyonlarda küçük bir değişiklik var: np.zeros(n + 1) yerine np.zeros((n + 1,) + np.shape(y0)) yazdık. y0 bir sayıysa np.shape(y0) boş demettir ve hiçbir şey değişmez; y0 \(m\) elemanlı bir diziyse y tablosu \((n + 1) \times m\) boyutlu olur. Böylece aynı fonksiyonlar ileride denklem sistemlerinde de çalışacak.
"""Başlangıç değer problemleri için Euler ve RK4 yöntemleri."""
import numpy as np
def euler(f, t0, y0, h, n):
"""y' = f(t, y), y(t0) = y0 problemini n Euler adımıyla çözer."""
t = t0 + h * np.arange(n + 1)
y = np.zeros((n + 1,) + np.shape(y0))
y[0] = y0
for i in range(n):
y[i + 1] = y[i] + h * f(t[i], y[i])
return t, y
def rk4(f, t0, y0, h, n):
"""y' = f(t, y), y(t0) = y0 problemini n RK4 adımıyla çözer."""
t = t0 + h * np.arange(n + 1)
y = np.zeros((n + 1,) + np.shape(y0))
y[0] = y0
for i in range(n):
k1 = f(t[i], y[i])
k2 = f(t[i] + h / 2, y[i] + h / 2 * k1)
k3 = f(t[i] + h / 2, y[i] + h / 2 * k2)
k4 = f(t[i] + h, y[i] + h * k3)
y[i + 1] = y[i] + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
return t, y
if __name__ == "__main__":
# sınama: y' = y, y(0) = 1 probleminde y(1) = e olmalı
for method in (euler, rk4):
t, y = method(lambda t, y: y, 0.0, 1.0, 0.25, 4)
print(f"{method.__name__:5s} y(1) ≈ {y[-1]:.8f}"
f" hata = {np.e - y[-1]:.2e}")Dosya doğrudan çalıştırıldığında (python odetools.py) sınama bloğu şunu basar:
Çıktı:
euler y(1) ≈ 2.44140625 hata = 2.77e-01
rk4 y(1) ≈ 2.71820994 hata = 7.19e-05Bundan sonra bu fonksiyonları kullanan kodlar from odetools import euler, rk4 satırını içerir; odetools.py dosyası çalıştırılan programla aynı klasörde olmalıdır. Şimdi hatayı adım uzunluğuyla ilişkilendiren kavramları tanımlayalım.
Tanım 13.4 (Global Hata ve Yöntemin Mertebesi) Bir yöntem \([t_0, T]\) aralığında \(n\) adımla, yani \(h = (T - t_0)/n\) adım uzunluğuyla uygulansın. \(E(h) = |y(T) - y_n|\) farkına yöntemin \(T\)’deki global hatası denir. \(f\) yeterince düzgün olduğunda \(h\)’den bağımsız bir \(C\) sabiti ve yeterince küçük her \(h\) için \(E(h) \le C h^p\) oluyorsa, yönteme \(p\)-inci mertebeden denir.
Yani mertebe, adımı küçültmenin ne kadar kazandırdığını söyler: \(p\)-inci mertebeden bir yöntemde adımı yarıya indirmek hatayı kabaca \(2^p\)’ye böler. Global hata, her adımda yapılan küçük hataların birikmiş sonucudur. Adım başına hata \(h^{p+1}\) mertebesindeyse ve \((T - t_0)/h\) adım atılıyorsa, toplam hata \(h^p\) mertebesinde kalır. Euler’de adım başına hata \(h^2\), RK4’te \(h^5\) mertebesinde olduğundan şu sonuç beklenir:
Teorem 13.1 (Euler ve RK4 Yöntemlerinin Mertebesi) \(f\) fonksiyonu, çözüm eğrisini içeren bir bölgede dördüncü mertebeye kadar sürekli kısmi türevlere sahip olsun. Bu durumda \([t_0, T]\) aralığında Euler yöntemi birinci mertebeden, RK4 dördüncü mertebedendir: bir \(C\) sabiti ve yeterince küçük her \(h\) için sırasıyla \(E(h) \le Ch\) ve \(E(h) \le Ch^4\) olur.
İspat, adım başına yapılan hatayı Taylor teoremiyle kestirmeye ve bu hataların sonraki adımlarda en fazla ne kadar büyüyebileceğini sınırlamaya dayanır; Euler’in adım başına hatasının \(\frac{h^2}{2}\,y''(\xi)\) olduğunu yukarıda gördük. Biz teoremi sayısal olarak doğrulayacağız. Bunun için tam çözümü bilinen, doğrusal olmayan bir problem seçelim: lojistik denklem \(y' = y(1 - y)\), \(y(0) = 0{,}1\). Bu ayrılabilen denklemin çözümü (bkz. Diferansiyel Denklemler)
\[y(t) = \frac{1}{1 + 9e^{-t}}\]
dir. \(T = 4\)’teki hatayı \(n = 8, 16, \ldots, 1024\) adımla hesaplayalım:
import numpy as np
from odetools import euler, rk4
def f(t, y):
return y * (1 - y)
def exact(t):
return 1 / (1 + 9 * np.exp(-t))
T = 4.0
ns = 2 ** np.arange(3, 11) # 8, 16, ..., 1024
hs = T / ns
err_e = np.array([abs(euler(f, 0, 0.1, T / n, n)[1][-1] - exact(T))
for n in ns])
err_r = np.array([abs(rk4(f, 0, 0.1, T / n, n)[1][-1] - exact(T))
for n in ns])
print(" n h Euler hata oran RK4 hata oran")
for k in range(len(ns)):
q_e = f"{err_e[k - 1] / err_e[k]:5.2f}" if k else " -"
q_r = f"{err_r[k - 1] / err_r[k]:5.2f}" if k else " -"
print(f"{ns[k]:5d} {hs[k]:.5f} {err_e[k]:.3e} {q_e}"
f" {err_r[k]:.3e} {q_r}")Çıktı:
n h Euler hata oran RK4 hata oran
8 0.50000 1.349e-02 - 5.828e-05 -
16 0.25000 5.489e-03 2.46 3.730e-06 15.62
32 0.12500 2.496e-03 2.20 2.365e-07 15.77
64 0.06250 1.192e-03 2.09 1.490e-08 15.88
128 0.03125 5.825e-04 2.05 9.348e-10 15.94
256 0.01562 2.880e-04 2.02 5.854e-11 15.97
512 0.00781 1.432e-04 2.01 3.663e-12 15.98
1024 0.00391 7.138e-05 2.01 2.286e-13 16.03“oran” sütunu, bir önceki satırın hatasının bu satırın hatasına bölümüdür. Euler’de oran \(2\)’ye, RK4’te \(16 = 2^4\)’e yaklaşıyor; teoremin söylediği tam olarak budur. \(n = 1024\)’te RK4’ün hatası \(2{,}3 \cdot 10^{-13}\)’tür. Adım daha da küçültülürse hata, kayan noktalı sayılarla hesap yapmanın getirdiği yuvarlama hatası düzeyine iner ve oradan aşağı inmez (bkz. Nümerik Analiz).
Mertebe en iyi log-log grafikte görünür. \(E \approx Ch^p\) ise \(\log E \approx \log C + p\log h\) olur ve noktalar eğimi \(p\) olan bir doğru üzerinde dizilir. Eğimi en küçük kareler doğrusuyla (np.polyfit, bkz. Polinomlar, İnterpolasyon ve Eğri Uydurma) ölçelim ve grafiği loglog ile çizelim. Kod bir önceki kodun devamıdır:
import matplotlib.pyplot as plt
p_e = np.polyfit(np.log(hs), np.log(err_e), 1)[0]
p_r = np.polyfit(np.log(hs), np.log(err_r), 1)[0]
print(f"Euler doğrusunun eğimi: {p_e:.3f}")
print(f"RK4 doğrusunun eğimi : {p_r:.3f}")
fig, ax = plt.subplots(figsize=(6, 4.5))
ax.loglog(hs, err_e, "o-", label="Euler")
ax.loglog(hs, err_r, "s-", label="RK4")
ax.loglog(hs, 0.05 * hs, "k--", lw=0.8, label="eğim 1")
ax.loglog(hs, 0.002 * hs**4, "k:", lw=0.8, label="eğim 4")
ax.set_xlabel("adım uzunluğu h")
ax.set_ylabel("t = 4'teki hata")
ax.legend()
fig.savefig("hata-adim.png", dpi=150)Çıktı:
Euler doğrusunun eğimi: 1.067
RK4 doğrusunun eğimi : 3.990Kod grafiği hata-adim.png dosyasına kaydeder; aşağıdaki şekil aynı grafiğin aynı verilerle çizilmiş hâlidir:
RK4’ün eğimi \(3{,}99\) ile tam beklendiği gibi. Euler’in eğimi \(1\)’den biraz büyük çıktı, çünkü büyük adımlarda hatanın \(h^2\)’li terimi de hissediliyor; tablodaki oranların küçük adımlarda \(2\)’ye oturması bunu doğruluyor.
13.5 SciPy ile Çözüm: solve_ivp
Kendi yazdığımız fonksiyonlar yöntemleri anlamak için idealdir, ama gerçek işte iki eksikleri var: adım uzunluğunu bizim seçmemiz gerekir ve sonucun ne kadar doğru olduğunu söylemezler. SciPy’nin scipy.integrate.solve_ivp fonksiyonu ikisini de kendisi halleder.
Tanım 13.5 (Uyarlamalı Adım Kontrolü) Her adımda o adımın hatasını kestirip adım uzunluğunu bu kestirime göre ayarlayan yöntemlere uyarlamalı adımlı (adaptive step size) yöntemler denir. Kestirilen hata istenen hata payının altındaysa adım kabul edilir ve sonraki adım büyütülebilir; üstündeyse adım küçültülerek yeniden atılır.
Yani yöntem, çözümün yavaş değiştiği yerlerde uzun, hızlı değiştiği yerlerde kısa adımlar atar. solve_ivp’nin varsayılan yöntemi RK45, her adımda aynı \(f\) değerlerinden hem dördüncü hem beşinci mertebeden bir sonuç üretir; ikisinin farkı adımın hatasının kestirimidir. Hata payı rtol (bağıl) ve atol (mutlak) parametreleriyle verilir: her bileşende adım başına hata kabaca atol + rtol * abs(y) değerinin altında tutulur. Varsayılan değerler rtol=1e-3 ve atol=1e-6’dır.
- Sağ tarafı
fun(t, y)biçiminde yazın: ilk argüman zaman, ikincisi bilinmeyenlerin dizisidir. Fonksiyon \(y'\) değerlerini bir liste ya da dizi olarak döndürür. - Aralığı
t_span=(t0, T), başlangıç değeriniy0=[...]listesi olarak hazırlayın; tek bir denklemde biley0tek elemanlı bir listedir. sol = solve_ivp(fun, t_span, y0, rtol=..., atol=...)çağrısını yapın. Belirli anlardaki değerler gerekiyorsat_evalya dadense_output=Trueekleyin.sol.successdeğerini denetleyin. Sonuçsol.t(zamanlar) vesol.y(her satırı bir bilinmeyen) dizilerindedir.
Örnek 13.3 (SciPy ile Tanjant Fonksiyonu) \(y' = 1 + y^2\), \(y(0) = 0\) problemini \([0;\ 1{,}5]\) aralığında solve_ivp ile varsayılan ayarlarla çözünüz ve sonucu tam çözüm \(y = \tan t\) ile karşılaştırınız.
Çözüm
Adım 1 ve 2. Tam çözümün \(\tan t\) olduğunu Diferansiyel Denklemler notlarında görmüştük (bkz. Diferansiyel Denklemler); çözüm \(t \to \pi/2 \approx 1{,}5708\) iken sonsuza gider. Sağ taraf f(t, y) fonksiyonudur, t_span=(0, 1.5) ve y0=[0.0]’dır.
Adım 3. Hata paylarını vermeden çağırıyoruz:
import numpy as np
from scipy.integrate import solve_ivp
def f(t, y):
return 1 + y**2
sol = solve_ivp(f, (0, 1.5), [0.0])
print(sol.success, sol.message)
print("y'nin şekli:", sol.y.shape, " f çağrısı:", sol.nfev)
for tk, yk in zip(sol.t, sol.y[0]):
print(f"t = {tk:.6f} y = {yk:10.6f} tan t = {np.tan(tk):10.6f}")Çıktı:
True The solver successfully reached the end of the integration interval.
y'nin şekli: (1, 10) f çağrısı: 74
t = 0.000000 y = 0.000000 tan t = 0.000000
t = 0.000100 y = 0.000100 tan t = 0.000100
t = 0.001100 y = 0.001100 tan t = 0.001100
t = 0.011100 y = 0.011100 tan t = 0.011100
t = 0.111100 y = 0.111559 tan t = 0.111559
t = 0.822584 y = 1.077292 tan t = 1.077281
t = 1.083655 y = 1.887806 tan t = 1.887783
t = 1.344726 y = 4.348962 tan t = 4.347783
t = 1.447468 y = 8.071618 tan t = 8.067294
t = 1.500000 y = 14.114791 tan t = 14.101420Adım 4. sol.success değeri True’dur. sol.y’nin şekli \((1, 10)\)’dur: tek bilinmeyen ve \(10\) zaman noktası, yani \(9\) adım. İlk adımlar \(0{,}0001\), \(0{,}001\), \(0{,}01\) ve \(0{,}1\) uzunluğundadır; yöntem önce ölçeği yoklar, sonra adımı büyütür. Sonuç ise beklenenden kaba: \(t = 1{,}5\)’te bulunan \(14{,}114791\) değeri \(\tan 1{,}5 = 14{,}101420\)’den \(0{,}013\) kadar sapıyor. Bağıl hata \(0{,}013/14{,}1 \approx 10^{-3}\)’tür, yani varsayılan rtol=1e-3 kadardır.
\(\blacksquare\)
solve_ivp’nin varsayılan rtol=1e-3 değeri yalnız üç basamaklık bir doğruluk hedefler; sonucun kaç basamağına ihtiyacınız varsa rtol ve atol’u açıkça verin. Argüman sırasına da dikkat edin: solve_ivp sağ tarafı fun(t, y) biçiminde ister. Eski scipy.integrate.odeint fonksiyonu ise fun(y, t) sırasını kullanır; odeint için yazılmış bir fonksiyon solve_ivp’ye verilirse \(t\) ile \(y\) yer değiştirir ve sonuç sessizce yanlış çıkar.
Hata payını sıkılaştırınca ne olduğunu görelim. Aynı problemi varsayılan RK45 ve sekizinci mertebeden DOP853 yöntemleriyle, atol=1e-12 sabit tutup farklı rtol değerleriyle çözüyoruz:
import numpy as np
from scipy.integrate import solve_ivp
def f(t, y):
return 1 + y**2
print("yöntem rtol adım f çağrısı y(1.5) hatası")
for method in ["RK45", "DOP853"]:
for rtol in [1e-3, 1e-6, 1e-9, 1e-12]:
sol = solve_ivp(f, (0, 1.5), [0.0], method=method,
rtol=rtol, atol=1e-12)
err = abs(sol.y[0, -1] - np.tan(1.5))
print(f"{method:6s} {rtol:.0e} {sol.t.size - 1:4d}"
f" {sol.nfev:9d} {err:.2e}")Çıktı:
yöntem rtol adım f çağrısı y(1.5) hatası
RK45 1e-03 9 74 1.33e-02
RK45 1e-06 24 248 6.20e-05
RK45 1e-09 80 494 7.14e-08
RK45 1e-12 289 1748 5.75e-11
DOP853 1e-03 8 134 3.36e-03
DOP853 1e-06 12 218 2.55e-06
DOP853 1e-09 24 506 4.74e-09
DOP853 1e-12 46 566 3.03e-12rtol bin kat küçüldükçe hata da kabaca \(200\) ila \(1600\) kat küçülüyor; bedeli daha çok adım ve daha çok \(f\) çağrısıdır. Yüksek doğrulukta sekizinci mertebeden DOP853 çok daha ekonomiktir: rtol=1e-12 için RK45 \(1748\), DOP853 ise \(566\) çağrıyla ve daha küçük hatayla bitirir. RK45’in rtol=1e-6 ile attığı \(24\) adım şöyle dağılıyor:
solve_ivp'nin y′ = 1 + y², y(0) = 0 problemini rtol=1e-6 ile çözerken attığı 24 adım. Solda adım noktaları y = tan t eğrisi üzerinde, sağda her adımın uzunluğu (sütunun genişliği adımın kapladığı aralık, yüksekliği adım uzunluğu). İlk dört adım, yöntemin ölçeği yoklarken attığı 0,0001, 0,001, 0,01 ve 0,1 uzunluğundaki adımlardır; sonra adımlar en fazla 0,21'e çıkar ve çözüm dikleştikçe küçülür.Çözüm \(t = \pi/2\)’deki dikey asimptotuna yaklaştıkça dikleşiyor; yöntem de adımlarını kısaltarak buna uyuyor. Sabit adımlı bir yöntem aynı doğruluğa ulaşmak için çözümün en dik yerinde gereken küçük adımı baştan sona kullanmak zorunda kalırdı. Sıkı bir hata payıyla solve_ivp, tam çözümü bilinmeyen problemlerde bir referans değer üretmek için de kullanılır.
Örnek 13.4 (Euler Kestirimini Doğrulamak) \(y' = t - y^2\), \(y(0) = 0\) probleminde \(y(2)\)’yi solve_ivp ile sıkı bir hata payıyla hesaplayınız ve daha önce bulduğumuz Euler sonuçlarının (bkz. Örnek 13.1) gerçek hatalarını bulunuz.
Çözüm
rtol=1e-12 ve atol=1e-14 ile bulunan değeri doğru kabul edip Euler sonuçlarından çıkarıyoruz:
from scipy.integrate import solve_ivp
from odetools import euler
def f(t, y):
return t - y**2
ref = solve_ivp(f, (0, 2), [0.0], rtol=1e-12, atol=1e-14).y[0, -1]
print(f"solve_ivp: y(2) ≈ {ref:.10f}")
for n in [100, 200, 400]:
t, y = euler(f, 0.0, 0.0, 2 / n, n)
print(f"Euler, n = {n}: hata = {y[-1] - ref:.2e}")Çıktı:
solve_ivp: y(2) ≈ 1.1935759753
Euler, n = 100: hata = 1.55e-04
Euler, n = 200: hata = 7.48e-05
Euler, n = 400: hata = 3.67e-05Referans değer \(y(2) \approx 1{,}1935759753\)’tür. \(n = 400\) için gerçek hata \(3{,}67 \cdot 10^{-5}\)’tir; ardışık farklardan kestirdiğimiz \(3{,}8 \cdot 10^{-5}\) buna çok yakındır. Kestirilen hatayı çıkararak bulduğumuz \(1{,}193575\) değeri ise referanstan yalnızca yaklaşık \(10^{-6}\) kadar farklıdır. Hatalar her satırda yarıya iniyor: Euler yöntemi birinci mertebedendir.
\(\blacksquare\)
Çözümün adım noktaları dışındaki değerleri de çoğu zaman gerekir. t_eval=[...] parametresi sonucu istenen anlarda verir. dense_output=True ise sol.sol adında, her \(t\)’de çağrılabilen bir fonksiyon döndürür; bu fonksiyon her adımın içinde yöntemin kendi ara değer polinomunu kullanır. Bir de events parametresi var.
Tanım 13.6 (Olay Fonksiyonu) solve_ivp’ye events parametresiyle verilen ve g(t, y) biçiminde çağrılan fonksiyona olay fonksiyonu (event function) denir. Çözüm boyunca \(g\bigl(t, y(t)\bigr)\) işaret değiştirdiğinde solve_ivp sıfır yerini bulup sol.t_events listesine yazar. Fonksiyonun terminal niteliği True ise çözüm ilk olayda durdurulur; direction niteliği \(+1\) ya da \(-1\) ise yalnız artarak ya da azalarak sıfırdan geçişler sayılır.
Yani olay fonksiyonu “çözüm şu koşulu ne zaman sağlıyor?” sorusunu yanıtlar: bir topun yere düştüğü an, bir sarkacın en alttan geçtiği an ya da bir nüfusun belli bir sayıya ulaştığı an. Sıfır yeri, adımın içindeki ara değer polinomu üzerinde Kök Bulma ve Optimizasyon bölümündekilere benzer bir kök bulma yöntemiyle bulunur. events parametresine birden çok olay fonksiyonu da verilebildiği için sol.t_events, her olay fonksiyonu için bir dizi içeren bir listedir. Tek olay fonksiyonunda olay anları sol.t_events[0] dizisindedir ve ilk olay anı sol.t_events[0][0]’dır. Olay anlarındaki çözüm vektörleri de aynı düzende sol.y_events listesinde durur: sol.y_events[0] dizisinin her satırı, bir olay anındaki çözüm vektörüdür. Tanjant probleminde çözümün \(10\)’a ulaştığı anı bulalım ve ara noktalardaki değerleri yoğun çıktıdan okuyalım:
import numpy as np
from scipy.integrate import solve_ivp
def f(t, y):
return 1 + y**2
def reach10(t, y):
return y[0] - 10 # y = 10 olduğunda sıfır
reach10.terminal = True # olay olunca çözümü durdur
sol = solve_ivp(f, (0, 1.5), [0.0], rtol=1e-10, atol=1e-12,
dense_output=True, events=reach10)
print(f"olay anı : {sol.t_events[0][0]:.12f}")
print(f"arctan 10: {np.arctan(10):.12f}")
print(f"son nokta: t = {sol.t[-1]:.12f}, y = {sol.y[0, -1]:.12f}")
t = np.array([0.25, 0.5, 0.75, 1.0, 1.25])
err = sol.sol(t)[0] - np.tan(t)
print(f"ara noktalardaki en büyük hata: {np.max(np.abs(err)):.1e}")Çıktı:
olay anı : 1.471127674330
arctan 10: 1.471127674304
son nokta: t = 1.471127674330, y = 10.000000000000
ara noktalardaki en büyük hata: 2.9e-10\(\tan t = 10\) denkleminin çözümü \(\arctan 10\)’dur ve olay anı onunla on ondalık basamakta uyuşuyor. terminal=True olduğundan çözüm \(t = 1{,}5\)’e gitmeden olay anında durdu. Ara noktalardaki en büyük hata \(2{,}9 \cdot 10^{-10}\)’dur: yoğun çıktı, istenen hata payına uygun doğruluktadır.
Örnek 13.5 (Sonlu Zamanda Patlayan Bir Çözüm) \(y' = y^2\), \(y(0) = 1\) problemini \([0, 2]\) aralığında solve_ivp ile çözmeye çalışınız ve ne olduğunu tam çözümle açıklayınız.
Çözüm
Tam çözüm. Denklem ayrılabilir: \(y^{-2}\,dy = dt\) integrallenince \(-1/y = t + c\) bulunur ve \(y(0) = 1\) koşulu \(c = -1\) verir. Çözüm
\[y(t) = \frac{1}{1 - t}\]
dir ve \(t \to 1^-\) iken sonsuza gider; \(t = 1\)’in ötesine uzatılamaz.
Sayısal çözüm.
from scipy.integrate import solve_ivp
def f(t, y):
return y**2
sol = solve_ivp(f, (0, 2), [1.0])
print(sol.status, sol.success)
print(sol.message)
print("son t:", sol.t[-1], " son y:", sol.y[0, -1])Çıktı:
-1 False
Required step size is less than spacing between numbers.
son t: 0.9999286400563746 son y: 480848617754195.6status değeri \(-1\) ve success değeri False’tur: çözüm başarısız oldu. Yöntem \(t = 1\)’e yaklaştıkça adımı küçülttü ve gereken adım, \(t\)’nin yakınındaki iki ardışık kayan noktalı sayı arasındaki aralıktan bile küçük olunca durdu. Son \(y\) değeri güvenilir değildir; o anda \(\frac{1}{1 - t} \approx 14\,000\) olmalıydı. Bu çıktıdan öğrendiğimiz asıl şey, çözümün \(t = 1\) civarında patladığıdır.
\(\blacksquare\)
solve_ivp çözemediği bir problemde istisna fırlatmaz; o ana kadar hesapladığı değerleri döndürür ve durumu sol.status, sol.success ve sol.message alanlarına yazar. Sonucu kullanmadan önce bu alanları mutlaka denetleyin.
13.6 Katı Denklemler
Bazı denklemlerde adım uzunluğunu doğruluk değil kararlılık belirler. \(\lambda > 0\) büyük bir sayı olmak üzere
\[y' = -\lambda\,(y - \cos t), \qquad y(0) = 0\]
problemini düşünelim. Bu birinci mertebeden doğrusal denklemin çözümü, integrasyon çarpanıyla (bkz. Diferansiyel Denklemler)
\[y(t) = \frac{\lambda^2}{\lambda^2 + 1}\left(\cos t + \frac{\sin t}{\lambda} - e^{-\lambda t}\right)\]
bulunur. \(e^{-\lambda t}\) terimi çok hızlı söner ve çözüm kısa sürede \(\cos t\)’ye yapışır; sonrasında çözüm yavaş ve düzgündür. Euler yöntemini \(\lambda = 50\) için \([0, 3]\) aralığında farklı adımlarla çalıştıralım:
import numpy as np
from odetools import euler
lam = 50.0
def f(t, y):
return -lam * (y - np.cos(t))
def exact(t):
a = lam**2 / (lam**2 + 1)
return a * (np.cos(t) + np.sin(t) / lam - np.exp(-lam * t))
for h in [0.01, 0.03, 0.0375, 0.04, 0.05]:
n = round(3 / h)
t, y = euler(f, 0.0, 0.0, h, n)
err = abs(y[-1] - exact(t[-1]))
print(f"h = {h:6.4f} h·λ = {h * lam:5.2f} y(3) hatası = {err:.3e}")Çıktı:
h = 0.0100 h·λ = 0.50 y(3) hatası = 9.839e-05
h = 0.0300 h·λ = 1.50 y(3) hatası = 2.956e-04
h = 0.0375 h·λ = 1.88 y(3) hatası = 3.926e-04
h = 0.0400 h·λ = 2.00 y(3) hatası = 9.996e-01
h = 0.0500 h·λ = 2.50 y(3) hatası = 3.677e+10\(h\lambda < 2\) olduğu sürece hata küçüktür ve kabaca \(h\) ile orantılıdır. \(h\lambda = 2\)’de hata \(1\) civarında takılır, \(h\lambda > 2\)’de ise sayılar patlar. Oysa tam çözüm \(|y| \le 1\) aralığında kalan zararsız bir fonksiyondur. Nedeni şu önermedir:
Önerme 13.1 (Euler Yönteminin Kararlılık Koşulu) \(\lambda > 0\) olsun. \(y' = -\lambda y\), \(y(0) = y_0 \ne 0\) problemine \(h\) adımlı Euler yöntemi uygulandığında \(y_i = (1 - h\lambda)^i\,y_0\) olur. Bu dizi ancak ve ancak \(0 < h < 2/\lambda\) ise sıfıra yakınsar; \(h > 2/\lambda\) ise \(|y_i| \to \infty\) olur.
İspat
Euler adımı \(y_{i+1} = y_i - h\lambda\,y_i\), yani \(y_{i+1} = (1 - h\lambda)\,y_i\)’dir; tümevarımla \(y_i = (1 - h\lambda)^i\,y_0\) bulunur. Bir \(q^i\) geometrik dizisi ancak ve ancak \(|q| < 1\) ise sıfıra gider ve \(|q| > 1\) ise mutlak değerce sonsuza gider. \(q = 1 - h\lambda\) için
\[|q| < 1 \iff -1 < 1 - h\lambda < 1 \iff 0 < h\lambda < 2\]
olur. \(h > 2/\lambda\) ise \(q < -1\), yani \(|q| > 1\)’dir.
\(\blacksquare\)
Tam çözüm \(y_0 e^{-\lambda t}\) her \(\lambda > 0\) için sönerken, Euler dizisi ancak \(h < 2/\lambda\) ise söner. Yukarıdaki denklem de doğrusal olduğundan, iki Euler dizisinin farkı her adımda tam olarak \(1 - h\lambda\) ile çarpılır. \(h\lambda > 2\) iken bu çarpanın mutlak değeri \(1\)’den büyüktür ve her adımda yapılan küçük hata katlanarak büyür; \(\lambda = 50\) için sınır \(h = 0{,}04\)’tür. Euler, RK4 ve RK45 gibi, yeni \(y_{i+1}\) değerini bilinen değerlerden doğrudan hesaplayan yöntemlere açık (explicit) yöntem denir.
Tanım 13.7 (Katı Denklem) Çözümünde çok hızlı sönen bir bileşenle yavaş değişen bir bileşenin bir arada bulunduğu ve bu yüzden açık yöntemlerin adım uzunluğunu doğruluğun değil kararlılığın sınırladığı denklemlere katı (stiff) denklem denir.
Yani katı bir denklemde hızlı bileşen kısa sürede söner ve sonrasında çözüm yavaş değişir; yine de açık yöntemler, sönmüş bileşen kararsızlaşmasın diye baştan sona küçük adım atmak zorundadır. Bu tür denklemler için kapalı (implicit) yöntemler geliştirilmiştir: her adımda yeni değeri bir denklemin çözümü olarak bulurlar ve büyük adımlarda da kararlı kalırlar. solve_ivp’de bunlar Radau ve BDF yöntemleridir; LSODA ise denklemin katı olup olmadığını kendisi sezip yöntem değiştirir. \(\lambda = 1000\) için dört yöntemi karşılaştıralım:
import numpy as np
from scipy.integrate import solve_ivp
lam = 1000.0
def f(t, y):
return -lam * (y - np.cos(t))
def exact(t):
a = lam**2 / (lam**2 + 1)
return a * (np.cos(t) + np.sin(t) / lam - np.exp(-lam * t))
print("yöntem adım f çağrısı y(10) hatası")
for method in ["RK45", "Radau", "BDF", "LSODA"]:
sol = solve_ivp(f, (0, 10), [0.0], method=method,
rtol=1e-6, atol=1e-9)
err = abs(sol.y[0, -1] - exact(10))
print(f"{method:6s} {sol.t.size - 1:5d} {sol.nfev:9d} {err:.1e}")Çıktı:
yöntem adım f çağrısı y(10) hatası
RK45 3587 21548 4.5e-07
Radau 81 624 2.1e-08
BDF 209 478 2.1e-09
LSODA 272 542 5.1e-09Çözüm ilk birkaç binde birlik süreden sonra \(\cos t\) kadar yavaş değişmesine karşın RK45 \(3587\) adım atıyor; ortalama adım \(10/3587 \approx 0{,}0028\)’dir. Adımı doğruluk değil kararlılık sınırlıyor. Radau \(81\), BDF \(209\) adımla ve daha küçük hatayla bitiriyor. Kapalı yöntemler her adımda \(f\)’nin Jacobi matrisine (bkz. Kök Bulma ve Optimizasyon) ihtiyaç duyar; buradaki gibi tek bir denklemde bu matris yalnızca \(\partial f/\partial y\) türevidir. jac parametresiyle verilmezse yöntemler onu sonlu farklarla kendileri hesaplar.
13.7 Denklem Sistemleri ve Yüksek Mertebeden Denklemler
Gerçek modellerin çoğunda birden çok bilinmeyen birbirini etkiler; bu yüzden tek bir denklem yerine bir denklem sistemi çözeriz.
Tanım 13.8 (Birinci Mertebeden Denklem Sistemi) \(\mathbf{y}(t) = \bigl(y_1(t), \dots, y_m(t)\bigr)\) bilinmeyen bir vektör fonksiyonu ve \(\mathbf{F}\colon \mathbb{R} \times \mathbb{R}^m \to \mathbb{R}^m\) verilmiş bir fonksiyon olsun.
\[\mathbf{y}' = \mathbf{F}(t, \mathbf{y}), \qquad \mathbf{y}(t_0) = \mathbf{y}_0\]
problemine birinci mertebeden bir denklem sistemi için başlangıç değer problemi denir.
Yani tek denklemdeki her şey aynen geçerlidir; yalnız \(y\) bir sayı değil bir vektör, \(f\) de vektör değerli bir fonksiyondur. Euler ve RK4 formülleri hiç değişmeden vektörler için de anlamlıdır ve solve_ivp zaten vektörlerle çalışır. Yüksek mertebeden bir denklem de böyle bir sisteme çevrilir.
- \(y^{(m)} = g\bigl(t, y, y', \dots, y^{(m-1)}\bigr)\) denkleminde yeni bilinmeyenler tanımlayın: \(u_1 = y\), \(u_2 = y'\), …, \(u_m = y^{(m-1)}\).
- Sistemi yazın: \(u_1' = u_2\), …, \(u_{m-1}' = u_m\) ve \(u_m' = g(t, u_1, \dots, u_m)\).
- Başlangıç koşullarını bir vektörde toplayın: \(\mathbf{u}(t_0)\)’ın bileşenleri \(y(t_0), y'(t_0), \dots, y^{(m-1)}(t_0)\)’dır.
Örnek 13.6 (Harmonik Salınıcı) \(y'' + y = 0\), \(y(0) = 1\), \(y'(0) = 0\) problemini birinci mertebeden bir sisteme çevirip solve_ivp ile \([0, 2\pi]\) aralığında çözünüz; sonucu tam çözüm \(y = \cos t\) ile karşılaştırınız.
Çözüm
Adım 1. \(u_1 = y\) ve \(u_2 = y'\) diyelim.
Adım 2. \(u_1' = y' = u_2\) ve \(u_2' = y'' = -y = -u_1\) olur; yani \(\mathbf{u}' = (u_2,\ -u_1)\).
Adım 3. \(\mathbf{u}(0) = (1, 0)\)’dır.
Kodda vektörü u, bileşenlerini okunaklı olsun diye y, v adlarıyla açıyoruz. Sonucu çeyrek periyotlarda istiyoruz:
import numpy as np
from scipy.integrate import solve_ivp
def oscillator(t, u):
y, v = u # u = (y, y')
return [v, -y]
t = np.pi / 2 * np.arange(5) # 0, π/2, π, 3π/2, 2π
sol = solve_ivp(oscillator, (0, 2 * np.pi), [1.0, 0.0], t_eval=t,
rtol=1e-10, atol=1e-12)
print("sol.y şekli:", sol.y.shape)
for tk, yk, vk in zip(sol.t, sol.y[0], sol.y[1]):
print(f"t = {tk:.4f} y = {yk:+.10f} y' = {vk:+.10f}")
print("cos ile en büyük fark:", np.max(np.abs(sol.y[0] - np.cos(t))))Çıktı:
sol.y şekli: (2, 5)
t = 0.0000 y = +1.0000000000 y' = +0.0000000000
t = 1.5708 y = -0.0000000000 y' = -1.0000000000
t = 3.1416 y = -1.0000000000 y' = +0.0000000000
t = 4.7124 y = +0.0000000000 y' = +0.9999999999
t = 6.2832 y = +0.9999999999 y' = -0.0000000000
cos ile en büyük fark: 7.59422524865272e-11sol.y artık iki satırlıdır: ilk satır \(y\), ikinci satır \(y'\) değerleridir. Değerler \(\cos t\) ve \(-\sin t\) ile \(10^{-10}\) duyarlıkla uyuşuyor; çıktıdaki -0.0000000000 gibi değerler, sıfıra çok yakın negatif sayılardır.
\(\blacksquare\)
Kendi yazdığımız fonksiyonlar da sistemlerde çalışır; tek koşul, sağ tarafın bir NumPy dizisi döndürmesidir, çünkü bir liste h / 2 gibi bir kayan noktalı sayıyla (float) çarpılamaz. Euler ve RK4’ü harmonik salınıcıda \(h = 2\pi/40\) ile bir tam tur çalıştıralım ve faz düzleminde \((y, y')\) noktasının orijine uzaklığını, \(r = \sqrt{y^2 + y'^2}\)’yi izleyelim. Tam çözümde \(r = \sqrt{\cos^2 t + \sin^2 t} = 1\) sabittir:
import numpy as np
from odetools import euler, rk4
def oscillator(t, u):
return np.array([u[1], -u[0]]) # NumPy dizisi döndürmeli
n = 40
h = 2 * np.pi / n # bir tam tur: t = 2π
for method in (euler, rk4):
t, u = method(oscillator, 0.0, np.array([1.0, 0.0]), h, n)
r = np.hypot(u[:, 0], u[:, 1]) # faz düzleminde orijine uzaklık
print(f"{method.__name__:5s} u şekli: {u.shape} "
f"y(2π) = {u[-1, 0]:.6f} r(2π) = {r[-1]:.6f}")
print("(1 + h²)^(n/2) =", (1 + h**2) ** (n / 2))Çıktı:
euler u şekli: (41, 2) y(2π) = 1.626114 r(2π) = 1.628225
rk4 u şekli: (41, 2) y(2π) = 0.999996 r(2π) = 0.999996
(1 + h²)^(n/2) = 1.6282250240481193RK4 tur sonunda \(r = 0{,}999996\) ile neredeyse başladığı yere döner. Euler’in yarıçapı ise \(1{,}628\)’e çıkmıştır ve bu sayı tam olarak \((1 + h^2)^{n/2}\)’dir. Nedeni şudur: \(v = y'\) yazarsak Euler adımı
\[\begin{pmatrix} y_{i+1} \\ v_{i+1} \end{pmatrix} = \begin{pmatrix} 1 & h \\ -h & 1 \end{pmatrix} \begin{pmatrix} y_i \\ v_i \end{pmatrix}\]
matris çarpımıdır. Bu matris, \(\sqrt{1 + h^2}\) sayısıyla bir dönme matrisinin çarpımıdır; her adım vektörü biraz döndürür ve boyunu \(\sqrt{1 + h^2}\) kat uzatır. \(n\) adımda boy \((1 + h^2)^{n/2}\) katına çıkar:
Euler yöntemi salınımlı sistemlerde enerjiyi yapay olarak büyütür; adım küçüldükçe bu etki azalır, ama her \(h > 0\) için vardır ve uzun sürelerde birikir. RK4 gibi yüksek mertebeden açık yöntemlerde de enerji sürüklenir, yalnızca çok daha yavaş: yukarıda RK4’ün yarıçapı bir turda \(0{,}999996\)’ya indi. Uzun süreli salınım hesaplarında bu yüzden ya hata payı sıkı tutulan yüksek mertebeden yöntemler ya da enerji hatası zamanla birikmeyecek biçimde kurulmuş simplektik (symplectic) yöntemler kullanılır.
13.8 Faz Düzlemi
İki bilinmeyenli bir sistemin çözümünü, zamanı gizleyip \((y_1, y_2)\) noktasının düzlemde çizdiği eğri olarak da düşünebiliriz. Bu bakış özellikle sağ tarafı zamana bağlı olmayan sistemlerde verimlidir.
Tanım 13.9 (Otonom Sistem) Sağ tarafı \(t\)’ye açıkça bağlı olmayan \(\mathbf{y}' = \mathbf{F}(\mathbf{y})\) sistemine otonom (autonomous) sistem denir.
Yani otonom bir sistemde hareket kuralı zamanla değişmez: aynı noktadan bugün de yarın da aynı biçimde çıkılır. Harmonik salınıcı, aşağıdaki Lotka–Volterra modeli ve sarkaç otonomdur.
Tanım 13.10 (Faz Düzlemi ve Yörünge) İki bilinmeyenli otonom bir \(x' = F(x, y)\), \(y' = G(x, y)\) sisteminin bir çözümü, \(t\) değiştikçe \(xy\)-düzleminde bir eğri çizer. Bu eğriye çözümün yörüngesi (trajectory), yörüngelerin çizildiği düzleme faz düzlemi (phase plane) denir. Tipik yörüngeleri bir arada gösteren resme de faz portresi (phase portrait) denir.
Yani faz düzlemi zamanı eksenlerden çıkarıp yalnız durumu gösterir. Harmonik salınıcının yörüngesi bir çemberdir: çemberin her noktası bir an, bir tam tur bir periyottur. \(F\) ve \(G\)’nin kısmi türevleri sürekli olduğunda varlık ve teklik teoremi sayesinde iki farklı yörünge kesişemez; kesişselerdi kesişim noktasından iki farklı çözüm geçerdi.
Tanım 13.11 (Denge Noktası) \(\mathbf{F}(\mathbf{y}^*) = \mathbf{0}\) eşitliğini sağlayan \(\mathbf{y}^*\) noktasına otonom \(\mathbf{y}' = \mathbf{F}(\mathbf{y})\) sisteminin denge noktası (equilibrium point) denir.
Yani denge noktasından başlayan çözüm hiç kıpırdamaz: \(\mathbf{y}(t) = \mathbf{y}^*\). Bir denge noktasının yakınındaki davranış hakkında \(\mathbf{F}\)’nin o noktadaki Jacobi matrisi \(J\) çok şey söyler (Jacobi matrisini Kök Bulma ve Optimizasyon bölümünde Newton yönteminde kullanmıştık). \(\mathbf{y} = \mathbf{y}^* + \mathbf{u}\) yazıp \(\mathbf{F}\)’yi doğrusallaştırınca küçük \(\mathbf{u}\) için \(\mathbf{u}' \approx J\mathbf{u}\) olur ve bu doğrusal sistemin davranışını \(J\)’nin özdeğerleri belirler (bkz. Lineer Cebir). Özdeğerler \(\pm i\omega\) biçiminde saf sanal sayılarsa doğrusal sistemin çözümleri \(2\pi/\omega\) periyotla döner. Asıl sistemin denge noktası çevresindeki yörüngeleri kapalıysa, ki aşağıdaki iki örnekte öyledir, noktaya çok yakın yörüngelerin periyodu da \(2\pi/\omega\)’ya yaklaşır. Harmonik salınıcıda \(J = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}\)’dir, özdeğerler \(\pm i\) ve periyot \(2\pi\)’dir.
13.9 Lotka–Volterra Av-Avcı Modeli
Faz düzleminin en ünlü örneklerinden biri, bir av türüyle onu yiyen bir avcı türün nüfuslarını birlikte modelleyen Lotka–Volterra denklemleridir.
Tanım 13.12 (Lotka–Volterra Modeli) \(a, b, c, d\) pozitif sabitler olmak üzere
\[x' = ax - bxy, \qquad y' = -cy + dxy\]
sistemine Lotka–Volterra modeli denir. Burada \(x(t) > 0\) av, \(y(t) > 0\) avcı nüfusudur.
Yani avcı yokken av üstel olarak çoğalır (\(ax\) terimi), av yokken avcı üstel olarak azalır (\(-cy\) terimi); iki türün karşılaşmaları (\(xy\) terimleri) avı azaltıp avcıyı çoğaltır. \(x' = x(a - by)\) ve \(y' = y(dx - c)\) yazılışından, sistemin \(x, y > 0\) bölgesindeki tek denge noktasının \(\left(\frac{c}{d}, \frac{a}{b}\right)\) olduğu görülür. Bu bölümde \(a = 1\), \(b = 0{,}5\), \(c = 0{,}75\), \(d = 0{,}25\) alacağız; denge noktası \((3, 2)\)’dir. Önce denge noktasındaki Jacobi matrisini
\[J = \begin{pmatrix} a - by^* & -bx^* \\ dy^* & -c + dx^* \end{pmatrix} = \begin{pmatrix} 0 & -1{,}5 \\ 0{,}5 & 0 \end{pmatrix}\]
yazıp özdeğerlerini NumPy ile Lineer Cebir bölümündeki gibi hesaplayalım:
import numpy as np
a, b, c, d = 1.0, 0.5, 0.75, 0.25
x_eq, y_eq = c / d, a / b # denge noktası (3, 2)
J = np.array([[a - b * y_eq, -b * x_eq],
[d * y_eq, -c + d * x_eq]]) # dengedeki Jacobi matrisi
print("denge noktası:", (x_eq, y_eq))
print(J)
print("özdeğerler :", np.linalg.eigvals(J))
print("2π / √(ac) :", 2 * np.pi / np.sqrt(a * c))Çıktı:
denge noktası: (3.0, 2.0)
[[ 0. -1.5]
[ 0.5 0. ]]
özdeğerler : [0.+0.8660254j 0.-0.8660254j]
2π / √(ac) : 7.255197456936871Özdeğerler \(\pm i\sqrt{ac} \approx \pm 0{,}866\,i\)’dir; çünkü \(J\)’nin izi sıfır, determinantı \(ac = 0{,}75\)’tir. Denge noktası çevresindeki küçük salınımların periyodu bu yüzden \(2\pi/\sqrt{ac} \approx 7{,}2552\) olmalıdır. Yörüngelerin gerçekten kapalı olduğunu görmek için bir korunan büyüklük kullanacağız.
Tanım 13.13 (Korunan Büyüklük) Bir \(V(\mathbf{y})\) fonksiyonu sistemin her \(\mathbf{y}(t)\) çözümü boyunca sabit kalıyorsa, yani \(\frac{d}{dt}V\bigl(\mathbf{y}(t)\bigr) = 0\) ise, \(V\)’ye sistemin korunan büyüklüğü (conserved quantity) ya da ilk integrali denir.
Yani her yörünge \(V\)’nin bir seviye eğrisinin üzerinde kalır. Fizikteki enerjinin korunumu bunun en bilinen örneğidir. Korunan büyüklük sayısal çözüm için de değerli bir denetimdir: \(V\) sayısal çözüm boyunca ne kadar değişiyorsa çözüm en az o ölçüde hatalıdır.
Önerme 13.2 (Lotka–Volterra Modelinin Korunan Büyüklüğü) \(x, y > 0\) bölgesinde
\[V(x, y) = dx - c\ln x + by - a\ln y\]
fonksiyonu Lotka–Volterra modelinin korunan büyüklüğüdür.
İspat
Zincir kuralıyla
\[\frac{d}{dt}V\bigl(x(t), y(t)\bigr) = \left(d - \frac{c}{x}\right)x' + \left(b - \frac{a}{y}\right)y'\]
olur. \(x' = x(a - by)\) ve \(y' = y(dx - c)\) yerine konursa
\[\frac{d}{dt}V = (dx - c)(a - by) + (by - a)(dx - c) = 0\]
bulunur.
\(\blacksquare\)
\(V\), iki konveks tek değişkenli fonksiyonun toplamıdır: \(dx - c\ln x\) ifadesi en küçük değerini \(x = c/d\)’de, \(by - a\ln y\) ifadesi \(y = a/b\)’de alır ve ikisi de \(0\)’a ya da sonsuza giderken sonsuza gider. Bu yüzden \(V\)’nin seviye eğrileri \((3, 2)\) noktasını çevreleyen kapalı eğrilerdir; her yörünge kapalıdır ve her çözüm periyodiktir.
Sistemi solve_ivp ile \((6, 2)\) başlangıcından çözelim. Parametreleri fonksiyonun içine gömmek yerine args parametresiyle veriyoruz: solve_ivp, args=(a, b, c, d) demetini hem sağ tarafa hem olay fonksiyonuna ek argüman olarak geçirir. Avın denge değeri \(3\)’ü yukarıdan aşağıya kestiği anları bir olay fonksiyonuyla yakalayıp periyodu ölçüyor, \(V\)’nin çözüm boyunca ne kadar değiştiğine de bakıyoruz:
import numpy as np
from scipy.integrate import solve_ivp
def lotka(t, z, a, b, c, d):
x, y = z
return [a * x - b * x * y, -c * y + d * x * y]
def V(x, y, a, b, c, d):
return d * x - c * np.log(x) + b * y - a * np.log(y)
def prey_down(t, z, a, b, c, d):
return z[0] - c / d # av, denge değerini yukarıdan aşağı keser
prey_down.direction = -1
p = (1.0, 0.5, 0.75, 0.25) # a, b, c, d
sol = solve_ivp(lotka, (0, 30), [6.0, 2.0], args=p, events=prey_down,
rtol=1e-10, atol=1e-12, dense_output=True)
for tk in [0, 2, 4, 6, 8]:
x, y = sol.sol(tk)
print(f"t = {tk}: av = {x:6.3f} avcı = {y:6.3f}")
v = V(sol.y[0], sol.y[1], *p)
print("V'nin en büyük sapması:", np.max(np.abs(v - v[0])))
print("olay anları:", np.round(sol.t_events[0], 6))
print("periyot :", np.round(np.diff(sol.t_events[0]), 6))Çıktı:
t = 0: av = 6.000 avcı = 2.000
t = 2: av = 1.817 avcı = 3.303
t = 4: av = 1.361 avcı = 1.447
t = 6: av = 3.402 avcı = 0.942
t = 8: av = 5.591 avcı = 2.695
V'nin en büyük sapması: 3.400024706223803e-11
olay anları: [ 1.350724 8.934098 16.517472 24.100846]
periyot : [7.583374 7.583374 7.583374]\(V\), \(30\) birimlik süre boyunca \(3{,}4 \cdot 10^{-11}\)’den fazla değişmiyor; çözüm hata payı ölçüsünde doğru. Olay anları arasındaki süre, yani periyot, \(7{,}583374\)’tür. Tablo modelin döngüsünü gösteriyor: av çokken avcı çoğalır (\(t = 2\)), avcı çoğalınca av azalır, av azalınca avcı da azalır (\(t = 4\)) ve avcı azalınca av yeniden çoğalır (\(t = 6\) ve \(t = 8\)).
Periyodun genliğe nasıl bağlı olduğunu görmek için birkaç yörüngeyi birlikte çizelim. Kod bir önceki kodun devamıdır; sol, lotka, prey_down ve p oradan gelir:
import matplotlib.pyplot as plt
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
t = np.linspace(0, 30, 601)
x, y = sol.sol(t)
ax1.plot(t, x, label="av x(t)")
ax1.plot(t, y, label="avcı y(t)")
ax1.set_xlabel("t")
ax1.legend()
for x0 in [3.5, 4.5, 6.0, 8.0]:
s = solve_ivp(lotka, (0, 30), [x0, 2.0], args=p, events=prey_down,
rtol=1e-10, atol=1e-12, dense_output=True)
T = s.t_events[0][1] - s.t_events[0][0]
tt = np.linspace(0, T, 400)
ax2.plot(*s.sol(tt))
print(f"x(0) = {x0}: periyot = {T:.4f}")
ax2.plot(3, 2, "ko")
ax2.set_xlabel("av x")
ax2.set_ylabel("avcı y")
fig.savefig("lotka-volterra.png", dpi=150)Çıktı:
x(0) = 3.5: periyot = 7.2684
x(0) = 4.5: periyot = 7.3556
x(0) = 6.0: periyot = 7.5834
x(0) = 8.0: periyot = 7.9976Kod iki grafiği lotka-volterra.png dosyasına kaydeder; aşağıdaki şekil aynı grafiklerin aynı verilerle çizilmiş hâlidir:
Yörünge denge noktasına yaklaştıkça periyot küçülüyor ve doğrusallaştırmanın verdiği \(7{,}2552\) değerine yaklaşıyor: \(x(0) = 3{,}5\) için periyot \(7{,}2684\)’tür.
13.10 Sarkaç
Bölümün son örneği, uzunluğu \(L\) olan ağırlıksız bir çubuğun ucundaki noktasal bir kütleden oluşan basit sarkaçtır. Sarkacın düşeyle yaptığı açı \(\theta(t)\) ise Newton’un ikinci yasası
\[\theta'' = -\frac{g}{L}\sin\theta\]
denklemini verir; \(g\) yerçekimi ivmesidir. Zaman birimini \(g/L = 1\) olacak biçimde seçersek denklem \(\theta'' = -\sin\theta\) olur. Açısal hızı \(\omega = \theta'\) ile gösterip denklemi sisteme çevirirsek
\[\theta' = \omega, \qquad \omega' = -\sin\theta\]
otonom sistemi elde edilir. Küçük açılarda \(\sin\theta \approx \theta\) olduğundan denklem harmonik salınıcıya, \(\theta'' = -\theta\) denklemine dönüşür ve periyot \(2\pi\) olur. Büyük açılarda ise durum değişir. Sistemin korunan büyüklüğü enerjidir:
\[E(\theta, \omega) = \frac{\omega^2}{2} - \cos\theta.\]
Gerçekten
\[\frac{d}{dt}E = \omega\,\omega' + \sin\theta\,\theta' = -\omega\sin\theta + \omega\sin\theta = 0\]
olur. \(\frac{\omega^2}{2}\) terimi kinetik, \(-\cos\theta\) terimi potansiyel enerjiye karşılık gelir. Enerjinin korunumu, periyot için kapalı bir formül verir.
Teorem 13.2 (Sarkacın Periyodu) \(0 < \theta_0 < \pi\) olsun. \(\theta'' = -\sin\theta\), \(\theta(0) = \theta_0\), \(\theta'(0) = 0\) sarkacının periyodu
\[T(\theta_0) = 4K(k), \qquad k = \sin\frac{\theta_0}{2}\]
dir. Burada
\[K(k) = \int_0^{\pi/2} \frac{d\varphi}{\sqrt{1 - k^2\sin^2\varphi}}\]
birinci türden tam eliptik integraldir.
İspat
Enerjinin korunumundan \(\frac{\omega^2}{2} - \cos\theta = -\cos\theta_0\), yani \(\omega^2 = 2(\cos\theta - \cos\theta_0)\)’dır. Sarkaç \(\theta_0\)’dan bırakılınca \(\omega < 0\) olur ve sarkaç \(\theta = 0\)’a varana kadar böyle kalır. Denklem \(\theta \mapsto -\theta\) ve \(t \mapsto -t\) dönüşümleri altında değişmediğinden hareketin dört çeyreği eşit sürer; \(\theta_0\)’dan \(0\)’a iniş bunlardan biridir. \(dt = d\theta/|\omega|\) olduğundan
\[\frac{T}{4} = \int_0^{\theta_0} \frac{d\theta}{\sqrt{2(\cos\theta - \cos\theta_0)}}\]
bulunur. \(\cos\theta = 1 - 2\sin^2\frac{\theta}{2}\) özdeşliğiyle \(2(\cos\theta - \cos\theta_0) = 4\left(k^2 - \sin^2\frac{\theta}{2}\right)\) olur. \(\sin\frac{\theta}{2} = k\sin\varphi\) değişken değiştirmesini yapalım: \(\theta\), \(0\)’dan \(\theta_0\)’a giderken \(\varphi\), \(0\)’dan \(\pi/2\)’ye gider. Ayrıca
\[\frac{1}{2}\cos\frac{\theta}{2}\,d\theta = k\cos\varphi\,d\varphi, \qquad \cos\frac{\theta}{2} = \sqrt{1 - k^2\sin^2\varphi}\]
ve \(\sqrt{4(k^2 - k^2\sin^2\varphi)} = 2k\cos\varphi\)’dir. Bunlar yerine konunca
\[ \begin{aligned} \frac{T}{4} &= \int_0^{\pi/2} \frac{1}{2k\cos\varphi} \cdot \frac{2k\cos\varphi\,d\varphi}{\sqrt{1 - k^2\sin^2\varphi}}\\[1mm] &= \int_0^{\pi/2} \frac{d\varphi}{\sqrt{1 - k^2\sin^2\varphi}} = K(k) \end{aligned} \]
elde edilir. İlk integralin \(\theta = \theta_0\)’daki tekilliği değişken değiştirmeyle ortadan kalktı.
\(\blacksquare\)
Yani \(\theta_0 \to 0\) iken \(k \to 0\), \(K(k) \to \pi/2\) ve \(T \to 2\pi\) olur: harmonik salınıcının periyodu. \(\theta_0 \to \pi\) iken ise \(k \to 1\) ve \(K(k) \to \infty\)’dur: tepeye yakın bırakılan sarkacın periyodu sınırsız büyür. Teoremi sayısal olarak sınayalım. Periyodu, \(\theta\)’nın aşağıdan yukarıya sıfırdan geçtiği ardışık iki anın farkı olarak ölçüyoruz; \(K\) için SciPy’nin scipy.special.ellipk fonksiyonunu kullanıyoruz:
import numpy as np
from scipy.integrate import solve_ivp
from scipy.special import ellipk
def pendulum(t, z):
theta, omega = z
return [omega, -np.sin(theta)]
def upward(t, z):
return z[0] # θ = 0 geçişleri
upward.direction = 1 # yalnız aşağıdan yukarıya geçişler
print(" θ0 T (sayısal) 4K(sin²(θ0/2)) T / 2π")
for theta0 in [0.1, 1.0, 2.0, 3.0]:
sol = solve_ivp(pendulum, (0, 60), [theta0, 0.0], events=upward,
rtol=1e-10, atol=1e-12)
te = sol.t_events[0]
T_num = te[1] - te[0]
T_ell = 4 * ellipk(np.sin(theta0 / 2) ** 2)
print(f"{theta0:3.1f} {T_num:13.10f} {T_ell:13.10f}"
f" {T_num / (2 * np.pi):.4f}")Çıktı:
θ0 T (sayısal) 4K(sin²(θ0/2)) T / 2π
0.1 6.2871145493 6.2871145493 1.0006
1.0 6.6999756642 6.6999756644 1.0663
2.0 8.3497529262 8.3497529269 1.3289
3.0 16.1555393359 16.1555393724 2.5712İki sütun en az sekiz anlamlı basamakta uyuşuyor. \(\theta_0 = 0{,}1\) radyanda (yaklaşık \(5{,}7°\)) periyot \(2\pi\)’den yalnız binde \(0{,}6\) uzundur ve harmonik yaklaşım çok iyidir. \(\theta_0 = 3\) radyanda (yaklaşık \(172°\)) ise periyot \(2\pi\)’nin \(2{,}57\) katıdır.
scipy.special.ellipk(m) fonksiyonu argüman olarak \(k\)’yi değil \(m = k^2\)’yi alır. Yukarıdaki kodda bu yüzden ellipk(np.sin(theta0 / 2) ** 2) yazdık; ellipk(np.sin(theta0 / 2)) yazmak sessizce yanlış bir periyot verirdi. Özel fonksiyonların hangi gösterimi kullandığını belgelerinden denetleyin.
Son olarak sarkacın faz portresini çizelim. Bir periyodu \(17\)’den kısa olan beş salınımı ve sarkacın tepeyi aşarak döndüğü iki hareketi çözüp her birinin enerjisini yazdırıyoruz. Dönen yörüngeleri \(\theta = -\pi\)’den başlatıp \(\theta = \pi\)’ye varınca bir olay fonksiyonuyla durduruyoruz; ters yönde dönen yörüngeler, denklemin simetrisi sayesinde \((\theta, \omega) \mapsto (-\theta, -\omega)\) yansımasıyla elde ediliyor:
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp
def pendulum(t, z):
return [z[1], -np.sin(z[0])]
def energy(theta, omega):
return omega**2 / 2 - np.cos(theta)
def at_right(t, z):
return z[0] - np.pi # sağ kenar: θ = π
at_right.terminal = True
opts = dict(rtol=1e-10, atol=1e-12, dense_output=True)
fig, ax = plt.subplots(figsize=(7, 4.5))
starts = [(th, 0.0) for th in [0.5, 1.0, 1.5, 2.0, 2.5]]
starts += [(-np.pi, om) for om in [0.8, 1.6]]
for z0 in starts:
if z0[1] == 0.0: # salınım: t = 17 bir tam periyottan uzun
s = solve_ivp(pendulum, (0, 17), z0, **opts)
else: # dönme: θ = π olana kadar
s = solve_ivp(pendulum, (0, 20), z0, events=at_right, **opts)
th, om = s.sol(np.linspace(0, s.t[-1], 400))
ax.plot(th, om, "C0" if z0[1] == 0.0 else "C1")
if z0[1] != 0.0:
ax.plot(-th, -om, "C1") # ters yönde dönen simetrik yörünge
E = energy(th, om)
print(f"başlangıç ({z0[0]:+.3f}, {z0[1]:.1f}): E = {E[0]:+.4f}"
f" sapma = {np.max(np.abs(E - E[0])):.1e}")
th = np.linspace(-np.pi, np.pi, 400)
ax.plot(th, 2 * np.cos(th / 2), "k--", th, -2 * np.cos(th / 2), "k--")
ax.set_xlabel("θ")
ax.set_ylabel("ω = θ'")
fig.savefig("sarkac-faz.png", dpi=150)Çıktı:
başlangıç (+0.500, 0.0): E = -0.8776 sapma = 5.5e-11
başlangıç (+1.000, 0.0): E = -0.5403 sapma = 2.6e-10
başlangıç (+1.500, 0.0): E = -0.0707 sapma = 8.8e-10
başlangıç (+2.000, 0.0): E = +0.4161 sapma = 5.0e-10
başlangıç (+2.500, 0.0): E = +0.8011 sapma = 5.6e-10
başlangıç (-3.142, 0.8): E = +1.3200 sapma = 8.4e-10
başlangıç (-3.142, 1.6): E = +2.2800 sapma = 1.7e-09Her yörüngede enerji sabit kalıyor; en büyük sapma \(1{,}7 \cdot 10^{-9}\)’dur. Kod grafiği sarkac-faz.png dosyasına kaydeder; aşağıdaki şekil aynı grafiğin aynı verilerle çizilmiş hâlidir:
Faz portresinde iki tür hareket görülüyor. \(E < 1\) ise sarkaç tepeye ulaşamaz ve orijin çevresindeki kapalı yörüngelerde salınır. \(E > 1\) ise sarkaç tepeyi aşarak hep aynı yönde döner; \(\theta\) sürekli artar ya da azalır ve şekilde bu yörüngeler \(\theta = \pm\pi\) kenarlarından girip çıkar. İki türü \(E = 1\) seviye eğrisi ayırır: \(\frac{\omega^2}{2} - \cos\theta = 1\) eşitliği \(\omega^2 = 2(1 + \cos\theta) = 4\cos^2\frac{\theta}{2}\) verir, yani bu ayırıcı eğri (separatrix) \(\omega = \pm 2\cos\frac{\theta}{2}\)’dir. \((0, 0)\) kararlı bir denge noktasıdır: Jacobi matrisinin özdeğerleri \(\pm i\)’dir ve yakınındaki yörüngeler onu çevreler. \((\pm\pi, 0)\) ise sarkacın tam tepede durduğu kararsız dengedir: oradaki Jacobi matrisi \(\begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}\)’nin özdeğerleri \(\pm 1\)’dir ve yörüngelerin bir kısmı noktaya yaklaşırken bir kısmı ondan uzaklaşır. \(\theta = \pi\) ile \(\theta = -\pi\) aynı fiziksel konumdur.
13.11 Alıştırmalar
Aşağıdaki alıştırmalarda bölümün yöntemlerini yeni problemlerde deneyeceğiz. Çözümlerdeki bazı kodlar odetools.py modülünü kullanır; modül, çalıştırılan programla aynı klasörde olmalıdır.
Alıştırma 13.1 (Doğrusal Bir Denklemde Euler Yöntemi) \(y' = t - y\), \(y(0) = 1\) probleminin tam çözümü \(y = t - 1 + 2e^{-t}\)’dir. \(y(1)\)’i Euler yöntemiyle \(h = 0{,}1\) ve \(h = 0{,}05\) alarak hesaplayınız ve iki hatanın oranını bulunuz.
Çözüm
Adım 1. \(f(t, y) = t - y\)’dir. Tam değer \(y(1) = 1 - 1 + 2e^{-1}\), yani yaklaşık \(0{,}73575888\)’dir.
Adım 2. \(h = 1/n\) olmak üzere \(n = 10\) ve \(n = 20\) adımla euler fonksiyonunu çağırıyoruz:
import numpy as np
from odetools import euler
def f(t, y):
return t - y
exact = 2 * np.exp(-1.0) # y(1) = 1 - 1 + 2e^(-1)
errs = []
for n in [10, 20]:
t, y = euler(f, 0.0, 1.0, 1 / n, n)
errs.append(abs(y[-1] - exact))
print(f"h = {1 / n:.2f} y(1) ≈ {y[-1]:.8f} hata = {errs[-1]:.3e}")
print(f"tam değer: {exact:.8f} hata oranı: {errs[0] / errs[1]:.3f}")Çıktı:
h = 0.10 y(1) ≈ 0.69735688 hata = 3.840e-02
h = 0.05 y(1) ≈ 0.71697184 hata = 1.879e-02
tam değer: 0.73575888 hata oranı: 2.044Adım 3. Hatalar \(0{,}0384\) ve \(0{,}0188\)’dir; oranları \(2{,}04 \approx 2\)’dir. Adım yarıya inince hata da yarıya iniyor: Euler yöntemi birinci mertebedendir.
Denetim. Bu doğrusal denklemde Euler dizisi kapalı biçimde de yazılabilir: \(y_{i+1} = (1 - h)\,y_i + h\,t_i\) adımından tümevarımla \(y_i = t_i - 1 + 2(1 - h)^i\) bulunur. \(h = 0{,}1\) için \(y_{10} = 2 \cdot 0{,}9^{10} \approx 0{,}69735688\) çıktıyla aynıdır. Tam çözümdeki \(e^{-t}\) çarpanının yerini \((1 - h)^{t/h}\) almıştır ve \(h \to 0\) iken \((1 - h)^{t/h} \to e^{-t}\)’dir.
\(\blacksquare\)
Alıştırma 13.2 (Doğrusal Bir Denklemde RK4) \(y' = t - y\), \(y(0) = 1\) probleminde \(y(1)\)’i RK4 ile \(h = 0{,}1\) ve \(h = 0{,}05\) alarak hesaplayınız. Hataların oranı, RK4’ün dördüncü mertebeden olmasıyla uyumlu mudur?
Çözüm
Adım 1. Tam değer yine \(y(1) = 2e^{-1}\)’dir. Dördüncü mertebeden bir yöntemde adımı yarıya indirmek hatayı kabaca \(2^4 = 16\)’ya bölmelidir.
Adım 2. Bu kez rk4 fonksiyonunu çağırıyoruz:
import numpy as np
from odetools import rk4
def f(t, y):
return t - y
exact = 2 * np.exp(-1.0)
errs = []
for n in [10, 20]:
t, y = rk4(f, 0.0, 1.0, 1 / n, n)
errs.append(abs(y[-1] - exact))
print(f"h = {1 / n:.2f} y(1) ≈ {y[-1]:.10f} hata = {errs[-1]:.3e}")
print(f"tam değer: {exact:.10f} hata oranı: {errs[0] / errs[1]:.2f}")Çıktı:
h = 0.10 y(1) ≈ 0.7357595488 hata = 6.665e-07
h = 0.05 y(1) ≈ 0.7357589223 hata = 3.995e-08
tam değer: 0.7357588823 hata oranı: 16.68Adım 3. Hatalar \(6{,}7 \cdot 10^{-7}\) ve \(4{,}0 \cdot 10^{-8}\)’dir; oran \(16{,}68\)’dir ve \(16\)’ya yakındır. Yanıt evettir. Aynı adım sayısıyla Euler’in hatası \(0{,}038\) idi; RK4’ünki yaklaşık \(58\,000\) kat küçüktür.
\(\blacksquare\)
Alıştırma 13.3 (Heun Yöntemi) Heun yöntemi önce bir Euler adımıyla \(\tilde{y} = y_i + h\,f(t_i, y_i)\) tahminini yapar, sonra iki uçtaki eğimin ortalamasıyla ilerler:
\[y_{i+1} = y_i + \frac{h}{2}\Bigl(f(t_i, y_i) + f(t_{i+1}, \tilde{y})\Bigr).\]
Yöntemi bir heun fonksiyonu olarak yazınız ve lojistik problemde (\(y' = y(1 - y)\), \(y(0) = 0{,}1\), \(T = 4\)) mertebesini adımı yarıya indirerek ölçünüz.
Çözüm
Adım 1. Fonksiyon euler’in döngüsüne ikinci bir eğim ekler. \(f\) yalnız \(t\)’ye bağlıysa adım, \(\int_{t_i}^{t_{i+1}} f\,dt\) için yamuk kuralı olur (bkz. Nümerik Analiz); bu yüzden yöntemin ikinci mertebeden olmasını bekleriz.
Adım 2. Tam değer \(y(4) = 1/(1 + 9e^{-4})\)’tür. \(n = 16, 32, \ldots, 256\) adımla hataları ve ardışık oranları hesaplıyoruz:
import numpy as np
def heun(f, t0, y0, h, n):
"""y' = f(t, y), y(t0) = y0 problemini n Heun adımıyla çözer."""
t = t0 + h * np.arange(n + 1)
y = np.zeros(n + 1)
y[0] = y0
for i in range(n):
k1 = f(t[i], y[i]) # baştaki eğim
k2 = f(t[i] + h, y[i] + h * k1) # Euler tahminindeki eğim
y[i + 1] = y[i] + h / 2 * (k1 + k2)
return t, y
def f(t, y):
return y * (1 - y)
exact = 1 / (1 + 9 * np.exp(-4.0))
prev = None
for n in [16, 32, 64, 128, 256]:
t, y = heun(f, 0.0, 0.1, 4 / n, n)
err = abs(y[-1] - exact)
ratio = f"{prev / err:.3f}" if prev else "-"
print(f"n = {n:3d} hata = {err:.3e} oran = {ratio}")
prev = errÇıktı:
n = 16 hata = 2.150e-03 oran = -
n = 32 hata = 5.402e-04 oran = 3.980
n = 64 hata = 1.356e-04 oran = 3.983
n = 128 hata = 3.400e-05 oran = 3.989
n = 256 hata = 8.512e-06 oran = 3.994Adım 3. Oranlar \(4 = 2^2\)’ye yaklaşıyor: Heun yöntemi ikinci mertebedendir. Adım başına iki \(f\) çağrısıyla Euler (\(1\)) ile RK4 (\(4\)) arasında durur; doğruluğu da öyle: \(n = 256\)’da hatası \(8{,}5 \cdot 10^{-6}\)’dır. Aynı \(n\) için Euler’in hatası \(2{,}9 \cdot 10^{-4}\), RK4’ünki \(5{,}9 \cdot 10^{-11}\)’dir.
\(\blacksquare\)
Alıştırma 13.4 (Aynı Doğruluk İçin Kaç Adım?) \(y' = y\), \(y(0) = 1\) probleminde \(y(1) = e\) değerini \(10^{-6}\)’dan küçük bir hatayla bulmak için Euler ve RK4 yöntemlerinin kaçar adım ve kaçar \(f\) çağrısı gerektirdiğini karşılaştırınız.
Çözüm
Adım 1. RK4 için \(n\)’yi birer birer artırıp hatanın \(10^{-6}\)’nın altına ilk düştüğü adım sayısını arıyoruz.
Adım 2. Euler için \(e - y_n \approx \frac{e}{2n}\) yaklaşıklığını kullanıyoruz (bkz. Örnek 13.2): \(\frac{e}{2n} < 10^{-6}\) koşulu \(n > \frac{e}{2} \cdot 10^6 \approx 1\,359\,141\) verir. Kestirimi, bu sayının iki yanındaki \(n = 1\,350\,000\) ve \(n = 1\,360\,000\) değerleriyle sınıyoruz:
import numpy as np
from odetools import euler, rk4
def f(t, y):
return y
n = 1
while abs(np.e - rk4(f, 0.0, 1.0, 1 / n, n)[1][-1]) >= 1e-6:
n += 1
err = abs(np.e - rk4(f, 0.0, 1.0, 1 / n, n)[1][-1])
print(f"RK4 : n = {n} hata = {err:.2e} f çağrısı = {4 * n}")
for n in [1_350_000, 1_360_000]:
err = abs(np.e - euler(f, 0.0, 1.0, 1 / n, n)[1][-1])
print(f"Euler: n = {n} hata = {err:.4e} f çağrısı = {n}")Çıktı:
RK4 : n = 13 hata = 7.44e-07 f çağrısı = 52
Euler: n = 1350000 hata = 1.0068e-06 f çağrısı = 1350000
Euler: n = 1360000 hata = 9.9937e-07 f çağrısı = 1360000Adım 3. RK4’e \(13\) adım, yani \(52\) çağrı yetiyor. Euler’de \(1\,350\,000\) adım yetmiyor, \(1\,360\,000\) adım yetiyor; kestirim doğru. Euler yaklaşık \(1{,}36\) milyon çağrı ister, yani RK4’ün yaklaşık \(26\,000\) katı. Yüksek mertebenin değeri budur.
\(\blacksquare\)
Alıştırma 13.5 (Lojistik Büyümede Yarı Doygunluk Anı) \(y' = 0{,}5\,y\left(1 - \frac{y}{100}\right)\), \(y(0) = 5\) lojistik denkleminin çözümünün \(50\) değerine ulaştığı anı solve_ivp’nin olay fonksiyonuyla bulunuz ve tam çözümden elde edilen değerle karşılaştırınız.
Çözüm
Adım 1: tam çözüm. \(y' = ry(1 - y/K)\), \(y(0) = y_0\) denkleminin çözümü
\[y(t) = \frac{K}{1 + \left(\frac{K}{y_0} - 1\right)e^{-rt}}\]
dir. \(r = 0{,}5\), \(K = 100\), \(y_0 = 5\) için \(y(t) = \frac{100}{1 + 19e^{-t/2}}\) olur. \(y = 50\) koşulu \(19e^{-t/2} = 1\), yani \(t = 2\ln 19 \approx 5{,}8889\) verir.
Adım 2: olay fonksiyonu. \(y - 50\) fonksiyonu aranan anda sıfır olur; terminal=True ile çözümü orada durduruyoruz:
import numpy as np
from scipy.integrate import solve_ivp
r, K = 0.5, 100.0
def logistic(t, y):
return r * y * (1 - y / K)
def half(t, y):
return y[0] - K / 2
half.terminal = True
sol = solve_ivp(logistic, (0, 50), [5.0], events=half,
rtol=1e-10, atol=1e-10)
print(f"sayısal : t = {sol.t_events[0][0]:.10f}")
print(f"2 ln 19 : t = {2 * np.log(19):.10f}")Çıktı:
sayısal : t = 5.8888779596
2 ln 19 : t = 5.8888779583Adım 3. İki değer sekiz ondalık basamakta uyuşuyor; aradaki \(1{,}3 \cdot 10^{-9}\)’luk fark verdiğimiz hata payı düzeyindedir. Nüfus, doygunluk değerinin yarısına \(t \approx 5{,}889\)’da ulaşır; büyüme hızı da bu anda en büyüktür.
\(\blacksquare\)
Alıştırma 13.6 (Sönümlü Salınıcı) \(y'' + 2y' + 5y = 0\), \(y(0) = 1\), \(y'(0) = 0\) problemini birinci mertebeden bir sisteme çevirip solve_ivp ile çözünüz. \(t = 1, 2, 3\)’teki değerleri tam çözüm \(y = e^{-t}\left(\cos 2t + \frac{1}{2}\sin 2t\right)\) ile karşılaştırınız.
Çözüm
Adım 1. \(u_1 = y\), \(u_2 = y'\) dersek \(u_1' = u_2\) ve \(u_2' = y'' = -2u_2 - 5u_1\) olur; başlangıç vektörü \((1, 0)\)’dır.
Adım 2. Tam çözümü denetleyelim. Karakteristik denklem \(r^2 + 2r + 5 = 0\)’ın kökleri \(-1 \pm 2i\)’dir, dolayısıyla \(y = e^{-t}(A\cos 2t + B\sin 2t)\) olur. \(y(0) = A = 1\) ve \(y'(0) = -A + 2B = 0\) koşulları \(B = \frac{1}{2}\) verir.
Adım 3. Sistemi t_eval ile istenen anlarda çözüyoruz:
import numpy as np
from scipy.integrate import solve_ivp
def damped(t, u):
y, v = u
return [v, -2 * v - 5 * y] # y'' = -2y' - 5y
def exact(t):
return np.exp(-t) * (np.cos(2 * t) + 0.5 * np.sin(2 * t))
t = np.array([1.0, 2.0, 3.0])
sol = solve_ivp(damped, (0, 3), [1.0, 0.0], t_eval=t,
rtol=1e-10, atol=1e-12)
for tk, yk in zip(sol.t, sol.y[0]):
print(f"t = {tk:.0f} sayısal = {yk:+.10f} tam = {exact(tk):+.10f}")Çıktı:
t = 1 sayısal = +0.0141640490 tam = +0.0141640489
t = 2 sayısal = -0.1396720846 tam = -0.1396720846
t = 3 sayısal = +0.0408484245 tam = +0.0408484245Adım 4. Değerler on ondalık basamakta (ilk satırda dokuz basamakta) uyuşuyor. Çözüm \(e^{-t}\) zarfı içinde salınarak sıfıra iner.
\(\blacksquare\)
Alıştırma 13.7 (Doğrusal Sistem ve Matris Üsteli) \(\mathbf{x}' = A\mathbf{x}\), \(\mathbf{x}(0) = (1, 0)\) sisteminde \(A = \begin{pmatrix} 0 & 1 \\ -2 & -3 \end{pmatrix}\)’tür. Kare bir \(A\) matrisinin üsteli (matrix exponential), \(e^x\)’in kuvvet serisinde \(x\) yerine \(tA\) yazılarak \(e^{tA} = \sum_{k \ge 0} \frac{t^k A^k}{k!}\) biçiminde tanımlanır. Çözümün \(t = 1\)’deki değerini solve_ivp ile hesaplayınız ve \(e^{A}\mathbf{x}(0)\) ile karşılaştırınız. Matris üsteli için scipy.linalg.expm fonksiyonunu kullanınız.
Çözüm
Adım 1. \(e^{tA}\) serisi her \(t\) için yakınsar ve terim terim türevlenebilir:
\[\bigl(e^{tA}\bigr)' = A + tA^2 + \frac{t^2}{2!}A^3 + \cdots = A\,e^{tA}.\]
\(e^{0A} = I\) olduğundan \(\mathbf{x}(t) = e^{tA}\mathbf{x}(0)\) fonksiyonu hem \(\mathbf{x}' = A\mathbf{x}\) denklemini hem de başlangıç koşulunu sağlar; yani sistemin çözümüdür. Bu yüzden \(t = 1\)’de \(e^{A}\mathbf{x}(0)\)’ı karşılaştırmalıyız. \(A\)’nın özdeğerleri \(-1\) ve \(-2\)’dir. Sistemin ilk bileşeni \(y'' + 3y' + 2y = 0\) denklemini sağlar ve \(y(0) = 1\), \(y'(0) = 0\) koşullarıyla \(y = 2e^{-t} - e^{-2t}\), \(y' = -2e^{-t} + 2e^{-2t}\) bulunur.
Adım 2. Sağ tarafı A @ x ile tek satırlık bir lambda olarak yazıyoruz:
import numpy as np
from scipy.integrate import solve_ivp
from scipy.linalg import expm
A = np.array([[0.0, 1.0],
[-2.0, -3.0]])
x0 = np.array([1.0, 0.0])
print("özdeğerler:", np.linalg.eigvals(A))
sol = solve_ivp(lambda t, x: A @ x, (0, 1), x0, rtol=1e-10, atol=1e-12)
print("solve_ivp :", sol.y[:, -1])
print("expm(A)x0 :", expm(A) @ x0)
e1, e2 = np.exp(-1.0), np.exp(-2.0)
print("formül :", np.array([2 * e1 - e2, -2 * e1 + 2 * e2]))
print("en büyük fark:", np.max(np.abs(sol.y[:, -1] - expm(A) @ x0)))Çıktı:
özdeğerler: [-1.+0.j -2.+0.j]
solve_ivp : [ 0.6004236 -0.46508832]
expm(A)x0 : [ 0.6004236 -0.46508832]
formül : [ 0.6004236 -0.46508832]
en büyük fark: 1.0809686479262837e-11Adım 3. NumPy özdeğerleri karmaşık sayı türünde döndürdü; sanal kısımları sıfırdır ve özdeğerler \(-1\) ile \(-2\)’dir. Üç yol aynı vektörü veriyor; solve_ivp ile expm arasındaki fark \(1{,}1 \cdot 10^{-11}\)’dir. İki özdeğer de negatif olduğundan çözüm sıfıra söner.
\(\blacksquare\)
Alıştırma 13.8 (Van der Pol Salınıcısının Limit Çevrimi) \(y'' - (1 - y^2)\,y' + y = 0\) van der Pol denkleminin çözümleri, başlangıç koşulundan bağımsız olarak aynı periyodik harekete yaklaşır. Denklemi \((y, y') = (0{,}5;\ 0)\) ve \((3;\ 0)\) başlangıçlarından \([0, 100]\) aralığında çözüp son tam turdaki periyodu ve genliği, yani \(y\)’nin en büyük değerini, ölçünüz.
Çözüm
Adım 1. \(u_1 = y\), \(u_2 = y'\) ile sistem \(u_1' = u_2\), \(u_2' = (1 - u_1^2)\,u_2 - u_1\) olur. Kodda katsayıyı mu = 1.0 adıyla genel bıraktık.
Adım 2. Periyodu, \(y\)’nin aşağıdan yukarıya sıfırdan geçtiği son iki anın farkı olarak ölçüyoruz; genliği de bu son tur boyunca yoğun çıktıdan okuyoruz:
import numpy as np
from scipy.integrate import solve_ivp
mu = 1.0
def vdp(t, u):
y, v = u
return [v, mu * (1 - y**2) * v - y]
def upward(t, u):
return u[0]
upward.direction = 1
for u0 in [[0.5, 0.0], [3.0, 0.0]]:
sol = solve_ivp(vdp, (0, 100), u0, events=upward,
rtol=1e-10, atol=1e-12, dense_output=True)
te = sol.t_events[0]
T = te[-1] - te[-2] # son tam tur
tt = np.linspace(te[-2], te[-1], 2001)
amp = np.max(sol.sol(tt)[0])
print(f"başlangıç {u0}: periyot = {T:.6f} genlik = {amp:.6f}")Çıktı:
başlangıç [0.5, 0.0]: periyot = 6.663287 genlik = 2.008619
başlangıç [3.0, 0.0]: periyot = 6.663287 genlik = 2.008619Adım 3. İki başlangıç da periyodu \(6{,}663287\), genliği \(2{,}008619\) olan aynı harekete varıyor. Faz düzleminde bu hareket, öteki yörüngeleri kendine çeken yalıtılmış bir kapalı yörüngedir ve limit çevrim (limit cycle) adını alır. Lotka–Volterra modelinden farkı buradadır: orada her başlangıç kendi kapalı yörüngesinde kalıyordu. Van der Pol denkleminde \(|y| < 1\) iken \(-(1 - y^2)\,y'\) terimi harekete enerji katar, \(|y| > 1\) iken enerji çeker; küçük salınımlar büyür, büyükler söner.
\(\blacksquare\)
Alıştırma 13.9 (Lotka–Volterra Modelinde Ortalama Nüfuslar) Lotka–Volterra modelinde (\(a = 1\), \(b = 0{,}5\), \(c = 0{,}75\), \(d = 0{,}25\)) bir periyot boyunca av ve avcı nüfuslarının zaman ortalamalarının yörüngeden bağımsız olarak denge değerlerine, \(c/d = 3\) ve \(a/b = 2\)’ye eşit olduğunu gösteriniz. Sonucu \(y(0) = 2\) ve \(x(0) = 4{,}5;\ 6;\ 8\) yörüngelerinde sayısal olarak doğrulayınız.
Çözüm
Adım 1: ispat. \(T\) periyot olsun. İkinci denklem \(\frac{y'}{y} = -c + dx\), yani \((\ln y)' = -c + dx\) biçiminde yazılabilir. \([0, T]\) üzerinde integralleyince, \(y(T) = y(0)\) olduğundan
\[0 = \ln y(T) - \ln y(0) = -cT + d\int_0^T x\,dt\]
bulunur; buradan \(\frac{1}{T}\int_0^T x\,dt = \frac{c}{d}\) çıkar. Aynı biçimde \((\ln x)' = a - by\) eşitliği \(\frac{1}{T}\int_0^T y\,dt = \frac{a}{b}\) verir.
Adım 2: sayısal doğrulama. Periyodu olay fonksiyonuyla ölçüp ortalamaları yoğun çıktı üzerinde scipy.integrate.quad ile (bkz. Sayısal Türev ve İntegral) hesaplıyoruz:
import numpy as np
from scipy.integrate import solve_ivp, quad
def lotka(t, z, a, b, c, d):
x, y = z
return [a * x - b * x * y, -c * y + d * x * y]
def prey_down(t, z, a, b, c, d):
return z[0] - c / d
prey_down.direction = -1
p = (1.0, 0.5, 0.75, 0.25)
for x0 in [4.5, 6.0, 8.0]:
sol = solve_ivp(lotka, (0, 30), [x0, 2.0], args=p, events=prey_down,
rtol=1e-10, atol=1e-12, dense_output=True)
t1, t2 = sol.t_events[0][:2] # bir tam periyot
T = t2 - t1
mx = quad(lambda t: sol.sol(t)[0], t1, t2)[0] / T
my = quad(lambda t: sol.sol(t)[1], t1, t2)[0] / T
print(f"x(0) = {x0}: T = {T:.4f} ort. x = {mx:.8f}"
f" ort. y = {my:.8f}")Çıktı:
x(0) = 4.5: T = 7.3556 ort. x = 3.00000000 ort. y = 2.00000000
x(0) = 6.0: T = 7.5834 ort. x = 3.00000000 ort. y = 2.00000000
x(0) = 8.0: T = 7.9976 ort. x = 3.00000000 ort. y = 2.00000000Adım 3. Periyotlar farklı olsa da ortalamalar sekiz ondalık basamakta \(3\) ve \(2\)’dir. Nüfuslar dengeden uzakta salınsa bile ortalamaları hep denge değerleridir.
\(\blacksquare\)
Alıştırma 13.10 (Sönümlü Sarkaç Kaç Tur Atar?) Sürtünmeli bir sarkacın denklemi \(\theta'' = -\sin\theta - 0{,}2\,\theta'\) olsun. Sarkaç en alttan (\(\theta(0) = 0\)) \(\theta'(0) = 4\) hızıyla itiliyor. Sarkaç tepe noktasından kaç kez geçer ve sonunda hangi açıda durulur?
Çözüm
Adım 1. Tepe noktası \(\theta = \pi\) konumudur; sarkaç döndükçe \(\theta\) artmaya devam eder, bu yüzden tepeden geçişler \(\theta = \pi, 3\pi, 5\pi, \ldots\) anlarıdır. Bu açıların hepsinde \(\cos\frac{\theta}{2} = 0\) olduğundan tek bir olay fonksiyonu, cos(theta / 2), hepsini yakalar.
Adım 2. Başlangıç enerjisi \(E = \frac{4^2}{2} - \cos 0 = 7\)’dir ve tepeyi aşmak için \(E > 1\) gerekir. Sürtünme enerjiyi yavaş yavaş tüketir; denklemi uzun bir süre boyunca çözüp geçişleri sayıyoruz. Geçiş anlarındaki açıları sol.y_events[0] dizisinden okuyoruz; bu dizinin her satırı bir geçiş anındaki \((\theta, \theta')\) vektörü olduğundan açılar ilk sütundadır:
import numpy as np
from scipy.integrate import solve_ivp
def damped_pendulum(t, z):
theta, omega = z
return [omega, -np.sin(theta) - 0.2 * omega]
def top(t, z):
return np.cos(z[0] / 2) # θ = π, 3π, 5π, ... iken sıfır
sol = solve_ivp(damped_pendulum, (0, 100), [0.0, 4.0], events=top,
rtol=1e-10, atol=1e-12)
print("tepeden geçiş anları:", np.round(sol.t_events[0], 4))
print("o anlardaki θ / π :", np.round(sol.y_events[0][:, 0] / np.pi, 6))
theta_end, omega_end = sol.y[:, -1]
print(f"θ(100) / 2π = {theta_end / (2 * np.pi):.6f}")
print(f"ω(100) = {omega_end:.1e}")Çıktı:
tepeden geçiş anları: [0.936 3.6575]
o anlardaki θ / π : [1. 3.]
θ(100) / 2π = 1.999975
ω(100) = -3.0e-05Adım 3. Sarkaç tepeden iki kez geçiyor: \(t \approx 0{,}94\)’te \(\theta = \pi\)’den, \(t \approx 3{,}66\)’da \(\theta = 3\pi\)’den. Sonra enerjisi \(1\)’in altına düşer ve \(\theta = 4\pi\) çevresinde sönümlü salınım yaparak durulur: \(\theta(100)/2\pi \approx 2\)’dir. Sarkaç iki tam tur atıp yine en alttaki konumda durulmuştur; \(4\pi\) açısı fiziksel olarak \(0\) ile aynı konumdur.
\(\blacksquare\)
Bu bölümde başlangıç değer problemlerini sayısal olarak çözmeyi öğrendik: Euler ve RK4 yöntemlerini kendimiz yazdık, mertebelerini log-log grafikte ölçtük, solve_ivp ile uyarlamalı adımlı ve katı denklemlere uygun yöntemleri kullandık, olay fonksiyonlarıyla belirli anları yakaladık ve faz düzleminde Lotka–Volterra modelini ve sarkacı inceledik. Buraya kadarki bütün modeller deterministikti: aynı başlangıç her seferinde aynı sonucu verdi. Bir sonraki bölüm Rastgele Sayılar ve Monte Carlo Yöntemleri, rastlantıyı hesaba katıyor: rastgele sayı üretmeyi, olasılık yasalarını simülasyonla görmeyi ve integralleri rastgele örneklemle hesaplamayı anlatıyor.