import numpy as np
from scipy.linalg import solve
import matplotlib.pyplot as plt
# Katsayı Matrisi
coefMat= np.array([\
[-2, 1, 0, 0, 0]\
, [ 1,-2, 1, 0, 0]\
, [ 0, 1,-2, 1, 0]\
, [ 0, 0, 1,-2, 1]\
, [ 0, 0, 0, 1,-2]])
# Sonuç Matrisi
solMat= np.array([-0.875,0.25,0.375,0.5,-8.375])
# Doğrusal Denklemi Çöz
solution=solve(coefMat, solMat)
# Başlangıç Değeri
y0= 1.
# Çözümü Başına Ekle
solution=np.insert(solution, 0, y0)
# Bitiş Değeri
ySon= 9.
# Çözümü Sonuna Ekle
solution=np.append(solution, ySon)
# Çözmek istediğimiz t değerleri
t= np.linspace(0, 3, 7)
# Çiz
plt.plot(t, solution, 'o-', label='Sonlu Farklar')
# --------------------------------------------
# Analitik Çözüm
# y''= t y(0)=1 y(3)=9
# y'= t^2/2 + C1
# y= t^3/6 + C1 t + C2
# y(0)=1 => C2=1
# y(3)=9 => 9=27/6 + 3C1 + 1=9 => C1=27/18
#=> y(t)= t^3/6 + 21t/18 + 1
# --------------------------------------------
plt.plot(t, t**3/6 + 21*t/18 + 1, 'r-', label='Analitik')
plt.legend()
plt.show()
SDP - Sonlu Farklar Yöntemi
Sonlu Farklar Yöntemi
Sonlu farklar yöntemini kullanarak sınır değer problemine ait diferansiyel denklem çözeceğiz. Bunun için Taylor açılımını kullanacağız. Taylor açılımında birinci türevi orta noktayı kullanarak yazalım.
\[ y'_{i} = \frac{y_{i+1}-y_{i-1}}{2h} \]
Farz edelim ki bu fonksiyonu \(g_{i}\) için yazmışız.
\[ g'_{i} = \frac{g_{i+1}-g_{i-1}}{2h} \]
Bu \(g_{i}\) fonksiyonu \(y'_{i}\) fonksiyonuna eşit olsun. Yani \(g_{i}=y'_{i}\) olur. Bu durumda
\[ \begin{align*} y''_{i} &= \frac{y'_{i+1}-y'_{i-1}}{2h} \\ &= \frac{1}{2h}\left(\frac{y_{i+2}-y_{i}}{2h}-\frac{y_{i}-y_{i-2}}{2h}\right) \end{align*} \]
Biraz düzenleyelim.
\[ \begin{align*} y''_{i} &= \frac{1}{4h^2}\left(y_{i+2}-y_{i}-y_{i}+y_{i-2}\right) \\ &= \frac{y_{i+2}-2y_{i}+y_{i-2}}{4h^2} \end{align*} \]
Eğer adım aralığını yarıya indirirsek yani \(h \rightarrow h/2\) yaparsak \(y_{i+2}\) noktasından \(y_{i+1}\) noktasına gelebiliriz. Benzer şekilde \(y_{i-2}\) noktasından \(y_{i-1}\) noktasına geliriz. O halde ikinci türeve ait olan denklem şöyle olur.
Bu işlem sadece sabit adım aralığı için geçerlidir.
\[ y''_{i} = \frac{y_{i+1}-2y_{i}+y_{i-1}}{h^2} \]
Bu denklemi kullanarak sınır değer problemlerini çözeceğiz. Yukarıdaki bağıntı, diferansiyel denklemi cebirsel denklem sistemine dönüştürür.
\[ y''_{i} = f\left(x,y_{i},\frac{y_{i+1}-y_{i-1}}{2h}\right)=\frac{y_{i+1}-2y_{i}+y_{i-1}}{h^2} \]
Alıştırma 1
Aşağıdaki denklemi adım adım çözelim.
\[ y''(t)=t, \quad y(0)=1, \quad y(3)=9 \]
Çözüm (Elle)
Adım aralığı \(h=1\) olsun. Çözeceğimiz \(t\)’nin değerleri \(0, 1, 2, 3\) olacaktır. Bu değerlerden \(y\) değerlerini bulalım.
\[ y''(t)=t=\frac{y_{i+1}-2y_{i}+y_{i-1}}{h^2} \]
Çözüme \(y(0)\)’dan başlayalım. Bunun için \(y(0) \equiv y_{i-1}\) almamız gerekecektir. Aksi taktirde köşe noktalar için denklemi çözemeyiz.
Aşağıdaki tabloyu \(h=1\) için yazalım.
\(y_{0}\) | \(y_{1}\) | \(y_{2}\) | \(y_{3}\) |
---|---|---|---|
\(y(0)\) | \(y(0+h)\) | \(y(0+2h)\) | \(y(0+3h)\) |
\(y(t=0)\) | \(y(t=1)\) | \(y(t=2)\) | \(y(t=3)\) |
1 | ? | ? | 9 |
Örneğin aşağıdaki tabloyu da \(h=0.5\) için yazalım.
\(y_{0}\) | \(y_{1}\) | \(y_{2}\) | \(y_{3}\) | \(y_{4}\) | \(y_{5}\) | \(y_{6}\) |
---|---|---|---|---|---|---|
\(y(0)\) | \(y(0.5)\) | \(y(1)\) | \(y(1.5)\) | \(y(2)\) | \(y(2.5)\) | \(y(3)\) |
1 | ? | ? | ? | ? | ? | 9 |
Şimdi yukarıdaki elde ettiğimiz denklemi yazalım.
\[ y''(t)=t=\frac{y_{i+1}-2y_{i}+y_{i-1}}{h^2} \]
\(i=1\) olmalı. Bundan dolayı \(t\) ilk değil bir sonraki adımından yani \(t+h=0+0.5\)’den itibaren başlatmalıyız. \[ 0.5 = \frac{y_{2}-2y_{1}+y_{0}}{h^2} \]
Bu denklemi \(h=0.5\) için adım adım yazalım.
\[ \begin{align*} 0.5\times (0.5)^{2} =& y_{2}-2y_{1}+y_{0}\\ 1\times (0.5)^{2} =& y_{3}-2y_{2}+y_{1}\\ 1.5\times (0.5)^{2} =& y_{4}-2y_{3}+y_{2}\\ 2\times (0.5)^{2} =& y_{5}-2y_{4}+y_{3}\\ 2.5\times (0.5)^{2} =& y_{6}-2y_{5}+y_{4} \end{align*} \]
Yukarıdaki denklemleri matris şeklinde yazalım. Yazarken \(y_{0}=1\) ve \(y_{6}=9\) yazacağız.
\[ \begin{align*} -0.875 =& y_{2}-2y_{1}\\ 0.25 =& y_{3}-2y_{2}+y_{1}\\ 0.375 =& y_{4}-2y_{3}+y_{2}\\ 0.5 =& y_{5}-2y_{4}+y_{3}\\ -8.375 =& -2y_{5}+y_{4} \end{align*} \]
Matris şeklinde yazalım.
\[ \begin{bmatrix} -2 & 1 & 0 & 0 & 0\\ 1 & -2 & 1 & 0 & 0\\ 0 & 1 & -2 & 1 & 0\\ 0 & 0 & 1 & -2 & 1\\ 0 & 0 & 0 & 1 & -2\\ \end{bmatrix} \begin{bmatrix} y_{1}\\ y_{2}\\ y_{3}\\ y_{4}\\ y_{5}\\ \end{bmatrix}= \begin{bmatrix} -0.875\\ 0.25\\ 0.375\\ 0.5\\ -8.375 \end{bmatrix} \]
Bu denklemi çözmeyi biliyoruz. Denklemi çözdüğümüzde \(y_{1}, y_{2}, y_{3}, y_{4}, y_{5}\) değerlerini bulmuş oluruz.
Çözüm (Kod)
Matris denklemini el ile yazmadan yani başlangıçta biz vermeden çözmeye çalışalım.
import numpy as np
from scipy.linalg import solve
import matplotlib.pyplot as plt
# Fonksiyon
def fonk(t):
return t
# Sınır Değerleri
y_t0=1
y_t3=9
# Adım aralığı
h=0.5
# Tum t değerleri
t0=0
tSon=3
# adet=(tSon-t0/h)+1
adet= int((tSon-t0)/h)+1
tTum=np.linspace(t0, tSon, adet)
# Yerine koyduğumuzda oluşan denklemler (tabloya bakınız)
# y''(t)=t = (y_{i+1}-2y_{i}+y_{i-1})/h^2
# t * h^2 = y_{i+1}-2y_{i}+y_{i-1}
# ilk denklem için t= h olacak. i=1 olacak
# h * h^2 = y_{2}-2y_{1}+y_{0}
# ilkDenklemSagTaraf= h*(h**2)-y_{0}
# Son denklem için t= tSon-h olacak. i=son-1 olacak
# (tSon-h) * h^2 = y_{son}-2y_{son-1}+y_{son-2}
# sonDenklemSagTaraf=(tSon-h)*(h**2)-y_{son}
ilkDenklemSagTaraf=fonk(h)*(h**2)-y_t0
sonDenklemSagTaraf=fonk(tSon-h)*(h**2)-y_t3
# Ara denklemleri oluştur.
# araDenklemlerSagTaraf[0]=ilkDenklemSagTaraf
araDenklemlerSagTaraf=np.array(ilkDenklemSagTaraf)
for i in range(1, adet-3):
araDenklemlerSagTaraf= np.append(araDenklemlerSagTaraf, fonk(tTum[i+1])*(h**2))
# Son denklemi ekle
sonucMatrisi= np.append(araDenklemlerSagTaraf,sonDenklemSagTaraf)
# Katsayı Matrisi
katsayiMat= np.zeros((adet-2, adet-2))
# y''(t)=t = (y_{i+1}-2y_{i}+y_{i-1})/h^2
# Denkleminde katsayılar hep aynı olur.
for i in range(adet-2):
katsayiMat[i, i]= -2
if not i == adet-3:
katsayiMat[i, i+1]= 1
if not i == 0:
katsayiMat[i, i-1]= 1
# Çöz
print(f"Katsayı Matrisi: \n{katsayiMat}")
print(f"Sonuç Matrisi: \n{sonucMatrisi}")
cozum=solve(katsayiMat, sonucMatrisi)
# İlk Çözümü Ekle
cozum=np.insert(cozum, 0, y_t0)
# Son Çözümü Ekle
cozum=np.append(cozum, y_t3)
# Çiz
plt.plot(tTum, cozum, 'o-', label='Sonlu Farklar')
plt.plot(tTum, tTum**3/6 + 21*tTum/18 + 1, 'r-', label='Analitik')
plt.legend()
plt.show()
Katsayı Matrisi:
[[-2. 1. 0. 0. 0.]
[ 1. -2. 1. 0. 0.]
[ 0. 1. -2. 1. 0.]
[ 0. 0. 1. -2. 1.]
[ 0. 0. 0. 1. -2.]]
Sonuç Matrisi:
[-0.875 0.25 0.375 0.5 -8.375]
Sonlu Farklar Yönteminin Çalışmadığı Örnek
Eğer oluşturduğumuz diferansiyel denklem sistemi doğrusal değilse, yukarıdaki algoritmayı kullanmayız. Aşağıdaki örneği inceleyelim.
\[ y''(t)=tyy'(t) \]
Çözmemiz gereken denklem
\[ y''(t)=tyy'(t)=\frac{y_{i+1}-2y_{i}+y_{i-1}}{h^2} \]
\(i=1\) için yazalım.
\[ ty_{1}y_{1}'=\frac{y_{2}-2y_{1}+y_{0}}{h^2} \]
Türevi de açıkça yazalım.
\[ ty_{1} \frac{y_{2}-y_{0}}{2h}=\frac{y_{2}-2y_{1}+y_{0}}{h^2} \]
Düzenleyelim.
\[ \frac{th}{2} (y_{1}y_{2}-y_{1}y_{0}) -y_{2}+2y_{1}-y_{0} = 0 \]
\(y_{1}y_{2}\) gibi terimlerin varlığından dolayı bu denklem doğrusal olmayan (nonlineer) denklemdir. Klasik yöntemlerle (LU ayrıştırma veya Gauss eleme) çözülemez.
Problemler
Problem 1
Aşağıdaki diferansiyel denklemi sonlu farklar yöntemi kullanarak çözebilir misinz?
\[ \frac{\hbar^{2}}{2 m}\frac{d^{2}}{dx^{2}}\psi(x) = E\psi(x), \qquad \psi(x=0)=0, \qquad \psi(x=L)=0 \text{.} \]
Çözümü \(n=1\) ve \(n=2\) için bulmaya çalışın. Hem bulduğunuz enerji değerlerini hem de dalga fonksiyonunun karesini (olasılık yoğunluğunu, \(|\psi(x)|^{2}\)) çizdirin. \(x\) değerleri \(0\) ile \(L\) arasında \(N=100\) adet olsun.
Denklemdeki katsayıları aşağıdaki gibi alın.
- \(\hbar=1\)
- \(m=1\) kg
- \(L=2\)
Enerji özdeğerinin gerçek değeri aşağıdaki gibi olmalıdır.
\(E=\frac{\pi^{2}\hbar^{2}}{2mL^{2}}n^{2}\) J.
Dalga fonksiyonunun analitik çözümü aşağıdaki gibi olmalıdır.
\[ \psi(x)=\sqrt{\frac{2}{L}}\sin\left(\frac{n\pi x}{L}\right) \]
Bu problem bir özdeğer problemidir. Diferansiyel operatörü \(\left(\frac{d^{2}}{dx^{2}}\right)\) doğrusal operatörü matris şeklinde yazmalısınız. Bu matrisi sonlu farklar yönteminde yazdık.
Problem 2
Aşağıdaki diferansiyel denklemi sonlu farklar yöntemini kullanarak çözün. Adım aralığını \(h=0.05\) alın. Analitik sonuçla karşılaştırın.
\[ y''(t) = t^{3}+t+5, \quad y(0)=1, \quad y(6)=545.8 \]
Kaynaklar
- Numerical Methods in Engineering with Python 3, Jaan Kiusalaas, Cambridge University Press, 2013, Syf:307
- Sonlu Farklar Yöntemi, Youtube Videosu
- Numerical Solutions of Schrodinger’s Equation, Neill Lambert, 2001