Transcript Bölüm 7 Adi Diferansiyel Denklemler
BÖLÜM 7 ADİ DİFERANSİYEL DENKLEMLER
1
BÖLÜM 7.1
Giriş
2
7.1 Giriş
İlk derslerde Newton’un 2. yasasına dayanılarak türetilen ve düşen bir paraşütçünün hızının zamana göre değişimini ifade eden aşağıdaki bağıntıdan söz edilmişti;
dv g c dt
m v
(PT7.1)
Burada
v
paraşütçünün hızı (
bağımlı değişken
),
t
de zaman (
bağımsız değişken
) idi. Denklemdeki diğer büyüklükler ise:
g
yerçekimi ivmesi,
c
direnç katsayısı,
m
de düşen paraşütçünün toplam kütlesi ollarak belirtilmişti. Bu şekilde bilinmeyen bir fonksiyon ve onun türevleri cinsinden ifade edilen denklemlere
diferansiyel denklemler
denir. Diferansiyel denklemler bazen oran denklemleri olarak da adlandırılır. Çünkü, denklemden de görüldüğü üzere diferansiyel denklemler, değişkenlerin değişim oranlarını değişkenlerin kendilerinin ve çeşitli parametrelerin fonksiyonları cinsinden ifade eden bir yapıya sahiptirler. Bir çok fiziksel olay, en iyi, değişim oranları şeklinde matematiksel olarak formüle edilebildiğinden, bu tür denklemlere mühendisliğin her alanında sıklıkla rastlanır ve bu nedenle diferansiyel denklemler mühendislikte önemli bir yere sahiptirler.
3
7.1 Giriş
Eğer diferansiyel denklem
tek bir bağımsız değişken
diferansiyel denklemlere
adi diferansiyel denklemler
içeriyorsa, bu tür (
ADD
) denir. Eğer
birden fazla bağımsız değişken
söz konusu ise, bu tür diferansiyel denklemlere de
kısmi türevli diferansiyel denklemler
(
KDD
) denir.
Diferansiyel denklemler mertebelerine göre de sınıflandırılırlar. Diferansiyel denklemdeki en yüksek türev mertebesi, diferansiyel denklemin mertebesi olarak adlandırılır. Örneğin, bir kütle yay sönüm elemanı sistemindeki kütlenin konumunu ifade eden aşağıdaki diferansiyel denklem ikinci mertebeden adi diferansiyel denklemdir.
m d
2 2
x
c dx
kx
0 (PT7.2)
dt dt
Bu denklemde
m
kütle,
c
sönüm katsayısı, yeni bir y değişkeni tanımlanırsa;
k
da yay katsayısıdır.
n. mertebeden bir diferansiyel denklem, n tane 1. dereceden diferansiyel şeklinde ifade edilebilir. Örneğin yukarıdaki denklemde aşağıdaki gibi
y
dx
(PT7.3)
dt
4
7.1 Giriş
ve bu değişkenin t ye göre diferansiyeli alınırsa;
dy dt
d dt
2 2
x
(PT7.4)
elde edilir. (PT7.3) ve (PT7.4) ifadeleri (PT7.2) de yerine yazılırsa;
m dy dt
cy
kx
0 (PT7.5) veya
dy dt
cy
m kx
(PT7.6) elde edilir. (PT7.3) ve (PT7.5) (veya (PT7.6)) denklemleri, (PT7.2) denklemine eşdeğer iki adet birinci mertebeden diferansiyel denklemden oluşan bir diferansiyel denklem takımıdır.
5
7.1 Giriş
Bir diferansiyel denklemin çözümü, bağımsız değişken ve parametrelerin, orijinal diferansiyel denklemi sağlayan özel bir fonksiyonudur. Bu kavramı aşağıdaki fonksiyon üzerinde inceleyelim;
y
0 .
5
x
4
4
x
3
10
x
2
8 .
5
x
1 (PT7.12)
Şekilde grafiği görülen yukarıdaki fonksiyon dördüncü dereceden bir polinomdur.
FIGURE PT7.3
Plots of (
a
)
y
versus
x
and (
b
)
dy/dx
versus
x
for the function 6
7.1 Giriş
(PT7.12) eşitliğinin x ‘e göre diferansiyeli alınırsa;
dy dx
2
x
3 12
x
2 20
x
8 .
5 (PT7.13) elde edilir. Elde edilen bu eşitlik orjinal fonksiyonun davranışını farklı bir şekilde ifade etmektedir. Orijinal fonksiyon her bir x değerine karşılık gelen y değerini gösterirken, elde edilen eşitlik her bir x değerine karşılık, bu noktada orijinal fonksiyonun sahip olduğu eğimi göstermektedir. Yandaki şekilde (Şekil PT7.3) orjinal fonksiyonun ve diferansiyelinin (eğiminin) x ‘e göre değişimi görülmektedir.
FIGURE PT7.3:
Plots of (
a
)
y
versus
x
and (
b
)
dy/dx
versus
x
for the function 7
7.1 Giriş
Dikkat edilirse, orijinal fonksiyonun yatay olduğu (eğimin sıfır olduğu) noktalarda diferansiyel sıfır değeri almaktadır. Yine diferansiyelin ekstremum yaptığı noktalarda orijinal fonksiyon mutlak değerce en büyük eğimlere sahip olmaktadır.
Yukarıda bir fonksiyondan hareket ederek buna karşılık gelen diferansiyel denklemin elde edilmesi incelendi. Halbuki bu konunun amacı tam tersine olarak, verilen bir diferansiyel denklemden hareketle, bu diferansiyel denkleme karşılık gelen orijinal fonksiyonun belirlenmesidir. Belirlenmek istenen orijinal fonksiyona diferansiyel denklemin çözümü adı verilir. FIGURE PT7.3
8
7.1 Giriş
Şimdi yukarıda elde ettiğimiz diferansiyel denklemden hareketle orijinal denklemi belirlemeye çalışalım. Bunun için (PT7.13) eşitliğinin integralini almamız gerekir;
dy
2
x
3
12
x
2
20
x
8 .
5
dx
İntegrasyon gerçekleştirilirse;
y
2
x
4
4
12
x
3
3
20
x
2
2
8 .
5
x
C y
0 .
5
x
4
4
x
3
10
x
2
8 .
5
x
C
(PT7.14)
elde edilir. Elde edilen çözüm, orijinal denklemle bir fark dışında özdeştir. Orijinal denklemin diferansiyeli alınırken yok olan
1
sabiti yerine integrasyon sonunda
sabiti C
sabiti gelmektedir.
C
sabitine
integrasyon
adı verilir. Böyle bir sabitin ortaya çıkması çözümün tek olmaması anlamına gelmektedir. Diğer bir ifadeyle sonsuz sayıda C değerine karşılık gelecek sonsuz sayıda çözüm söz konusudur. 9
7.1 Giriş
Yandaki şekilde diferansiyel denklemi sağlayan olası altı çözüm görülmektedir.
Dolayısıyla tek bir çözüm elde edilmek isteniyorsa, diferansiyel denklemle birlikte bazı ilave şartların da bilinmesi gerekir. Birinci mertebeden bir diferansiyel denklem için bu şart da bir tanedir ve
başlangıç şartı
(ya da
değeri
) adı verilir. Örneğin x=0 için y=1 şeklindeki bir başlangıç şartı bilinirse, diferansiyel denklemi sağlayan tek çözümü belirlemek mümkün olur.
FIGURE PT7.4
Six possible solutions for the integral of
dy
. Each conforms to a different value of the constant of integration
C
.
10
7.1 Giriş
Bu değerler (PT7.14) denkleminde yerine yazılırsa; 1 0 .
5 4 4 3 10 2 8 .
5
C
(PT7.15) elde edilir. Buradan C =1 olarak belirlenir. Belirlenen değer (PT7.14)’de yerine yazılırsa, diferansiyel denklemi ve başlangıç şartını sağlayan tek çözüm;
y
0 .
5
x
4
4
x
3
10
x
2
8 .
5
x
1 (PT7.16)
şeklinde belirlenmiş olur. Başlangıç şartının bilinmesiyle çözüm fonksiyonu sanki başlangıç noktasından sabitlenerek tek bir çözümün elde edilmesini sağlamaktadır.
11
7.1 Giriş
Fiziksel problemlere karşılık gelen diferansiyel denklemler için başlangıç değerleri önemli bir anlama sahiptir. Örneğin düşen paraşütçü probleminde başlangıç şartının sıfır olması, paraşütçünün düşey doğrultudaki başlangıç hızının sıfır olması anlamına gelmektedir.
n. mertebeden bir diferansiyel denklem n tane birinci mertebeden diferansiyel denklemden oluşan denklem takımına karşılık geldiğine göre, n. mertebeden diferansiyel denklemin çözümü için n tane şarta ihtiyaç duyulacağı açıktır. Eğer bu şartların tümü başlangıç değeri için veriliyorsa, bu tür problemlere
başlangıç değer problemleri
, diğer noktalar için de verilmiş ise
sınır değer problemleri
adını alır.
12
7.2 Euler Yöntemi
Bu bölümde;
dy
f
(
x
,
y
)
dx
yapısındaki
adi diferansiyel denklemlerin
sayısal yöntemlerle çözümü incelenecektir. Bu yöntemler
New value
old value
slope
x
step size
şeklinde formüle edilebilecek genel bir yapıya sahiptir. Bu eşitliğe göre eğimi ve y
i
eski değerinden yararlanılarak, h adımı sonundaki y
i+1
yeni değeri belirlenebilir (Şekil 25.1);
y i
1
y i
h
(25.1)
FIGURE 25.1
Graphical depiction of a onestep method.
13
7.2 Euler Yöntemi
Bu formül adım adım uygulanarak belirli bir aralıktaki çözüm adım adım belirlenmiş olur. Bütün tek adımlı yöntemlerde aynı temel düşünce kullanılmaktadır. Farklı yöntemler, eğimin belirlenmesindeki farklılıklardan kaynaklanmaktadır.
Eğimin belirlenmesinde en basit yol, adımın başındaki eğimi doğrudan kullanmaktır. Adımın başındaki eğimin adım boyunca geçerli olduğu varsayımına dayalı olan bu yaklaşım
Euler yöntemi
olarak adlandırılır.
FIGURE 25.1
Graphical depiction of a onestep method.
14
7.2 Euler Yöntemi
Birinci türev x oluşturur;
i
noktasında fonksiyonun eğimi için doğrudan bir tahmin
f
(
x i
,
y i
)
Burada f (x
i ,y i
), diferansiyel denklemin x
i
, y
i
noktasındaki değerinden başka bir şey değildir. Bu eğim değeri (25.1) eşitliğinde yerine yazılırsa;
FIGURE 25.2
Euler’s method.
y i
1
y i
f
(
x i
,
y i
)
h
(25.2)
Bu formül,
Euler
(
Euler
-
Cauchy
olarak adlandırılır.
veya
nokta
-
eğim
)
yöntemi
15
7.2 Euler Yöntemi
Example 25.1
h=0.5 adımı ile, x=0 da y=1 başlangıç koşulu altında Euler yöntemiyle çözünüz. Aşağıdaki diferansiyel denklemi x=0 ile x=4 aralığında,
dy dx
2
x
3 12
x
2 20
x
8 .
5 Not: Analitik çözüm
y
0 .
5
x
4
4
x
3
10
Çözüm:
x
0
x
0
,
y
1
x y
2
(
8 .
5
x
1 (PT7.16) 0 )
y
(
x
0
)
y
0
y i
1
y i
f
(
x i
,
y i
)
h y
1
y
0
f
(
x
0
,
y
0
)
h y
( 0 .
5 )
y
( 0 )
f
( 0 , 1 ) * 0 .
5 16
7.2 Euler Yöntemi
x=0 ve y(0)=1 için (0,1) noktasında çözüm fonksiyonunun eğimi;
f
( 0 , 1 ) 2 3 12 2 20 8 .
5 8 .
5 şeklinde belirlenir. Buna göre adımın sonundaki y değeri, Euler yönteminden;
y
( 0 .
5 ) 1 8 .
5 * 0 .
5 5 .
25 bulunur. Halbuki x=0.5 deki analitik çözüm ise
y
0 .
5 4 4 3 10 2 8 .
5 1 3 .
21875 şeklindedir. İki çözüm arasındaki hata oranı;
t
3.21875
100
2.03125
100 3.21875
t
63 .
1 %
kadar olmaktadır.
17
7.2 Euler Yöntemi
İkinci adım (2h=2*0,5=1) için çözüm;
y
( 1 )
y
( 0 .
5 )
f
( 0 .
5 , 5 .
25 ) * 0 .
5 5 .
25 [ 2 3 12 2 20 8 .
5 ] * 0 .
5
5 .
875
şeklinde elde edilir.
x=1.0 deki analitik çözüm 3.0 olup, hata oranı;
t
95 .
8 %
olmaktadır. 18
25.1 EULER YÖNTEMİ
Tablo 25.1 başlangıç şartı x=0’da y=1 için integralinin kesin ve yaklaşık değerlerinin karşılaştırılması. Yaklaşık değerler, 0.5’lik bir adım büyüklüğünde Euler yöntemi kullanılarak hesaplanmıştır. Yerel hatalar, tek bir adım boyunca oluşan hataya göre olup Taylor serisi açılımı ile hesaplanmıştır. Genel hata, önceki ve şimdiki adımlar arasındaki farktır.
19
7.2 Euler Yöntemi
Dikkat edileceği üzere sayısal çözüm analitik çözümün genel karakteristiğini yakalamakta fakat hatanın mertebesi oldukça büyük olmaktadır. Hatanın büyüklüğü adım miktarının büyük olmasından kaynaklanmaktadır. Adım miktarı küçültüldükçe hata mertebesi de azalacaktır.
FIGURE 25.3
Comparison of the true solution with a numerical solution using Euler’s method for the integral of
dy
from
x
0 to
x
4 with a step size of 0.5. The initial condition at
x
0 is
y
1.
20
7.2 Euler Yöntemi
Örneğin adım miktarı yarıya düşürüldüğünde aşağıdaki çözüm elde edilmektedir.
FIGURE 25.4
(
a
) Comparison of two numerical solutions with Euler’s method using step sizes of 0.5 and 0.25.
21
7.3 Heun Yöntemi
Euler yönteminde hatanın temel bir kaynağı aralığın başlangıcındaki eğimin, tüm aralık boyunca geçerli olduğunun varsayılmasıdır. Bu yetersizliğin olumsuz etkisini azaltmak için iki basit düzeltme yapılabilir. Bunlardan biri, aralığın sadece başlangıcındaki eğimi kullanmak yerine, biri aralığın başlangıcındaki, diğeri de aralığın sonundaki olmak üzere iki noktadaki eğimin ortalamasının aralık boyunca geçerli olduğunun varsayılmasıdır. Bu yaklaşım Heun yöntemi olarak adlandırılmaktadır (Şek. 25.9). Euler yönteminde aralığın başındaki eğim;
y i
f
(
x i
,
y i
) (25.12)
şeklinde belirlenmekteydi. Bu değer Heun yönteminde y
i+1
lineer ektrapolasyonunda kullanılmaktadır; değerinin 22
7.3 Heun Yöntemi
y i
0 1
y i
f
(
x i
,
y i
)
h
(25.13) Standart Euler yönteminden farklı olarak (25.13) eşitliği ile belirlenen
y 0 i+1
değeri son değer değil, bir ara tahmin değeridir ve bu nedenle ‘0’ üst indisi ile gösterilmiştir. (25.13) eşitliği bu nedenle “
tahmin denklemi
(
predictor equation
)” olarak adlandırılır ve aralığın sonundaki eğim değerinin;
y
i
1
f
(
x i
1 ,
y i
0 1 ) (25.14) şeklinde hesaplanabilmesi için bir tahmin değeri oluşturmaktadır.
FIGURE 25.9
Heun yönteminin grafik olarak açıklanması. (
a
) Deneme (
b
) düzeltme.
23
7.3 Heun Yöntemi
y
Böylece (25.12) ve (25.14) eğimleri, aralık için ortalama eğimin belirlenmesinde kullanılabilir;
y i
2
y
i
1
f
(
x i
,
y i
)
2
f
(
x i
1
,
y i
0 1
)
Bu ortalama eğim artık y yararlanılarak lineer bir ekstrapolasyonla
y i+1 i
den in belirlenmesinde kullanılabilir;
y i
1
y i
y
h
y i
y i
2
y i
1
h
y i
f
(
x i
,
y i
)
2
f
(
x i
1
,
y i
0 1
)
h
Bu son denklem “
düzeltme denklemi
adlandırılır. (
corrector equation
)” olarak 24
7.3 Heun Yöntemi
Böylece Heun yöntemi bir “
tahmin etme-düzeltme
(
predictor-corrector approach
)” yöntemi olarak tanımlanabilir; Tahmin etme (Şek. 25.9a):
y i
0 1
y i
f
(
x i
,
y i
)
h
(25.15) Düzeltme (Şek. 25.9b);
y i
1
y i
f
(
x i
,
y i
) 2
f
(
x i
1 ,
y i
0 1 )
h
(25.16) Dikkat edileceği gibi (25.16) eşitliğinin her iki tarafında da y kullanılabilir.
i+1
bulunduğundan, bu eşitlik iteratif bir yapıdadır ve eski değerler kullanılarak daha iyileştirilmiş yeni değerlerin belirlenmesinde
Denklem
(25.15)
Denklem
(25.16) 25
7.3 Heun Yöntemi
Tabii ki bu iteratif yaklaşımın analitik çözümün bulunmasını sağlamasını bekleyemeyiz. İyileştirme ancak yaklaşımın sahip olduğu belirli bir kesme hatası ile meydana gelecektir. Bu durum aşağıda bir örnekle incelenecektir.
İterasyonun sonlandırılmasında daha önceden açıklandığı üzere aşağıdaki yaklaşık hata oranı hesabı kullanılacaktır;
a
y i j
1
y i j
1
y i j
1 1
100 % (25.17)
burada y
j-1 i+1
ve y
j i+1
değerleri sırasıyla bir önceki ve son iterasyon adımında belirlenen değerlerdir.
26
7.3 Heun Yöntemi
Example 25.5
y’ =4 e 0.8x
– 0.5y diferansiyel denklemini x=0 dan x=4 e kadar h=1 adım miktarı ile integre ediniz. Başlangıç şartı: x=0 da y=2.
Çözüm:
Sayısal çözümü gerçekleştirmeden önce analitik çözümün
y
4 1 .
3
e
0 .
8
x
e
0 .
5
x
2
e
0 .
5
x
(E25.5.1) şeklinde olduğunu hatırlatalım ve bu çözümü Tablo 25.2 deki analitik değerlerin hesaplanmasında kullanıldığını belirtelim.
İlk olarak (x
0
,y
0
) daki eğim değeri bulunabilir;
y
0
f
(
x
0 ,
y
0 ) 4
e
0 .
8 *( 0 ) 0 .
5 * ( 2 ) 3 Şimdi tahmin denklemini [(25.15) eşitliği] kullanılarak y’nin 1’deki ilk tahmini çözümü belirleyebiliriz;
y
1 0
y
0
f
(
x
0 ,
y
0 )
h
2 ( 3) * 1 5 27
7.3 Heun Yöntemi
Belirlenen bu değer normal olarak standart Euler yöntemiyle belirlenecek olan değerdir. Tablo 25.2 deki analitik değer ile karşılaştırılırsa bulunan yaklaşık çözümün %25.3 gerçek hata oranına sahip olduğunu görebiliriz.
Tablo 25.2
x y true
0 2.0000000
1 6.1946314
2 14.8439219
3 33.6771718
4 75.3389626
1
y Heun
2.0000000
6.7010819
16.3197819
37.1992489
83.3377674
İterasyonlar
t
(%) 0.00
8.18
9.94
10.46
10.62
15 ( bu sütunu dikkate almayınız )
y Heun
2.0000000
6.3608655
t
(%) 0.00
2.68
15.3022367
34.7432761
77.7350962
2.68
3.17
3.18
28
7.3 Heun Yöntemi
Şimdi y
i+1
deki çözümü iyileştirmek için y
0 1
eğimin belirlenmesinde kullanılabilir; değeri aralığın sonundaki
y
1
f
(
x
1
,
y
1 0
)
4
e
0 .
8 *( 1 )
0 .
5 * ( 5 )
6 .
402164
Belirlenen bu değer ile aralığın başındaki eğim kullanılarak x=0 dan 1 e kadar olan aralık boyunca ortalama eğim belirlenebilir;
y
y
i
y
i
1 2 3 6 .
402164 2 4 .
701082 Bu eğim değeri, gerçek ortalama eğim değeri olan 4.1946 a daha yakın bir değerdir. Şimdi bu eğim değeri düzeltme denkleminde [(25.16) eşitliği] yerine yazılarak x=1 deki iyileştirilmiş yaklaşık çözüm belirlenebilir;
y
1
y
0
y
h
2
( 4 .
701082 ) * 1
6 .
701082
Bu yaklaşık çözüm %8.18 gerçek hata oranına sahip olup, iterasyon yapılmaksızın belirlenen bu değer Euler yöntemi ile elde edilenden daha yüksek bir doğruluğa sahiptir.
29
7.3 Heun Yöntemi
Şimdi bulunan değerin (25.16) eşitliğinin sağ tarafına yeniden konularak hesabın tekrarlanmasıyla y
1
çözümü iteratif olarak iyileştirilebilir; ) , )
y i
1
i
( ,
i i
2
i
1
y i
0 1
h
e
0.8*1 0.5*6.701082)]
y
1 2 Elde edilen son çözüm % 1.31 yaklaşık hata oranına sahiptir. Bu çözüm kullanılarak iterasyona devam edilirse;
e
0.8*1 0.5*6.275811)]
y
1 2 Bu kez yaklaşık hata oranı %3.03 olmaktadır. Bu da iteratif yaklaşımlarda bazen hatanın büyüyebileceğini bize göstermektedir.
Not:
yukarıdaki hesaplamalar iteratif hesaplamalardır. Aşağıdaki tabloda Heun yöntemine göre hesaplanan 2,3 ve 4. değerleri hesaplayınız. Bu ödevdir ve sorumlusunuz.
30
7.3 Heun Yöntemi
Bu tür hata artışları özellikle büyük adım değerlerin kullanılması durumlarında meydana gelebilir ve bu gibi durumlar, ard arda gelen iterasyonlar sonucunda iteratif yaklaşımlarda mutlaka yakınsamanın meydana geleceği şeklinde genel bir sonuca ulaşmamızı engeller. Bununla birlikte yeterince küçük adımlar kullanıldığında, ard arda gelen iterasyonlar sonucu çözüm belirli bir değere yakınsayabilir. Yukarıdaki örnekte 15 iterasyon sonucunda 6.360865 değerine %2.68 yaklaşık hata oranı ile ulaşılmaktadır. Tablo 25.2 de her bir adım için tek iterasyon ve 15 iterasyon sonucu elde edilen değerlerin karşılaştırılması;
Tablo 25.2
x
0 1 2 3 4
y true
2.0000000
6.1946314
14.8439219
33.6771718
75.3389626
1
y Heun
2.0000000
6.7010819
16.3197819
37.1992489
83.3377674
İterasyonlar
t
(%) 0.00
8.18
9.94
10.46
10.62
15 ( bu sütunu dikkate almayınız )
y Heun
2.0000000
6.3608655
15.3022367
34.7432761
77.7350962
t
(%) 0.00
2.68
2.68
3.17
3.18
31
Runga-Kutta Yöntemleri
7.4 Orta-Nokta Runga-Kutta Yöntemi
Euler yönteminin iyileştirilmesi için yapılabilecek bir diğer iyileştirme Orta-nokta R-K yöntemi olarak adlandırılan, önce aralığın orta noktasının belirlenmesi için Euler yöntemi kullanılmakta, sonra da bu noktanın eğiminin aralık boyunca geçerli olduğu varsayımıyla aralık sonundaki değer belirlenmektedir.
y i
1 / 2
y
i
1 / 2
y i f
f
(
x i
1 / (
x i
, 2
,
y i
)
h
2 (25.25)
y i
1 / 2
) (25.26)
y i
1
y i
f
(
x i
1 / 2
,
y i
1 / 2
)
h
(25.27)
FIGURE 25.12 : Orta nokta yönteminin grafiksel açıklanması.
(
a
) Eq. (25.25) and (
b
) Eq. (25.27).
32
7.4 Orta-Nokta Runga-Kutta Yöntemi
y i
1
y i
f
(
x i
1 / 2
,
y i
1 / 2
)
h
(25.27)
Eğim olarak aralığın orta noktasının eğimi kullanıldığından, yöntemin Euler yöntemine göre daha iyi sonuçlar vereceği tahmin edilebilir.
FIGURE 25.12
Graphical depiction of the midpoint method.
(
a
) Eq. (25.25) and (
b
) Eq. (25.27).
33
7.4.
4.Mertebeden Klasik Runga-Kutta Yöntemi
Çeşitli farklı Runge-Kutta yöntemleri vardır. Ancak bunların tümü aşağıdaki gibi ortak bir genel formda yazılabilir;
y i
1
i
( ,
i i
burada
(x
i
,y
olarak düşünülebilir. Artım fonksiyonu aşağıdaki genel yapıda yazılabilir;
i
,h)
artım fonksiyonu
olarak adlandırılır ve
temsili eğim
( ,
i i
a k
1 1
a k
2 2
a k n n
(25.29)
burada a ‘lar katsayılar olup, k ‘lar da çözüm fonksiyonunun çeşitli noktalarındaki eğimleri olarak düşünülebilir. Bunların belirlenmesi konusuna burada girilmeyecektir. Ancak uygulamada en çok kullanılan R-K yöntemi olan 4. mertebeden R-K yöntemi için katsayıların verilmesi ile yetinilecektir.
34
7.4.
4.Mertebeden Klasik Runga-Kutta Yöntemi
En çok kullanılan R-K yöntemi 4. derecedendir. Aşağıdaki form en yaygın olarak kullanılan 4. mertebeden klasik R-K yöntemidir ve artım fonksiyonu ile birlikte
y i
1
y i
1 6
k
1 2
k
2 2
k
3
k
4
h
(25.40)
k
şeklinde yazılabilir. Buradaki k katsayıları aşağıdaki gibi belirlenmiştir; 1
f
(
x i
,
y i
) (25.40a)
k
2
f x i
1 2
h
,
y i
1 2
k
1
h
(25.40b)
k
3
k
4
f x i
1 2
h
,
y i
1 2
k
2
h
(25.40c)
f
x i
h
,
y i
k
3
h
(25.40d)
FIGURE 25.15
4. Dereceden R-K yöntemleriyle karşılaştırıldığında eğim tahminlerinin grafiksel olarak açıklanması.
35
7.4.
4.Mertebeden Klasik Runga-Kutta Yöntemi
denklemleri yukarıdaki gibi verilir. Yalnızca x’in fonksiyonu olan ADD’ler için klasik 4. dereceden R-K yöntemi Simpson’un 1/3 kuralına benzerdir. Ayrıca,, integral aralığı için iyileştirilmiş bir eğim elde etmek amacıyla geliştirilmiş olan çoklu tahmin açısından, 4. dereceden R-K yöntemi Heun yöntemiyle benzerdir. Aşağıdaki şekilde gösterildiği gibi, k’ların her biri bir eğimi göstermektedir. Bu durumda aşağıdaki eşitlik iyileştirilmiş bir eğim elde etmek için bunların bir ağırlıklı ortalamasıdır.
y i
1 1 6
k
1 2
k
2 2
k
3
k
4
h
FIGURE 25.15
Graphical depiction of the slope estimates comprising the fourth-order RK method.
36
7.4.
4.Mertebeden Klasik Runga-Kutta Yöntemi
Example 25.7
Aşağıdaki diferansiyel denklemleri x=0 ile 0.5 aralığında integre etmek için klasik
4.mertebeden R-K yöntemi
ni kullanınız.
(a)
f
(
x
,
y
) 2
x
3 12
x
2 20
x
8 .
5
(b)
h =0.5 adımı ile y (0)=1 başlangıç şartı altında.
f
(
x
,
y
)
4
e
0 .
8
x
0 .
5
y
h =0.5 adımı ile y (0)=2 başlangıç şartı altında.
Çözüm: (a)
(25.40a) ila (25.40d) eşitlikleri kullanılarak k katsayıları hesaplanırsa: k 1 =8.5, k 2 =4.21875, k 3 =4.21875 ve k 4 =1.25 elde edilir ( slayt 38 ). Bu değerler (25.40) eşitliğinde yerlerine yazılırsa;
y
1 6 3.21875
*0.5
0.5
x
4
4
x
3
10
x
2
8.5
x f
(0.5)
3.2188
37
k
1
f
(
x
,
y
) 2
x
3 12
x
2 20
x
8 .
5 ( ,
i i
)
f
(0,1) 2*0 3 12*0 2
k
2
f x i
1 2 8.5
i
1 2
k h
1 =
f
0 1 2 * 0.5, 1 1 *8.5* 0.5
2
f
(0.25, 2.125)
k
2
f
(0.25, 2.125) 2*0.25
3 12*0.25
2
k
3
f x i
1 2 4.2188
i
1 2
k h
2
f
0 1 2 *0.5, 1 1 * 4.2188*0.5
2
f
(0.25,1.0547)
k
3
f
(0.25, 2.0547) 2*0.25
3 12*0.25
2
k
4
i
i
k h
3
f k
4
f
(0.5,1.0273) 2*0.5
3 12*0.5
2 4.2188
f(0.5, 2.1094)
38
7.4.
4.Mertebeden Klasik Runga-Kutta Yöntemi
(b)
Bu kez aralığın başındaki eğim;
k
1
f
(
x i
,
y i
)
f
( 0 , 2 ) 4
e
0 .
8 *( 0 ) 0 .
5 * ( 2 ) 3 olur. Bu değer kullanılarak orta noktadaki (0.5*h=0.25) y değeri ve eğim aşağıdaki gibi hesaplanabilir;
y
( 0 .
25 ) 2 3 * ( 0 .
25 ) 2 .
75
k
2
f x i
1 2
i
1 2
k h
1
f
(0
1 2 0.5, 2
1 2 *3*0.5)
f
(0.25, 2.75)
4
e
0.8*(0.25)
0.5*(2.75)
3.510611
Bu eğim değeri kullanılarak, orta noktada yeni bir eğim değeri hesaplanabilir;
y
( 0 .
25 ) 2 3 .
510611 * ( 0 .
25 ) 2 .
877653
k
3
f x i
1 2
i
1 2
k h
2
f
0 1 2 *0.5, 2 1 *3.510611*0.5
2
f
(0.25, 2.877653) 4
e
0.8*(0.25) 0.5*(2.877653) 3.446785
39
7.4.
4.Mertebeden Klasik Runga-Kutta Yöntemi
Bu eğim kullanılarak aralığın sonundaki y değeri ve eğimi;
k
4
i
i
k h
3
f
4
e
0.8*(0.5)
0.5*(3.723392)
4.105603
şeklinde belirlenir. Bu dört eğimden yararlanılarak, ortalama bir eğim değeri ve bu eğim değeri yoluyla aralığın sonundaki son y değeri aşağıdaki gibi hesaplanır; 1 6
k
1 2
k
2 2
k
3
k
4 1 6 3 2 * ( 3 .
510611 ) 2 * ( 3 .
446785 ) 4 .
105603 3 .
503399
y
( 0 .
5 ) 2 3 .
503399 * ( 0 .
5 ) 3 .
751699 Bu değer analitik değer olan 3.751152 ‘e çok yakındır.
40
7.5 Diferansiyel Denklem Sistemleri (25.4)
Mühendislik problemlerinin bir çoğunda tek bir diferansiyel denklemden ziyade bir diferansiyel denklem takımının çözümü gerekir. Böyle bir denklem takımı aşağıdaki gibi gösterilebilir;
dy
1
dx dy
2
dx dy n dx
f
1 (
x
,
y
1 ,
y
2 , ,
y n
)
f
2 (
x
,
y
1 ,
y
2 , ,
y n
) (25.42)
f n
(
x
,
y
1 ,
y
2 , ,
y n
) Bu şekildeki bir denklem sisteminin çözülebilmesi için n adet başlangıç şartına ihtiyaç vardır.
Bu tür denklem takımlarının çözümü, tek bir denklemi çözmek için kullanılan yöntemlerin n denklem için genişletilmesiyle elde edilir. Burada her bir adım için tüm denklemler sıra ile çözüldükten sonra bir sonraki adıma geçilerek çözüm sürdürülür.
41
7.5.1 Euler Yöntemi
Example 25.9
Aşağıdaki diferansiyel denklem takımını x=0 ile x=2 aralığında verilen şartlar altında Euler yöntemi ile çözünüz.
x=0 için y 1 =4, ve y 2 =6. Adım büyüklüğü h=0.5,
dy
1
dx
0 .
5
y
1
dy dx
2 4 0 .
3
y
2 0 .
1
y
1
Çözüm:
Euler yöntemi her bir değişken için formüle edilirse;
y i
1
y y
1 ,
i
1
i
f y
1 ,
i
(
x i
,
y i
)
h f
1
(
x i
,
y
1
( 0 .
5 )
y
1
4
( 0 )
f
1
0 .
5
y
1 ,
i
,
y
2 ,
i
)
h
* ( 0 ,
y
1
( 4 ) ( 0 ),
* 0 .
5
y
2
( 0 )) 3 *
h
Not: Aşağıda y 2 hesaplanırken, y 1 ‘in yeni değeri kullanılmayacaktır.
42
7.5.1 Euler Yöntemi
y i
1
y y
2 ,
i
1
i
y f
(
x i
,
2 ,
i f y i
)
h
2
(
x i
,
y
1 ,
i
,
y
2 ,
i
)
h
y
2
( 0 .
5 )
6
y
2
( 0 )
4
f
2
0 .
3 ( 0 , *
y
1
( 6 ) ( 0 ),
y
2
0 .
1 * ( 0 )) ( 4 ) *
*
h
0 .
5
6 .
9
İşlemler benzer tarzda devam ettirilirse ( lütfen devam ettirelim );
x
0
y
4
1 y
6
2
0.5
1.0
3 2.25
6.9
7.715
1.5
2.0
1.6875
1.265625
elde edilir.
8.44525
9.094087
43
7.5.2 Runge-Kutta Yöntemi
Herhangi bir R-K yöntemi Euler yöntemindekine benzer mantıkla kolaylıkla uygulanabilir. Burada dikkat edilmesi gereken durum, belirli bir noktada önce bütün denklemler için eğimlerin belirlenip, sonra y değerlerinin belirlenmesidir. Örneğin dördüncü mertebeden R-K yöntemi için, önce bütün denklemler için aralığın başındaki k değerleri (k
1
‘ler) bulunur, sonra aralığın orta noktasında bağımlı değişkenler hesaplanır. Benzer işlemler diğer k ‘lar için de tekrarlanır.
Example 25.10
Yukarıdaki örnekteki denklemleri dördüncü mertebeden R-K yöntemi ile çözünüz.
Çözüm:
Her bir denklem için aralığın başındaki eğimler belirlenirse;
k
1
k
1 , 1
f
(
x i
,
y i
)
f
1
(
x i
,
y
1
i
,
y
2
i
)
k
1,2
0 .
5 * ( 4 )
2
2 ( ,
i
1
i
,
y
2
i
)
f
1
( 0 , 4 , 6 )
f
2 (0, 4, 6) 44
7.5.2 Runge-Kutta Yöntemi
Bu değerlerden yararlanılarak aralığın orta noktasındaki y
1
ve y
2
değeri;
y
2
y
1
k
1,1
k
1,2
h
2
h
2 0.5
3.5
2 0.5
6.45
2 şeklinde belirlenir. Şimdi de bu değerler kullanılarak orta noktadaki eğimlerin ilk değerini hesaplamak için kullanılabilir;
k
2,1
f x
1
(
i
1/ 2
,
y
0.5*(3.5)
,
y
1.75
)
f
1
(0.25,3.5, 6.45)
k
2,2
f x
2
(
i
1/ 2
,
y
,
y
)
f
2
(0.25,3.5, 6.45)
45
7.5.2 Runge-Kutta Yöntemi
Şimdi bu değerler kullanılarak, aralığın orta noktasındaki ikinci y değerleri belirlenebilir;
y
1,
i y
2,
i
k
2,1
k
2,2
h
2
h
2 0.5
3.5625
2 0.5
6.42875
2 Son değerlerden yararlanılarak, aralık ortasındaki eğimlerin ikinci değeri belirlenir;
k
3,1
f x
1
(
i
1/ 2
,
y
,
y
)
f
1
(0.25,3.5625, 6.42875)
0.5*(3.5625)
1.78125
k
3,2
f x
2
(
i
1/ 2
,
y
,
y
)
f
2
(0.25,3.5625, 6.42875)
Bu değerler yardımıyla da aralık sonundaki y değerleri belirlenir; 46
y
1,
i
k h
3,1
y
2,
i
k h
3,2
7.5.2 Runge-Kutta Yöntemi
3.109375
6.857563
ve bu değerler son noktadaki eğimleri hesaplamak için kullanılabilir;
k
4,1
f x
1
(
i
1
,
y
,
y
0.5*(3.109375) )
f
1
(0.5,3.109375, 6.857563) 1.554688
k
4,2
f x
2
(
i
1
,
y
,
y
)
f
2
(0.5,3.109375, 6.857563)
Her bir denklem için belirlenen tüm eğim değerleri ve aralık başındaki y değerleri kullanılarak aralık sonunda bulunmak istene y değerleri hesaplanabilir; 47
7.5.2 Runge-Kutta Yöntemi
y
y
1,
i
1
h
y
1,
i
1 6
k
1,1
2
k
2,1
2
k
3,1
k
4,1
h
6 1
3.115234
y
y
2,
i
2
h
y
2,
i
1 6
k
1,2
2
k
2,2
2
k
3,2
k
4,2
h
1 6
6.857670
İşlemlere benzer şekilde devam edilirse, elde edilen sonuçlar aşağıdaki tabloda görülmektedir.
48
7.5.2 Runge-Kutta Yöntemi
Lütfen bu değerleri işlemleri yaparak elde edelim;
x
0
y
4
1 y
6
2
0.5
1.0
3.115234
2.426171
6.857670
7.632106
1.5
2.0
1.889523
1.471577
8.326886
8.946865
49