BDP - Runge-Kutta Yöntemleri

Yazar

Taygun Bulmus

Yayınlanma Tarihi

26 Ağustos 2024

Runge-Kutta Yöntemleri

Euler yöntemi için Taylor serisini birinci dereceden türevli terimden kesmiştik. Eğer Taylor serisinde daha yüksek mertebelerden türevler kullanılırsa, bu yöntemlere Runge-Kutta yöntemleri adı verilir.

Dikkat

Yüksek mertebeden türevler kullanmak demek, yüksek mertebe diferansiyel denklem çözümü anlamına gelmez. Örneğin, birinci dereceden diferansiyel denklemi, Taylor serisinin üçüncü dereceden teriminden keserek çözebiliriz. Runge-Kutta yöntemleri, Euler yönteminin bir genelleştirilmesidir.

İkinci Derece Runge-Kutta Yöntemi

Ana denklemi yazalım.

\[ y'(x)=f(y,x) \]

Euler yöntemini hatırlayalım. \(y'(x)=f(y,x)\) bağıntısını da kullanalım.

\[ y(x+h)= y(x)+hy'(y,x) = y(x)+hf(y,x) \]

Taylor serisinin ikinci mertebesinden kesmeden önce en genel formülü yazalım. Bu formül Runge-Kutta yöntemlerinin genel formülüdür.

\[ y(x+h)= y(x)+c_{0}hF(y,x) + c_{1}hf(x+ph,y+qhF(y,x)) \]

Burada \(c_{0},c_{1},p,q\) birer katsayıdır. Taylor serisine dönelim.

\[ \begin{align*} y(x+h) &= y(x) + y'(x)h + \frac{1}{2}y''(x)h^{2}\\ &= y(x) + F(y,x)h + \frac{1}{2}f'(y,x)h^{2} \end{align*} \]

Burada \(f'(y,x)\) aşağıdaki gibi yazılır.

\[ \begin{align*} f'(y,x)&= \frac{\partial f}{\partial x} + \frac{\partial f}{\partial y}f \\ &= \frac{\partial f}{\partial x} + \frac{\partial f}{\partial y}y' \end{align*} \]

Vektörize edilmiş \(F\) için ise,

\[ F'(y,x)= \frac{\partial F}{\partial x} + \sum_{i=0}^{n-1}\frac{\partial F}{\partial y_{i}}F_{i}(y,x) \]

şeklinde yazılır. Taylor açılımından elde ettiğimiz denklemde yerine koyarsak,

\[ y(x+h)= y(x) + F(y,x)h + \frac{h^{2}}{2}\frac{\partial F}{\partial x} + \sum_{i=0}^{n-1}\frac{\partial F}{\partial y_{i}}F_{i}(y,x) \]

ifadesini elde ederiz. Şimdi genelleştirilmiş formülü, yani \(y(x+h)= y(x)+c_{0}hF(y,x) + c_{1}hF(x+ph,y+qhF(y,x))\) terimindeki \(F(x+ph,y+qhF(y,x))\) terimi açalım.

\[ F(x+ph,y+qhF(y,x)) = F(y,x) + \frac{\partial F}{\partial x}ph + \sum_{i=0}^{n-1}\frac{\partial F}{\partial y_{i}}qhF_{i}(y,x) \]

  • Yukarıdaki açılımı genelleştirilmiş formülde yerine koyalım.

\[ y(x+h) = y(x) + (c_{0}+c_{1})F(y,x)h+ c_{1}\left[\frac{\partial F}{\partial x}ph + qh \sum_{i=1}^{n-1}\frac{\partial F}{\partial y_{i}}F_{i}(y,x)\right] \]

Burada elde ettiğimiz formül ile Taylor serisinde elde ettiğimiz formülü karşılaştıralım. Eğer

\[ c_{0}+c_{1}=1 \quad \text{ve} \quad c_{1}p=c_{1}q=\frac{1}{2} \]

olursa, iki formül de aynı olur. Buradaki \(c_{0}, c_{1},q,p\) parametreleri farklı farklı seçilebilir. Bazı farklı seçimler için özel isimlendirme yapılır.

İsim \(c_{0}\) \(c_{1}\) p q
Değiştirilmiş (Modified) Euler 0 1 1/2 1/2
Heun Yöntemi 1/2 1/2 1 1
Ralston Yöntemi 1/3 2/3 3/4 3/4
  • Tüm bu sınıflandırma 2. derece Runge-Kutta yöntemleri altındadır.
  • Bu yöntemler birbirlerinden üstün değildir.
  • Birinci derece Runge-Kutta yöntemi olan Euler yönteminde ise \(c_{0}=1\) ve diğer terimler sıfırdır.

Genelleştirmek adına değiştirilmiş Euler motdunu aşağıdaki gibi yazabiliriz.

\[ \vec{y}(x+h) = \vec{y}(x) + \vec{F}\left[x+\frac{h}{2},\vec{y}+\frac{h}{2}\vec{F}(x,\vec{y}) \right] \]

Runge-Kutta metodları için \(\vec{K}_{i}\) vektörleri tanımlanır.

\[ \vec{K}_{0} = h\vec{F}(x,\vec{y}) \]

  • En popüler Runge-Kutta diferansiyel denklem çözme yöntemi dördüncü dereceden Runge-Kutta yöntemidir.

Dördüncü Derece Runge-Kutta Yöntemi

Dördüncü Derece Runge-Kutta yöntemi için \(\vec{K}_{i}\) vektörleri tanımlanır.

\[ \begin{align*} \vec{K}_{0} &= h\vec{F}(x,\vec{y})\\ \vec{K}_{1} &= h\vec{F}\left(x+\frac{h}{2},\vec{y}+\frac{1}{2}\vec{K}_{0}\right)\\ \vec{K}_{2} &= h\vec{F}\left(x+\frac{h}{2},\vec{y}+\frac{1}{2}\vec{K}_{1}\right)\\ \vec{K}_{3} &= h\vec{F}\left(x+h,\vec{y}+\vec{K}_{2}\right) \end{align*} \]

Sonuç olarak bir sonraki adım aşağıdaki gibi yazılır.

\[ \vec{y}(x+h)= \vec{y}(x) + \frac{1}{6}\left(\vec{K}_{0}+2\vec{K}_{1}+2\vec{K}_{2}+\vec{K}_{3}\right) \]

Katsayı tablosu oluşturalım.

\(K_{0}\) \(K_{1}\) \(K_{2}\) \(K_{3}\)
\(hF(y,x)\) \(hF\left(y+\frac{1}{2}K_{0}, x+\frac{h}{2}\right)\) \(hF\left(y+\frac{1}{2}K_{1}, x+\frac{h}{2}\right)\) \(hF\left(y+K_{2}, x+h\right)\)

Katsayıları kullanarak bir sonraki adım hesaplanır.

\[ y(x+h)= y(x) + \frac{1}{6}(K_{0}+2K_{1}+2K_{2}+K_{3}) \]

Python kodunu yazalım. add_coz_rk4_sistem adlı fonksiyon bilYonMod.py modülünde tanımlıdır.

Not

add_coz_rk4_sistem fonksiyonunu inceleyiniz.

Alıştırma 1

Birbirine bağlı iki adet diferansiyel denklemi çözmek için Euler yöntemini, 4. Mertebe Runge-Kutta ve scipy.integrate.odeint yöntemlerini kullanınız. Toplamda 10 adım kullanın.

Denklem sistemi aşağıdaki gibidir.

\[ \begin{align*} \frac{d}{dx}y(v(x), x) &= v(x) \\ \frac{d}{dx}v(x) &= x \end{align*} \]

Analitik çözümü yazabiliriz. Önce \(v(x)\)’i çözelim.

\[ \frac{d}{dx}v(x) = x \Rightarrow v(x) = \frac{x^{2}}{2} + C \]

Başlangıç koşulunu, \(v(0)=0\) koşulunu kullanarak \(C\) sabitini bulalım.

\[ v(0) = 0 \Rightarrow C = 0 \Rightarrow v(x) = \frac{x^{2}}{2} \]

Şimdi \(y(x)\)’i çözelim.

\[ \begin{align*} \frac{d}{dx}y(v(x), x) =& v(x) \\ y(x) =& \int v(x)dx \\ =& \int \frac{x^{2}}{2}dx = \frac{x^{3}}{6} + C \end{align*} \]

Başlangıç koşulunu, \(y(0)=1\) koşulunu kullanarak \(C\) sabitini bulalım.

\[ y(0) = 1 \Rightarrow C = 1 \Rightarrow y(x) = \frac{x^{3}}{6} + 1 \]

Analitik çözüm aşağıdaki gibi olur.

\[ y(x)= \frac{x^{3}}{6} + 1 \quad \text{ve} \quad v(x) = \frac{x^{2}}{2} \]

Çözüm

################################################
## Modül yolunu varsayılan yol olarak ekleme ve modülü içe aktarma
import os
import sys
# Bu dosyanın bulunduğu dizini al
current_dir = os.path.abspath('')
# 3 üst dizine çık
module_dir = os.path.join(os.path.abspath(os.path.join(current_dir, os.pardir, os.pardir, os.pardir)), 'moduller')
# moduller dizinini yol olarak ekle
sys.path.append(module_dir)
# bilYonMod.py modülünü içe aktar
import bilYonMod as bym
################################################

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as spInt
# Fonksiyon
def fonk_yVek_x(yVek, x):
    dydx = yVek[1]
    return np.array([dydx, x])
# Başlangıç koşulları
y0 = 1
v0 = 0
x0 = 0
xSon = 10
n = 10
# Çöz, Euler
xTumEuler, yTum_VekEuler = bym.add_coz_euler_sistem(fonk_yVek_x, x0, xSon, np.array([y0, v0]), n)
# Çöz, 4. Mertebe Runge-Kutta
xTumRK4, yTum_VekRK4 = bym.add_coz_rk4_sistem(fonk_yVek_x, x0, xSon, np.array([y0, v0]), n)
# Çöz, scipy.integrate.odeint
xTumScipy= np.linspace(x0, xSon, n)
yTumScipy= spInt.odeint(fonk_yVek_x, np.array([y0, v0]), xTumScipy)
# Analitik çözüm
yAnalitik = xTumEuler**3/6 + 1
vAnalitik = xTumEuler**2/2
# Çiz
fig, ax = plt.subplots(1, 2)
# Sol grafik
ax[0].plot(xTumEuler, yTum_VekEuler[0], lw='5', label='y(x) Euler', color='red')
ax[0].plot(xTumRK4, yTum_VekRK4[0], lw='3', label='y(x) RK4', color='green')
ax[0].plot(xTumScipy, yTumScipy[:,0], lw='1', label='y(x) Scipy', color='blue')
ax[0].plot(xTumEuler, yAnalitik, label='y(x) Analitik', color='black', linestyle='--')
ax[0].set_xlabel('x')
ax[0].set_ylabel('y(x)')
ax[0].legend()
# Sağ grafik
ax[1].plot(xTumEuler, yTum_VekEuler[1], lw='5', label='v(x) Euler', color='blue')
ax[1].plot(xTumRK4, yTum_VekRK4[1], lw='3', label='v(x) RK4', color='green')
ax[1].plot(xTumScipy, yTumScipy[:,1], lw='1', label='v(x) Scipy', color='red')
ax[1].plot(xTumEuler, vAnalitik, label='v(x) Analitik', color='cyan', linestyle='--')
ax[1].set_xlabel('x')
ax[1].set_ylabel('v(x)')
ax[1].legend()
plt.show()

Alıştırma 2

Aşağıdaki diferansiyel denklem sistemini tüm bildiğiniz çözüm yöntemleri ile çözün  [1].

\[ \begin{align*} m_{1}\ddot{x}_{1}+b_{1}\dot{x}_{1}+ k_{1}(x_{1} - L_{1})-k_{2}(x_{2}-x_{1}-L_{2})=0 \\ m_{2}\ddot{x}_{2}+b_{2}\dot{x}_{2}+ k_{2}(x_{2}- x_{1}-L_{2})=0 \end{align*} \]

  • \(m_{1}= 1\) kg, \(m_{2}=1.5\) kg: Cisimlerin kütlesi
  • \(k_{1}= 8\) N/m, \(k_{2}=40.0\) N/m: Yayların esneklik katsayısı
  • \(b_{1}= 0.8\) Ns/m, \(b_{2}=0.5\) Ns/m: Sürtünme katsayısı
  • \(L_{1}= 0.5\) m, \(L_{2}=1.0\) m: Yayların kuvvet yokkenki uzunluğu

Başlangıç koşulları:

  • \(x_{1}=0.5\) m
  • \(v_{1}=0.0\) m/s
  • \(x_{2}=2.25\) m
  • \(v_{2}=0.0\) m/s

Çözüm

################################################
## Modül yolunu varsayılan yol olarak ekleme ve modülü içe aktarma
import os
import sys
# Bu dosyanın bulunduğu dizini al
current_dir = os.path.abspath('')
# 3 üst dizine çık
module_dir = os.path.join(os.path.abspath(os.path.join(current_dir, os.pardir, os.pardir, os.pardir)), 'moduller')
# moduller dizinini yol olarak ekle
sys.path.append(module_dir)
# bilYonMod.py modülünü içe aktar
import bilYonMod as bym
################################################

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as spInt
# Global değişkenler
m1 = 1.0
m2 = 1.5
k1 = 8.0
k2 = 40.0
L1 = 0.5
L2 = 1.0
b1 = 0.8
b2 = 0.5
# Fonksiyonlar
def fonkVek(yVek, t):
    x1= yVek[0]
    v1= yVek[1]
    x2= yVek[2]
    v2= yVek[3]
    return np.array([v1\
        , (-b1 * v1 - k1 * (x1 - L1) + k2 * (x2 - x1 - L2)) / m1 \
        , v2 \
        , (-b2 * v2 - k2 * (x2 - x1 - L2)) / m2])
# Başlangıç koşulları
yVek0= np.array([0.5, 0.0, 2.25, 0.0])
t0=0
tSon=50
n=100
# Çöz, Euler
tTumEuler, yTum_VekEuler = bym.add_coz_euler_sistem(fonkVek, t0, tSon, yVek0, n)
# Çöz, 4. Mertebe Runge-Kutta
tTumRK4, yTum_VekRK4 = bym.add_coz_rk4_sistem(fonkVek, t0, tSon, yVek0, n)
# Çöz, scipy.integrate.odeint
tTumScipy= np.linspace(t0, tSon, n)
yTumScipy= spInt.odeint(fonkVek, yVek0, tTumScipy)
# Çiz
fig, ax = plt.subplots(2, 2)
# Sol üst grafik
ax[0,0].plot(tTumEuler, yTum_VekEuler[0], lw='5', label='x1(t) Euler', color='red')
ax[0,0].plot(tTumRK4, yTum_VekRK4[0], lw='3', label='x1(t) RK4', color='green')
ax[0,0].plot(tTumScipy, yTumScipy[:,0], lw='1', label='x1(t) Scipy', color='blue')
ax[0,0].set_xlabel('t')
ax[0,0].set_ylabel('x1(t)')
ax[0,0].legend()
# Sağ üst grafik
ax[0,1].plot(tTumEuler, yTum_VekEuler[1], lw='5', label='v1(t) Euler', color='blue')
ax[0,1].plot(tTumRK4, yTum_VekRK4[1], lw='3', label='v1(t) RK4', color='green')
ax[0,1].plot(tTumScipy, yTumScipy[:,1], lw='1', label='v1(t) Scipy', color='red')
ax[0,1].set_xlabel('t')
ax[0,1].set_ylabel('v1(t)')
ax[0,1].legend()
# Sol alt grafik
ax[1,0].plot(tTumEuler, yTum_VekEuler[2], lw='5', label='x2(t) Euler', color='red')
ax[1,0].plot(tTumRK4, yTum_VekRK4[2], lw='3', label='x2(t) RK4', color='green')
ax[1,0].plot(tTumScipy, yTumScipy[:,2], lw='1', label='x2(t) Scipy', color='blue')
ax[1,0].set_xlabel('t')
ax[1,0].set_ylabel('x2(t)')
ax[1,0].legend()
# Sağ alt grafik
ax[1,1].plot(tTumEuler, yTum_VekEuler[3], lw='5', label='v2(t) Euler', color='blue')
ax[1,1].plot(tTumRK4, yTum_VekRK4[3], lw='3', label='v2(t) RK4', color='green')
ax[1,1].plot(tTumScipy, yTumScipy[:,3], lw='1', label='v2(t) Scipy', color='red')
ax[1,1].set_xlabel('t')
ax[1,1].set_ylabel('v2(t)')
ax[1,1].legend()
plt.show()

Runge-Kutta-Fehlberg (RK45) Yöntemi

  • Runge-Kutta-Fehlberg (RK45) yöntemi, dördüncü ve beşinci mertebe Runge-Kutta yöntemi kullanarak çözüm elde eder.
  • RK45, değişken adım aralığı (adaptive step size) kullanır.
  • RK45’in girdilerinden (input) biri rölatif (relative) hata değeridir. RK4 ile elde edilen sonuç ile RK5 ile elde edilen sonuç arasındaki fark rölatif hatadan büyükse adım aralığı küçültülür ve o adım tekrar hesaplanır.
  • Bunun gibi metotlara değişken (adaptive) adım aralığı yöntemleri denir.
  • Runge-Kutta metotları açık (explicit) yöntemlerdir.

Scipy Paketleri Hakkında

  • Scipy içerisinde başlangıç değer problemi çözümleri için birçok fonksiyon var  [2].
  • Eski api olarak odeint hala kullanılabiliyor ancak scipy’ın gelecekteki versiyonlarında çıkarılma ihtimali çok yüksek.
  • odeint fonksiyonunda çağırılan \(f\) fonksiyonunda \(f(y,t)\) olmak zorunda. solve_ivp’de ise bu durum tam tersi. Bu farklılığı tfirst=True parametresi ile değiştirebilirsiniz - 1.1.0 versiyonundan itibaren.
  • odeint ile ilgili bilgi için  [3] referansına bakabilirsiniz.

Problemler

Problem 1

Aşağıdaki matris denklemlerini yazdığımız rk4 fonksiyonu ile çözün.

\[ \begin{align*} \frac{d}{dt} \begin{bmatrix} \rho_{11}(t) & \rho_{12}(t) \\ \rho_{21}(t) & \rho_{22}(t) \end{bmatrix} = \begin{bmatrix} -\rho_{11}(t) & -2\rho_{12}(t) \\ -3\rho_{21}(t) & -4\rho_{22}(t) \end{bmatrix} \end{align*} \]

Problem 2

scipy.integrate.solve_ivp fonksiyonunu kullanarak aşağıdaki denklemi çözün. Çözüm yöntemi olarak method='RK45' ve method='LSODA' kullanın.

\[ \frac{d}{dx}y(x) = \sin(5x) \]

Başlangıç koşulu: \(y(-3)=4.5\). Çözümü -3 ile 3 arasında çizdirin.

Aynı denklemi sağ taraf \(\sin(15x)\) ve \(\sin(25x)\) olacak şekilde çözün.

Problem 3

Sönümlü harmonik salınıcının Lagranjiyen’i aşağıdaki gibidir.

\[ L(x,v) = \frac{1}{2}mv^{2} - \frac{1}{2}kx^2 \]

Burada \(v\) hız olup konumun \(x\) zamana göre birinci türevine eşittir. Lagrange hareket denklemleri aşağıdaki gibi yazılır.

\[ \left(\frac{d}{dt}\frac{\partial L}{\partial v} \right) - \frac{\partial L}{\partial x} = \frac{\partial F^{dis}}{\partial v} \]

Yukarıda verilen formülden diferansiyel denklem setini oluşturun. Burada \(F^{dis}=-\frac{1}{2}bv^{2}\) sönümleyici kuvvettir. \(t=0-50\) arasındaki değerler için ve \(n=1000\) adımda, - Runge-Kutta-4 (veya 5(4) mertebe için Runge-Kutta methodu) yöntemi ile çözün. - x ile t’nin oluşturduğu grafiği çizdirin.

Başlangıç koşulları ve sabitler aşağıdaki gibidir.

  • \(m=1\) kg
  • \(k=1\) N/m
  • \(b=0.1\) Ns/m
  • \(x(t=0)=1\) m
  • \(v(t=0)=0\) m/s

Problem 4

Aşağıdaki diferansiyel denklemi Runge Kutta 4(5) yöntemini kullanarak \(N=100\) adımda çözün. Analitik sonuçla karşılaştırın.

\[ y''(t) = t^{3}+t+5, \quad y(0)=1, \quad y'(0)=5 \]

Problem 5

Düşey düzlemde \(l=1\) m uzunluğunda bir ipin ucunda birim kütleli bir cisim salınsın. Bu cisim denge noktasından \(\theta(t=0s)=\pi/9\) rad açı ile harekete başlıyor. \(g=9.81\) m/\(s^{2}\). Cisim bırakıldığı anda açısal hızı \(\omega(t=0s)=0\) rad/s. Bu cismin hareket denklemini \(t=(0,10)\) arasında

  1. Euler yöntemiyle çözün.
  2. RK4 (veya RK45 veya LSODA) yöntemiyle çözün.
  3. Analitik çözümü de dikkate alın ve tüm çözümleri üst üste bindirerek çizdirin.

Euler metodunda adım sayısını \(n=10000\) alın.

\[ \frac{d^{2}\theta}{dt^{2}}=-\frac{g}{l}\sin \theta \]

Analitik çözüm: \(\theta(t) = A\cos(t\sqrt{g/l})+B\sin(t\sqrt{g/l})\)