Ders 2
Gerçek dünyada çoğu ODE sayısal (numerical) yöntemlerle çözülür. Bilgisayarınızda bir ODE'yi grafiklettirdiğiniz zaman da aslında arka planda bilgisayar o denklemi sayısal olarak çözmekte ve sonucu grafiklettirmektedir. Bir başlangıç değerli (initial value) probleminin formunu yazalım:
$$ y' = f(x,y) $$
$$ y_1(x_o) = y_o $$
Problemde ilk satır ODE, ikinci satır bu ODE'nin başlangıç değeri, $x_0$ ve $y_0$ sabit değerler.
Sayısal olarak (mesela Euler yöntemiyle) bu denklemi çözmek ne demektir? Alttaki çizime bakalım, $y_0$'dan değerinden başlıyoruz, bu noktada $x_0,y_0$ noktasındaki eğimi $y'$ ile hesaplıyoruz, ve bu eğim bize $y$'nin olacağı bir sonraki yeri söylüyor. Bu eğim ile yukarı ya da aşağı çıkıyoruz, ve bunu devam ettiriyoruz, ta ki bir sonuca gelinceye kadar.
Peki bahsedilen ODE bağlamında bir yere "gitmek" ne demektir? Takip ettiğimiz, cevap olarak odaklandığımız $y$ fonksiyonudur. Unutmayalım ki bu fonksiyonun tam hali $y(x)$, yani $y$, $x$'in bir fonksiyonu, $y$'nin türevi $y$'nin $x$'e göre türevi demek. Türev eğim demektir, eğim fonksiyonun o noktadaki kabaca, yaklaşıksal bir yönüdür. O yönü takıp edersek o fonksiyonu aşağı yukarı takip ediyoruz demektir.
Grafikte $h$ basamak mesafesi, yani $x$ üzerinde yaptığımız sabit zıplama mesafesi. Bu kordinatta hangi aralıkla zıplıyoruz? 0.1 mi, 1 mi, 5 mi? Bunun seçimini biz yapıyoruz.
Euler Denklemleri bir adım için şöyle tanımlıdır:
$$ x_{n+1} = x_n + h $$
$$ y_{n+1} = y_n + hA_n $$
$$ A_n = f(x_n, y_n) $$
Bu basamaklar özyineli (recursive) olarak tanımlanır, bir sonraki adım, bir önceki adımın değerlerini kullanır.
Örnek
$$ y' = x^2 - y^2 $$
$$ y_1(0) = 1 $$
$$ h = 0.1 $$
Üstteki formül temel (elementary) fonksiyonlar kullanılarak çözülemez. O yüzden Euler'in yöntemi gibi bir sayısal çözüm burada uygun olur.
\begin{array}{ccccc} n & $x_n$ & $y_n$ & $A_n$ & $hA_n$ \\ \hline 0 & 0 & 1 & -1 & -0.1 \\ \hline 1 & .1 & .9 & -.80 & -0.08 \\ \hline 2 & .2 & .82 & & \end{array}
$y$ için eriştiğimiz sonuç .82 değeridir. Şimdi şunu soralım: Bu cevap çok yukarıda mı çok aşağıda bir cevap mı? Pür sayısal sonuçlarda karşılaşılan bir problem budur, gerçek cevabı analitik olarak bilmediğimiz için ona ne kadar yaklaşıp yaklaşmadığımız. Cevabı geometrik olarak verelim. Eğer çözüm bir düz çizgi olsaydı, Euler metodu her bu çizgi üzerinde hep doğru cevabı veriyor olurdu.
Eğer çözüm dışbükey (concave) olsaydı, üstteki gibi Euler metodu çok aşağı düşecekti. İlk adımda fazla aşağı inecekti, ve sonra bu hatadan dönemeyecek, hep esas fonksiyona uzak kalacaktı. İçbükey olunca benzer şekilde, ama fazla yukarıda kalacaktı.
Peki elimizde bir analitik cevap olmadığına göre cevabın düşbükey (convex) mi içbükey mi (concave) olup olmadığını nereye bakarak anlayacağız? Calculus tekrar hızır gibi imdada yetişiyor. İkinci türevi hatırlayalım: Eğer $y" > 0$ ise birinci türev sürekli artıyor demektir, yani $y$ dişbükeydir. Eğer $y" < 0$ ise tam tersi. Fakat hala bir problem var, analitik fonksiyon yok ise ikinci türevi nasıl hesaplayacağız? Cevap: Diferansiyel fonksiyonun kendisini kullanarak.
$y' = x^2 - y^2$'nin türevini alırsak,
$y" = 2x - 2yy'$ sonucunu elde ederiz (türev alırken zincirleme kanununu kullandığımıza dikkat).
O zaman başlangıç noktası $(0,1)$ de $y"$ nedir? $y'(0) = -1$, $y"= 2 \cdot 0 - 2\cdot 1 = 2$. Demek ki çözüm başlangıçta dişbükey, ve Euler çözümünü uzun süreli takip etmezsek çözümün çok altında kalabiliriz.
Tabii ki çözüm bir dış bir iç olarak sürekli değişen, dalgalı bir yapıda olabilir, bu da mümkün. Burada asıl göstermek istediğimiz diferansiyel denklemin kendisini kullanarak çözüm hakkında analitik hiçbir şey bilmeden onun hakkında nasıl bilgi edinebileceğimizi görmektir.
Hata Analizi
Hata analizi Euler'in çözüme ne kadar uzak kaldığının hesabıdır, yani $e$ sayısını hesaplamaktır. Bu değerin mutlak değeri (absolute value) kullanılır.
Daha iyi sonuçlar için daha küçük $h$ basamakları kullanılabilir, o zaman sonuca daha yakın kalabiliriz. O zaman $e$'nin $h$'ye bağlı olduğunu söyleyebiliriz. Formülsel olarak bu ifade suna benzer:
$$ |e \sim c_1 h| $$
Buna ifadeye göre Euler metotu birinci derece bir metottur denir, bu derecenin ODE'nin derecesiyle alakası yok, $h$'nin üstteki formülde hangi üstel formde olduğuyla alakalı. Birincil derecede bir ilişki mesela basamağı yarısına indirince hatayı yarısına indirirmek demektir.
Euler metodundan daha bir yöntem bulmak demek, eğimi daha iyi hesaplayan bir yöntem bulmak demektir. Eğer hatada rol oynayan en önemli faktör eğim olduğuna göre, daha az hata için daha iyi eğim hesaplamak mantıklı olacaktır.
Daha iyi eğim nasıl hesaplanır? Diyelim ki tek bir zıplama yerine iki kere zıpladık. Dişbükey durumda birinci zıplamada çok aşağı, ikincide biraz daha yukarı gidiyor olurduk, o zaman bunların ortalamasını alırsak, daha iyi bir eğim elde edebilirdik.
$$ x_{n+1} = x_n + h $$
$$ \hat{y}_{n+1} = y_n + h A_n $$
$$ B_n = f(x_{n+1},\hat{y}_{n+1}) $$
$$ y_{n+1} = y_n + h(\frac{A_n+B_n}{2}) $$
Niye şapkalı $y$ yani $\hat{y}$ kullandık? Çünkü ortalama hesaplamak için aslında n+1 noktasında geçici bir değer hesaplıyoruz, bu geçiciliği göstermek için $\hat{y}$ ifadesini kullandık.
Bu metot Heun, Geliştirilmiş (Improved) Euler, Değiştirilmiş (Modified) Euler, RK2 gibi isimlerle anılır. RK2 Runga-Kutta'nın kısaltması, '2' ibaresi bu yöntemin ikinci dereceli bir metot olmasıdır. Yani
$$ e \sim c_2 h^2 $$
Basamağı yarısına indirmek hatayı dörtte birine indirmek demektir. O zaman niye bu metot her yerde kullanılmıyor? Çünkü RK2 ile eğim iki kere hesaplanıyor, Euler ile bir kere, yani sayısal kod iki kat daha fazla çalışmak, yani daha yavaşlamak zorundadır.
RK4 te var, bu dördüncü seviyede bir metot. Eğim şöyle hesaplanır:
$$ \frac{A_n + 2B_n + 2C_n + D_n}{6} $$
Sayısal Yöntemlerde Bazı Tehlikeli Noktalar
1) Hoca ödevde bizim keşfetmemizi istedi
2) Şu basit denkleme bakalım: $y' = y^2$. Değişkenleri ayıralım ve analitik cevabı $y=\frac{1}{c}$. O zaman bu denklemi sayısal olarak çözerken alttaki grafik takip ediliyor olacak.
Diyelim ki $y(0) = 1$'den başladık ve $y(2)$'i bulacağız. Sağa doğru yavaş yavaş gidiyoruz ama problem, bu fonksiyon $y(1)$ değerinde sonsuza gidiyor. Demek ki adım adım sağa giden sayısal çözüm o noktayı hiçbir zaman aşamayacaktır, sonsuzlukta kaybolacaktır. Bu tehlikeli noktayı önceden tahmin edemez miydik? Hayır. Üstteki diferansiyel denklemin her çözümünün kendine has bir eşsizsel (singularity) noktası vardır ve bundan sadece kendisi haberdadır.
Ekler
Euler yöntemininin çoğu kaynakta ODE sistemini "entegre ettiği" söylenir; hakikaten de bir entegrasyon yapıyoruz, yani $dy$ değişim formülünü $y$'ye uygulayarak adım adım ilerliyoruz, $y$'nin değerlerini bulmuş oluyoruz, bu bir entegrasyon işlemidir. Euler fonksiyonu ile analitik fonksiyonları sayısal entegre etmiş oluyoruz.
Taylor Açılımı ile ODE çözümü [1, sf. 395]
Range-Kutta ile atılan adımlar üzerinde oynama yapılarak ilerleme sağlanmaya uğraşılıyor. Fakat aslında Taylor açılımı kullanarak ta ODE çözmek mümkün. Bu metot çoğunlukla teşvik edilmiyor çünkü bazı denklemlerin 2., 3. türevini sembolik olarak almak çok zor. Fakat otomatik türev ile, bkz [2], bu türevler rahatlıkla alınabilir. Altta gösterme amaçlı olarak sembolik alınabilen türevli ODE gösteriyoruz.
ODE bir $y(t)' = f(x,y)$ sistemidir,
$$ y(t)' = t + y , \quad y(1) = 2$$
Bu sistemin çözümü biliniyor,
$$ y(t) = -t -1 + 4e^{-1}e^t$$
Eğer çözüm bilinmeseydi, sayısal olarak çözmek gerekecekti. $a$ yakınında $g(x)$'in Taylor açılımını analitik olarak biliyoruz,
$$ g(x) = g(a) + g'(a)(x-a) + g"(a) \frac{(x-a)^2}{2}$$
Eğer $x$'i $h$ adım atarak $a$'den gelinen nokta olarak belirtirsek
$$ x=a+h$$
$$ g(a+h) = g(a) + g'(a)h + g"(a)\frac{h^2}{2}$$
Bu açılım özyineli bir hesabı ima ediyor; yani $h$ adımı sonrasında elde edilen yeni $g$'nin bir önceki adımdaki değerler ile nasıl bağlantılı olduğunu gösteriyor. O zaman $g$ yerine üstteki $f$ tanımını geçirirsek,
$$ y_{k+1} = y_k + h \frac{d}{dt} f(t_k,y_l) + \frac{h^2}{2} \frac{d^2}{dt} f(t_k,y_l)$$
$$ = y_k + hf(t_k,y_k) + \frac{h^2}{2} \frac{d}{dt} f(t_k,y_l)$$
Aslında Euler yöntemi olarak bilinen yöntem üstteki formülde ilk iki terimin tutulup gerisinin atılmasıyla elde edilmiştir. O zaman daha iyi hesap için üstteki tüm terimleri (hatta açılımı devam ettirerek daha bile fazlasını) kullanabilirdik. 3. terimdeki diferansiyel hesabı için, tam diferansiyel (total derivative) tanımından,
$$ df = \frac{\partial f}{\partial t} dt + \frac{\partial f}{\partial y} dy $$
Yani
$$ \frac{df}{dt} = \frac{\partial f}{\partial t} + \frac{\partial f}{\partial y} \frac{dy}{dt} $$
$y'(t)=dy/dt$ olduğuna göre $t + y$'ye eşittir, onu üstte yerine koyarız, kısmi türevleri de hesapladıktan sonra,
$$ \frac{df}{dt} = 1 + t + y $$
O zaman
$$ = y_k + hf(t_k,y_k) + \frac{h^2}{2} (1 + t_k + y_k) $$
İlk önce normal Euler gösterelim,
import pandas as pd
y1 = 2; t = 1.; N = 5; h = 0.1
vals = []
for i in range(N):
yreal = -t-1+4*(1./np.e)*np.exp(t)
vals.append([t, y1, yreal])
y1 = y1 + (t + y1)*h
t = t + h
df = pd.DataFrame(vals)
df.columns = ['t','euler', 'analitik']
df = df.set_index('t')
print df
df.plot()
plt.savefig('2_6.png')
euler analitik
t
1.0 2.0000 2.000000
1.1 2.3000 2.320684
1.2 2.6400 2.685611
1.3 3.0240 3.099435
1.4 3.4564 3.567299
Oldukça büyük bir sapma var. Şimdi tüm Taylor terimlerini kullanalım,
import pandas as pd
y2 = 2; t = 1.; N = 5; h = 0.1
vals = []
for i in range(N):
yreal = -t-1+4*(1./np.e)*np.exp(t)
vals.append([t, y2, yreal])
y2 = y2 + (t + y2)*h + (1+t+y2)*(h**2/2)
t = t + h
df = pd.DataFrame(vals)
df.columns = ['t','taylor', 'analitik']
df = df.set_index('t')
print df
df.plot()
plt.savefig('2_7.png')
taylor analitik
t
1.0 2.000000 2.000000
1.1 2.320000 2.320684
1.2 2.684100 2.685611
1.3 3.096931 3.099435
1.4 3.563608 3.567299
Kaynaklar
[1] Ackleh, Classical and Modern Numerical Analysis Theory
[2] Bayramlı, Bilgisayar Bilim, Yapay Zeka, Otomatik Türev Almak
Yukarı