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