4 SOLUSI PERSAMAAN DIFERENSIAL BIASA dengan

Seri Mata Kuliah : PEMODELAN dan MATEMATIKA TERAPAN

SOLUSI PERSAMAAN DIFERENSIAL BIASA dengan HARGA AWAL dalam PEMODELAN dan MODEL MATEMATIS

4

4.1. Solusi Persamaan Diferensial Biasa Secara umum, problem persamaan diferensial biasa selalu melibat-kan harga awal, yang dapat ditulis sebagai berikut :

y ' = f ( x , y ), y ( x 0 ) = y 0 x0 ≤ x ≤ xN Secara numerik, solusi yang seringkali diterapkan dalam problem ini adalah dengan metode eksplisit. Dalam hal ini, solusi dari problem di atas adalah berada dalam interval

[x 0 , x N ] yang dibagi secara tetap (equidistance) sebanyak N buah panel: h =

sehingga :

xi = x0 + i h,

xN − x0 N

i = 0 ,1, 2 , K , N

Jika y (x ) adalah “solusi eksak” dari PDB di atas, maka dengan melakukan ekspansi dengan “deret Taylor” dengan sisanya akan diperoleh:

( x i +1 − x i ) 2 y ( x i +1 ) = y ( x i ) + ( x i +1 − x i ) y' ( x i ) + y" ( ξi ), 2! x i ≤ ξi ≤ x i +1 substitusi dari PDB di atas dan dengan pendekatan metode eksplisit pada persamaan di atas akan didapatkan :

h2 y ( x i +1 ) = y ( x i ) + h f ( x i , y ( x i )) + f ' ( ξi , y ( ξ i )) 2!

Property of Setijo Bismo

Halaman (1)


4.2. Metode EULER Metode yang paling sederhana dalam menterjemahkan deret Taylor diatas, yaitu dengan dengan cara “pemotongan” term kedua (atau di atasnya). Dan bila dinotasikan

u i ≈ y ( x i ) , maka : u i +1 = u i + h f ( x i , u i ),

i = 0 ,1 , K , N − 1

u 0 = y0 Dari persamaan di atas, dapat diketahui bahwa untuk menghitung u i +1 diperlukan informasi tentang harga-harga x i dan u i , yang disebut sebagai harga awal. Akurasi dari metode ini adalah dalam “order satu” (first-order approximation) :

ei +1 = 0( h 1 )

Gambar 1. Representasi METODE EULER.

Untuk memperkecil “galat pemotongan lokal” pada setiap tahap (panel), ukuran h dapat dibuat sekecil mungkin.


Halaman (2)


4.3. Metode RUNGE-KUTTA (order 2) Metode ini termasuk algoritma eksplisit yang melibatkan evaluasi fungsi f di antara

xi

dan

x i +1 . Formula umum dari metode ini adalah sebagai berikut :

u i +1 = u i +

ν

∑ ωj K j

j =1

dengan :

K j = h f (x i + c j h , u i + c1 = 0 Perlu dicatat, bila

ν = 1, ω = 1,

dan

j −1

∑ a jl K l )

l =1

K 1 = h f (x i , u i )

maka formula

yang akan diperoleh adalah : METODE EULER. Hal ini berarti bahwa metode EULER adalah order terendah dari METODE RUNGE-KUTTA. Untuk mendapatakan formula dengan akurasi yang lebih tinggi, parameterparameter

ω , c , dan a

dapat diubah dengan tetap melakukan ekspansi solusi eksak

melalui deret Taylor. Sebagai contoh, bila kita inginkan

ν = 2,

maka pertama kali kita ekspansikan

solusi eksak di atas dengan deret Taylor :

h2 y ( x i +1 ) = y ( x i ) + h f ( x i , y ( x i )) + f ' ( x i , y ( x i )) + 0( h 3 ) 2! sedangkan

f ' ( x i , y ( x i )) dapat dituliskan sebagai : df i ∂f ∂f dy = i + i dx ∂x ∂y dx


x =x i

= ( f x + f y f )i

Halaman (3)


Jika persamaan terakhir di atas disubstitusikan kedalam persamaan di atasnya dengan pemotongan term di atasnya, maka akan diperoleh :

u i +1 Untuk ekspansi formula

h2 = u i + h fi + ( f x + f y f )i 2!

Kj

pada posisi ke- i , perlu dicatat bahwa harga K 1

merupakan bentuk paling sederhana, sebagai berikut :

K 1 = h f ( x i , u i ) = h fi sehingga harga K 2 dapat dihitung melalui formula berikut :

K 2 = h f ( x i + c2 h , u i + a 21 K 1 ) Jika diketahui bahwa sembarang dua fungsi berdekatan dengan

xi

η dan φ yang lokasinya berturut-turut

dan u i , maka akan diperoleh :

f ( η, φ) ≈ f ( x i , u i ) + ( η − x i ) f ( x i , u i ) + ( φ − u i ) f y ( x i , u i ) Dengan menggunakan persamaan di atas untuk K 2 , maka akan didapatkan bentuk persamaan untuk menghitung K 2 :

K 2 = h ( fi + c2 h f x + a 21 K 1 f y ) atau

K 2 = hf i + h 2 ( c2 f x + a 21 f y f )i Bila disubstitusikan persamaan terakhir dan persamaan K 1 di atas kedalam formula dasar Runge-Kutta ( u i +1 ) di atas, maka akan diperoleh :

u i +1 = u i + ω1 h fi + ω2 h fi + ω2 h 2 c2 ( f x )i + a 21 ω2 h 2 ( f y f )i


Halaman (4)


Jika dibandingkan persamaan

u i +1

u i +1

yang terakhir ini dengan persaamaan

sebelumnya, maka akan diperoleh :

ω1 + ω2 = 1 ,0 ω2 c2 = 0 ,5 ω2 a 21 = 0 ,5 Algoritma METODE RUNGE-KUTTA dilengkapi dengan cara memilih salah satu di antara paramater-parameter ω1 , ω2 , c2 atau a 21 , sedangkan parameter lainnya ditetapkan dengan menggunakan formula-formula di atas. Jika dipilih c2 = 0 ,5 , maka skema metode RK menjadi METODE TITIK TENGAH (midpoint method) adalah :

u i +1 = u i + h f ( x i + 0 ,5 h , u i + 0 ,5 h fi ),

i = 0 ,1 , K , N − 1

u 0 = y0 Dari persamaan di atas, dapat diketahui bahwa untuk menghitung ‘harga baru’, u i +1 , diperlukan besaran-besaran x i dan u i : • x i dan u i merupakan ‘harga awal’ • harga f i dihitung sebagai fungsi dari x i dan u i : f i = f ( x i , u i ) • h adalah ‘lebar panel’ atau jarak antara u i u i +1 Representasi grafik dari formula metode titik tengah di atas adalah sebagai berikut :

Gambar 2. Representasi METODE TITIK TENGAH. Property of Setijo Bismo

Halaman (5)


Sedangkan jika dipilih c2 = 1 , maka akan diperoleh formulasi kelandaian ratarata (average slope) dari metode Runge-Kutta order-2 :

u i +1 = u i +

h [fi + f ( x i + h , u i + h fi )], 2

i = 0 ,1 , K , N − 1

u 0 = y0 Representasi grafik dari metode kelandaian rata-rata di atas dapat disajikan seperti pada halaman berikut :

Gambar 2. Representasi METODE Kelandaian rata-rata RUNGE-KUTTA order-2.

Kedua skema Metode RK di atas memiliki akurasi order-2, karena pemotongan deret Taylor dilakukan setelah

0(h 2 ) .

Jika diinginkan akurasi pada order-p, maka harus diambil harga

ν

yang sebesar

mungkin. Namun hal ini hampir tidak mungkin atau memerlukan usaha/ pekerjaan yang besar (?).

Tabel 1. hubungan order-p dengan ν .


p

2

3

4

5

6

…

ν

2

3

4

6

8

… Halaman (6)


4.3.1. Peranan harga h dalam solusi numerik Untuk metode RK order-2, harga h sangat berpengaruh dalam peroleh atau solusi numeris dari model yang dimaksudkan, seperti dapat dilihat pada gambar berikut :

4.3.2. Ketelitian beberapa metode numeris Ketelitian dari beberapa metode numeris yang umum digunakan dapat dilihat pada gambar berikut :


Halaman (7)


4.4. Metode RUNGE-KUTTA order tinggi 1. Metode RUNGE-KUTTA-GILL: Metode Runge-Kutta-Gill (RKG) tergolong dalam keluarga metode RK order-4, yang memiliki 4 (empat) buah ‘konstanta perhitungan antara’ yang dikombinasikan dengan konstanta-konstanta lain (a, b, c, dan d) sebagai keluarga bilangan emas (golden numbers). Algoritma ringkas dari metode RKG ini dapat dituliskan seperti di bawah ini :

1 6

1 3

u i +1 = u i + ( K 1 + K 4 ) + (b K 2 + d K 3 ) K 1 = h f (x i ,u i ) K 2 = h f ( x i + 12 h , u i + 12 K 1 ) K 3 = h f ( x i + 12 h , u i + a K 1 + b K 2 ) K 4 = h f (x i + h ,u i + c K 2 + d K 3 ) a = c =

2 −1 , 2 2 − , 2

b =

2− 2 2

d = 1+

2 2

untuk : i = 0, 1, 2,…,N-1 dan harga awal : u 0 = y0


Halaman (8)


2. Metode RUNGE-KUTTA-MERSON : Metode Runge-Kutta-Merson (RKM) tergolong dalam keluarga metode RungeKutta order-4, namun memiliki ketelitian sampai order-5. Keistimewaan ini dimungkinkan karena metode RKM memiliki 5 (lima) buah ‘konstanta perhitungan antara’ yang berperan untuk memprediksi harga solusi yang diinginkan pada 2 (dua) keadaan sedemikian rupa sehingga ‘galat pembulatan’ dapat diminimisasi sampai order-5. Formulasi ringkas dari metode RKM ini dapat dituliskan seperti di bawah ini : 1 2

3 2

1 6

2 3

u i +1 = u i + K 1 − K 3 + 2 K 4 1 6

u i +1 = u i + K 1 + K 4 + K 5 K 1 = h f (x i ,u i ) K 2 = h f ( x i + 13 h , u i + 13 K 1 ) K 3 = h f ( x i + 13 h , u i + 16 K 1 + 16 K 2 ) K 4 = h f ( x i + 12 h , u i + 18 K 1 + 38 K 3 ) K 5 = h f ( x i + h , u i + 12 K 1 − 32 K 3 + 2 K 4 ) untuk :

i = 0, 1, 2,…,N-1 dan harga kondisi awal :

u 0 = y0


Halaman (9)


3. Metode RUNGE-KUTTA-FEHLBERG : Sama halnya dengan metode RKM, metode Runge-Kutta-Fehlberg (RKF45) juga tergolong dalam keluarga metode Runge-Kutta order-4, namun memiliki ketelitian sampai order-5. Ketelitian yang tinggi ini dimungkinkan karena metode RKF45 memiliki 6 (enam) buah ‘konstanta perhitungan antara’ yang berperan untuk meng-update solusi sampai order-5. Formulasi ringkas dari metode RKM ini adalah :

K 1 = h f (x i ,u i ) K 2 = h f ( x i + 14 h , u i + 14 K 1 ) K 3 = h f ( x i + 38 h , u i + 323 K 1 + 329 K 2 ) K 4 = h f ( x i + 12 h , u i + 1932 K 1 − 7200 K 2 + 7296 K3) 13 2197 2197 2197 845 K 5 = h f ( x i + h , u i + 439 K 1 − 8 K 2 + 3680 K 3 − 4104 K4) 216 513

K 6 = h f ( x i + 12 h , u i − 278 K 1 + 2 K 2 − 3544 K 3 + 1859 K 4 − 11 K5) 2565 4104 40 • Formula ‘update’ order-4 :

u i +1 = u i +

25 216

K1 +

1408 2565

K1 +

6656 12825

K3 +

2197 4104

1 5

K4 − K5

• Formula order-5 :

uˆ i +1 = u i +

16 135

K3 +

28561 56437

K4 −

9 50

K5 +

2 55

K6

• Galat ‘pembabatan’ order-4 :

uˆ i +1 − u i +1 = untuk :


1 360

K1 −

128 4275

K3 −

2197 75240

K4 +

1 50

K5 +

2 55

K6

i = 0, 1, 2,…,N-1 dan u 0 = y 0

Halaman (10)


Contoh-contoh Integrasi Numerik : Persamaan Tunggal :

dy = − 25 y dt Solusi Eksak (Analitis) :

y = e − 25 t Tabel 2. Hasil perhitungan persamaan tunggal dengan metode RK order-2. x

Nilai Solusi Eksak (y)

Contoh 1 (RKSR)

Contoh 2 (RKMP)

Contoh 3 (RKG)

0.00

1.0000E+00

1.0000E+00

1.0000E+00

1.0000E+00

0.20

6.7379E-03

6.7415E-03

6.7415E-03

6.7379E-03

0.40

4.5400E-05

4.5448E-05

4.5448E-05

4.5400E-05

0.60

3.0590E-07

3.0639E-07

3.0639E-07

3.0590E-07

0.80

2.0612E-09

2.0655E-09

2.0655E-09

2.0612E-09

1.00

1.3888E-11

1.3925E-11

1.3925E-11

1.3888E-11


Hasil Perhitungan Numerik (ui+1)

Halaman (11)


Contoh 1 : {Program Solusi Persamaan Diferensial Biasa dengan Metode RUNGE-KUTA 'Slope Rata-rata'} Type Real = Extended; Real01 = Array [0..1] of Real; Procedure F(x,Y : Real; Var DY : Real); Begin DY := -25*Y; End; Procedure DRK2SR(XV : Real01; YI : Real; Var YF : Real; NP : Integer; Eps : Real); {----------------------------------------------------PROGRAM : 2-nd ORDER 'MEAN SLOPE' RUNGE-KUTTA METHOD FOR ORDINARY DIFFERENTIAL EQUATION (M.E. Davis, p13) F : Function F(x,y) to be integrated N : Number of differential equations XV : Vector of xv[0]=initial and xv[1]=final YI,YF : Values of Y-initial and Y-final H : Step length -----------------------------------------------------} Var H,X,YP : Real; K1,K2 : Real; Begin H := (xv[1] - xv[0])/NP; X := XV[0]; YF := YI; Repeat YP := YF; F(X,YP,K1); YF := YP + H*K1; F(X+H,YF,K2); YF := YP + H*(K1 + K2)/2; X := X + H; Until (ABS(XV[1]-X) <= EPS); End; Var I,NP : Integer; XV : Real01; Yi,Yf,Eps,h : Real; Begin Eps := 1.0E-6; NP := 200; xv[0] := 0; xv[1] := 0.2; Yi := 1.0; Writeln(xv[0]:0:4,' ',Yi:13,' ',exp(-25*xv[0]):13); For I := 1 to 5 do Begin DRK2SR(xv,Yi,Yf,NP,Eps); Writeln(xv[1]:0:4,' ',Yf:13,' ',exp(-25*xv[1]):13); xv[0] := xv[1]; xv[1] := xv[1] + 0.2; Yi := Yf; End; Readln; End.


Halaman (12)


Contoh 2 : {Program Solusi Persamaan Diferensial Biasa dengan Metode RUNGE-KUTTA 'Midpoint'} Type Real = Extended; Real01 = Array [0..1] of Real; Real50 = Array [1..50] of Real; Procedure F(x,Y : Real; Var DY : Real); Begin DY := -25*Y; End; Procedure DRK2MP(XV : Real01; YI : Real; Var YF : Real; NP : Integer; Eps : Real); {----------------------------------------------------PROGRAM : 2-nd ORDER 'MIDPOINT' RUNGE-KUTTA METHOD FOR ORDINARY DIFFERENTIAL EQUATION (M.E. Davis, p13) F : Function F(x,y) to be integrated XV : Vector of xv[0]=initial and xv[1]=final YI,YF : Value of Y-initial and Y-final NP : Number of panels H : Step length -----------------------------------------------------} Var H,X,YP : Real; I : Integer; K1,K2 : Real; Begin H := (xv[1] - xv[0])/NP; X := XV[0]; YF := YI; Repeat YP := YF; F(X,YP,K1); YF := YP + H*K1/2; F(X+H/2,YF,K2); YF := YP + H*K2; X := X + H; Until (ABS(XV[1]-X) <= EPS); End; Var I,NP : Integer; XV : Real01; Yi,Yf : Real; Eps : Real; Begin Eps := 1.0E-6; NP := 200; xv[0] := 0; xv[1] := 0.2; Yi := 1.0; Writeln(xv[0]:0:4,' ',Yi:13,' ',exp(-25*xv[0]):13); For I := 1 to 5 do Begin DRK2MP(xv,Yi,Yf,NP,Eps); Writeln(xv[1]:0:4,' ',Yf:13,' ',exp(-25*xv[1]):13); xv[0] := xv[1]; xv[1] := xv[1] + 0.2; Yi := Yf; End; Readln; End.


Halaman (13)


Contoh 3 : {Program Solusi Persamaan Diferensial Biasa dengan Metode RUNGE-KUTA-GILL} Type Real = Extended; Real01 = Array [0..1] of Real; Procedure F(x,Y : Real; Var DY : Real); Begin DY := -25.0*Y; End; Procedure DRKGIL(XV : Real01; YI : Real; Var YF : Real; NP : Integer; Eps : Real); {----------------------------------------------------PROGRAM : 4-th ORDER RUNGE-KUTTA-GILL METHOD FOR ORDINARY DIFFERENTIAL EQUATION (M.E. Davis, p13) F : Function F(x,y) to be integrated N : Number of differential equations XV : Vector of xv[0]=initial and xv[1]=final YI,YF : Value of Y-initial and Y-final H : Step length -----------------------------------------------------} Var a,b,c,d,H,X,YP : Real; I : Integer; K1,K2,K3,K4 : Real; Begin a := (Sqrt(2)-1)/2; b := (2-Sqrt(2))/2; c := -Sqrt(2)/2; d := 1 + Sqrt(2)/2; H := (xv[1] - xv[0])/NP; X := XV[0]; YF := YI;


Repeat YP := YF; F(X,YP,K1); YF := YP + H*K1/2; F(X+H/2,YF,K2); YF := YP + H*(a*K1+b*K2); F(X+H/2,YF,K3); YF := YP + H*(c*K2+d*K3); F(X+H,YF,K4); X := X + H; YF := YP + H*((K1+K4)/6 + (b*K2+d*K3)/3); Until (ABS(XV[1]-X) <= EPS); End; Var I,NP : Integer; XV : Real01; Yi,Yf : Real; Eps : Real; Begin Eps := 1.0E-6; NP := 200; xv[0] := 0; xv[1] := 0.2; Yi := 1.0; Writeln(xv[0]:0:3,' ',Yi:13,' ',exp(-25*xv[0]):13); For I := 1 to 5 do Begin DRKGIL(xv,Yi,Yf,NP,Eps); Writeln(xv[1]:0:3,' ',Yf:13,' ',exp(-25*xv[1]):13); xv[0] := xv[1]; xv[1] := xv[1] + 0.2; Yi := Yf; End; Readln; End.

Halaman (14)


Persamaan Jamak (Sistem PDB) : dy  3 ,21  = − 0 ,1744 exp  *  y dx T  dT *  3 ,21  = 0 ,06984 exp  *  y dx T  Tabel 3. Hasil perhitungan persamaan jamak dengan metode RK order-2.

Hasil Perhitungan Numerik untuk Sistem PDB : x

Contoh 4 : RK-Slope rata2

Contoh 5 : RK-Midpoint

y

T

y

T

0.000

1.000000

1.000000

1.000000

1.000000

0.100

0.700397

1.119979

0.700430

1.119966

0.200

0.529232

1.188523

0.529259

1.188512

0.300

0.413765

1.234763

0.413787

1.234754

0.400

0.329942

1.268331

0.329959

1.268324

0.500

0.266512

1.293732

0.266525

1.293726

0.600

0.217224

1.313469

0.217235

1.313465

0.700

0.178223

1.329088

0.178232

1.329084

0.800

0.146954

1.341610

0.146961

1.341607

0.900

0.121638

1.351748

0.121644

1.351745

1.000

0.100989

1.360017

0.100993

1.360015


Halaman (15)


Contoh 4 : {Program Solusi Sistem Persamaan Diferensial Biasa dengan Metode RUNGE-KUTA 'Slope Rata-rata'} Type Real = Extended; Real01 = Array [0..1] of Real; Real50 = Array [1..50] of Real; Procedure F(N : Integer; x : Real; Y : Real50; Var DY : Real50); Begin DY[1] := -0.1744*exp(3.21/Y[2])*Y[1]; DY[2] := 0.06984*exp(3.21/Y[2])*Y[1]; End; {$I DRK2SR} Var

I,N,NP : Integer; XV : Real01; Yi,Yf : Real50; Eps : Real;

Begin Eps := 1.0E-6; N := 2; NP := 20; xv[0] := 0; xv[1] := 0.1; Yi[1] := 1.0; Yi[2] := 1.0; Writeln(xv[0]:0:3,' ',Yi[1]:0:6,' For I := 1 to 10 do Begin DRK2SR(N,xv,Yi,Yf,NP,Eps); Writeln(xv[1]:0:3,' ',Yf[1]:0:6,' xv[0] := xv[1]; xv[1] := xv[1] + 0.1; Yi := Yf; End; Readln; End.


',Yi[2]:0:6);

',Yf[2]:0:6);

Halaman (16)


Contoh 5 : {Program Solusi Sistem Persamaan Diferensial Biasa dengan Metode RUNGE-KUTTA 'Midpoint'} Type Real = Extended; Real01 = Array [0..1] of Real; Real50 = Array [1..50] of Real; Procedure F(N : Integer; x : Real; Y : Real50; Var DY : Real50); Begin DY[1] := -0.1744*exp(3.21/Y[2])*Y[1]; DY[2] := 0.06984*exp(3.21/Y[2])*Y[1]; End; {$I DRK2MP} Var I,N,NP : Integer; XV : Real01; Yi,Yf : Real50; Eps : Real; Begin Eps := 1.0E-6; N := 2; NP := 20; xv[0] := 0; xv[1] := 0.1; Yi[1] := 1.0; Yi[2] := 1.0; Writeln(xv[0]:0:3,' ',Yi[1]:0:6,' For I := 1 to 10 do Begin DRK2MP(N,xv,Yi,Yf,NP,Eps); Writeln(xv[1]:0:3,' ',Yf[1]:0:6,' xv[0] := xv[1]; xv[1] := xv[1] + 0.1; Yi := Yf; End; Readln; End.


',Yi[2]:0:6);

',Yf[2]:0:6);

Halaman (17)


Tabel 4. Perbandingan hasil perhitungan sistem PDB jamak antara metode-metode order-4 : RK-Gill dan RK-Merson.

Hasil Perhitungan Numerik untuk Sistem PDB :

x

Contoh 6 : RK-Gill

Contoh 7 : RK-Merson

y

T

y

T

0.000

1.000000

1.000000

1.000000

1.000000

0.100

0.700372

1.119989

0.700372

1.119989

0.200

0.529209

1.188532

0.529209

1.188532

0.300

0.413745

1.234771

0.413746

1.234771

0.400

0.329925

1.268337

0.329925

1.268337

0.500

0.266497

1.293738

0.266497

1.293738

0.600

0.217212

1.313474

0.217212

1.313474

0.700

0.178213

1.329092

0.178213

1.329092

0.800

0.146945

1.341613

0.146945

1.341613

0.900

0.121631

1.351751

0.121631

1.351751

1.000

0.100982

1.360020

0.100982

1.360020


Halaman (18)


Contoh 6 : {Program Solusi Sistem Persamaan Diferensial Biasa dengan Metode RUNGE-KUTA-GILL} Type Real = Extended; Real01 = Array [0..1] of Real; Real50 = Array [1..50] of Real; Procedure F(N : Integer; x : Real; Y : Real50; Var DY : Real50); Begin DY[1] := -0.1744*exp(3.21/Y[2])*Y[1]; DY[2] := 0.06984*exp(3.21/Y[2])*Y[1]; End; {$I DRKGIL} Var I,N,NP : Integer; XV : Real01; Yi,Yf : Real50; Eps : Real; Begin Eps := 1.0E-6; N := 2; NP := 10; xv[0] := 0; xv[1] := 0.1; Yi[1] := 1.0; Yi[2] := 1.0; Writeln(xv[0]:0:3,' ',Yi[1]:0:6,' For I := 1 to 10 do Begin DRKGIL(N,xv,Yi,Yf,NP,Eps); Writeln(xv[1]:0:3,' ',Yf[1]:0:6,' xv[0] := xv[1]; xv[1] := xv[1] + 0.1; Yi := Yf; End; Readln; End.


',Yi[2]:0:6);

',Yf[2]:0:6);

Halaman (19)


Contoh 7 : {Program Solusi Sistem Persamaan Diferensial Biasa dengan Metode RUNGE-KUTA-MERSON} Type Real = Real01 Real02 Real50

Extended; = Array [0..1] of Real; = Array [0..2] of Real; = Array [1..50] of Real;

Procedure F(N : Integer; x : Real; Y : Real50; Var DY : Real50); Begin DY[1] := -0.1744*exp(3.21/Y[2])*Y[1]; DY[2] := 0.06984*exp(3.21/Y[2])*Y[1]; End; {$I DRKMER} Var I,MSG,N,NP : Integer; xv : Real01; s : Real02; Yi,Yf : Real50; Eps : Real; Begin Eps := 1.0E-4; N := 2; NP := 100; xv[0] := 0; xv[1] := 0.1; s[0] := (xv[1] - xv[0])/NP; s[1] := s[0]/64; Yi[1] := 1.0; Yi[2] := 1.0; Writeln(xv[0]:0:3,' ',Yi[1]:0:6,' ',Yi[2]:0:6); For I := 1 to 10 do Begin DRKMER(N,xv,Yi,Yf,s,Eps,MSG); Writeln(xv[1]:0:3,' ',Yf[1]:0:6,' ',Yf[2]:0:6,' ',MSG); xv[0] := xv[1]; xv[1] := xv[1] + 0.1; Yi := Yf; End; Readln; End.


Halaman (20)


Tabel 5. Perbandingan hasil perhitungan sistem PDB jamak antara metode-metode order-4 : RK-Gill dan RK-Fehlberg. Hasil Perhitungan Numerik untuk Sistem PDB : Contoh 6 : RK-Gill

x

Contoh 8 : RK-Fehlberg y T

y

T

0.000

1.000000

1.000000

1.000000

1.000000

0.100

0.700372

1.119989

0.700372

1.119989

0.200

0.529209

1.188532

0.529209

1.188532

0.300

0.413745

1.234771

0.413745

1.234771

0.400

0.329925

1.268337

0.329925

1.268337

0.500

0.266497

1.293738

0.266497

1.293738

0.600

0.217212

1.313474

0.217212

1.313474

0.700

0.178213

1.329092

0.178213

1.329092

0.800

0.146945

1.341613

0.146945

1.341613

0.900

0.121631

1.351751

0.121631

1.351751

1.000

0.100982

1.360020

0.100982

1.360020


Halaman (21)


Contoh 8 : {Program Solusi Sistem Persamaan Diferensial Biasa dengan Metode RUNGE-KUTA-FEHLBERG} Type Real = Real01 Real02 Real50


Procedure F(N : Integer; x : Real; Y : Real50; Var DY : Real50); Begin DY[1] := -0.1744*exp(3.21/Y[2])*Y[1]; DY[2] := 0.06984*exp(3.21/Y[2])*Y[1]; End; {$I DRKF45} Var I,MSG,N,NP : Integer; xv : Real01; s : Real02; Yi,Yf : Real50; Eps : Real; Begin Eps := 1.0E-5; N := 2; NP := 100; xv[0] := 0; xv[1] := 0.1; s[0] := (xv[1] - xv[0])/NP; s[1] := s[0]/64; Yi[1] := 1.0; Yi[2] := 1.0; Writeln(xv[0]:0:3,' ',Yi[1]:0:6,' ',Yi[2]:0:6); For I := 1 to 10 do Begin DRKF45(N,xv,Yi,Yf,s,Eps,MSG); Writeln(xv[1]:0:3,' ',Yf[1]:0:6,' ',Yf[2]:0:6); xv[0] := xv[1]; xv[1] := xv[1] + 0.1; Yi := Yf; End; Readln; End.


Halaman (22)


Solusi PDB Order-2 menggunakan Metode RKG dan RKM dengan ‘Teknik Shooting’ : Persamaan Diferensial Biasa Order-2 non-linier berikut :

d2y dx

2

= −1 − ( x 2 +1 ) y

y ( −1) = 0 dan y (1) = 0 PDB di atas memiliki informasi tentang harga-harga fungsi pada x = −1 dan x = 1 , akan tetapi tidak memiliki informasi yang memadai tentang harga-harga awalnya untuk turunan pertama dan kedua (karena merupakan PDB order-2). Salah satu solusi yang mungkin dilakukan adalah dengan metode ‘trial and error’, sehingga PDB order-2 di atas dapat diubah menjadi Sistem PDB berikut :

 dy 1  dx = y 2   dy 2 = − 1 − ( x 2 + 1 ) y 1  dx

 y1 ( −1) = 0    y ( −1) = input  2

Strategi solusi PDB yang mungkin dilakukan dapat dijelaskan sebagai berikut :

• Karena PDB tunggal di atas berbentuk order-2, maka dimisalkan suatu sistem PDB baru dengan 2 (dua) buah variabel terikat, yaitu dalam y1 dan y 2 , • Harga awal (kondisi awal) dari y1 diketahui, namun untuk y 2 tidak sehingga dalam hal ini perlu diberikan dengan cara ‘trial and error’ • Integrasi PDB dapat dilakukan untuk informasi yang diketahui, dalam hal ini antara x = −1 sampai x = 1 .


Halaman (23)


Tabel 6. Perbandingan hasil perhitungan sistem PDB order-2 dengan ‘teknik shooting’ antara metode RK-Gill dan RK-Merson. x

Contoh 9 : Metode RKG

Contoh 10 : Metode RKM

y1

y2

y1

y2

-1.00

0.000000

1.736465

0.000000

1.736465

-0.90

0.168104

1.620549

0.168104

1.620549

-0.80

0.323231

1.478228

0.323231

1.478228

-0.70

0.463112

1.316726

0.463112

1.316726

-0.60

0.586136

1.141981

0.586136

1.141981

-0.50

0.691223

0.958637

0.691223

0.958637

-0.40

0.777693

0.770133

0.777693

0.770133

-0.30

0.845157

0.578843

0.845157

0.578843

-0.20

0.893419

0.386258

0.893419

0.386258

-0.10

0.922394

0.193192

0.922394

0.193192

0.00

0.932054

-0.000000

0.932054

-0.000000

0.10

0.922394

-0.193192

0.922394

-0.193192

0.20

0.893419

-0.386259

0.893419

-0.386259

0.30

0.845157

-0.578843

0.845157

-0.578843

0.40

0.777693

-0.770133

0.777693

-0.770133

0.50

0.691223

-0.958637

0.691223

-0.958637

0.60

0.586136

-1.141981

0.586136

-1.141981

0.70

0.463112

-1.316727

0.463112

-1.316727

0.80

0.323231

-1.478228

0.323231

-1.478228

0.90

0.168104

-1.620549

0.168104

-1.620549

1.00

-0.000000

-1.736465

-0.000000

-1.736465


Halaman (24)


Contoh 9 : {Program Solusi Sistem PDB 'turunan kedua' (order 2) dengan 'Shooting' menggunakan Metode RUNGE-KUTA-GILL PDB order 2 : d2Y/dx2 = -1 - (x^2 + 1)*Y x(-1) = 0 dan x(1) = 0} Type Real = Extended; Real01 = Array [0..1] of Real; Real50 = Array [1..50] of Real; Procedure F(N : Integer; x : Real; Y : Real50; Var DY : Real50); Begin DY[1] := Y[2]; DY[2] := -1 - (Sqr(x) + 1)*Y[1]; End; {$I DRKGIL} Var I,N,NP : Integer; XV : Real01; Yi,Yf : Real50; Eps : Real; Begin Eps := 1.0E-6; N := 2; NP := 10; xv[0] := -1; xv[1] := -0.9; Yi[1] := 0.0; Yi[2] := 1.736465; Writeln(xv[0]:0:3,' ',Yi[1]:0:6,' ',Yi[2]:0:6); For I := 1 to 20 do Begin DRKGIL(N,xv,Yi,Yf,NP,Eps); Writeln(xv[1]:0:3,' ',Yf[1]:0:6,' ',Yf[2]:0:6); xv[0] := xv[1]; xv[1] := xv[1] + 0.1; Yi := Yf; End; Readln; End.


Halaman (25)


Contoh 10 : {Program Solusi Sistem PDB 'turunan kedua' (order 2) dengan 'Shooting' menggunakan Metode RUNGE-KUTA-MERSON PDB order 2 : d2Y/dx2 = -1 - (x^2 + 1)*Y x(-1) = 0 x(1) = 0} Type Real = Real01 Real02 Real50


Procedure F(N : Integer; x : Real; Y : Real50; Var DY : Real50); Begin DY[1] := Y[2]; DY[2] := -1 - (Sqr(x) + 1)*Y[1]; End; {$I DRKMER} Var I,MSG,N,NP : Integer; xv : Real01; s : Real02; Yi,Yf : Real50; Eps : Real; Begin Eps := 1.0E-4; N := 2; NP := 10; xv[0] := -1.0; xv[1] := -0.9; s[0] := (xv[1] - xv[0])/NP; s[1] := s[0]/64; Yi[1] := 0.0; Yi[2] := 1.736465; Writeln(xv[0]:0:3,' ',Yi[1]:0:6,' ',Yi[2]:0:6); For I := 1 to 20 do Begin DRKMER(N,xv,Yi,Yf,s,Eps,MSG); Writeln(xv[1]:0:3,' ',Yf[1]:0:6,' ',Yf[2]:0:6); xv[0] := xv[1]; xv[1] := xv[1] + 0.1; Yi := Yf; End; Readln; End.


Halaman (26)


Aplikasi Solusi PDB tunggal Order-2 dengan ‘Teknik Substitusi Variabel’ (Rice & Duong Do, hal. 229-230) menggunakan Metode Runge-Kutta-Gill : Persamaan Diferensial Biasa Order-2 non-linier berikut :

y ′′ + y ′ = − sin ( x ) + cos ( x ) y(0 ) = 0 ;

y ′(0) = 1

Harga-harga (kondisi) awal dari PDB tunggal di atas dapat memadai, jika persamaan tersebut dikonversikan menjadi formula baku permisalan Sistem PDB berikut (untuk order-2) :  dy 1  dx = 1   dy 2 = y3  dx   dy 3  dx = − sin ( y1 ) + cos( y1 ) - y 3 

dy 1   y1 = x ⇒ dx = 1   y2 = y  dy  y3 = 2 dx 

 y1 (0) = 0    y 2 (0) = 0    y 3 (0) = 1

Langkah-langkah permisalan di atas dapat dilakukan berdasarkan algoritma berikut :

• Permisalan dimulai pada ‘variabel bebas’ x , untuk y1 sehingga diketahui harga turunannya (lihat kolom tengah), • Permisalan kedua adalah untuk ‘variabel terikat’ y , untuk y 2 , • Permisalan ketiga (yang terakhir) adalah untuk ‘turunan variabel terikat’ ( dy dx = dy 2 dx ) sebagai y 3 , • Sistem PDB baru yang diperoleh adalah dengan cara menyusun turunan-turunan dari variabel-variabel permisalan di atas ( y1 dan y 2 ) yang digabungkan dengan penyusunan ulang PDB tunggal di atas untuk dy 3 dx (lihat kolom kiri) • Harga-harga atau kondisi awal dari sistem PDB di atas berturut-turut merupakan harga-harga pada x = 0 (dalam hal ini y1 (0) ), variabel terikat pada saat x = 0 (dalam hal ini y 2 (0) ), dan turunan pertama variabel terikat pada saat x = 0 (dalam hal ini y 3 (0) ). Property of Setijo Bismo

Halaman (27)


• Solusi yang diinginkan adalah harga-harga y 2 sebagai fungsi x . Perbandingan harga-harga solusi numerik dengan metode RKG dari PDB di atas ( y 2 ) dengan solusi eksak disajikan pada tabel di bawah ini : Tabel 7. Perbandingan solusi numerik sistem PDB order-2 dengan metode RK-Gill dengan solusi analitis.

x 0.000 0.157 0.314 0.471 0.628 0.785 0.942 1.100 1.257 1.414 1.571 1.728 1.885 2.042 2.199 2.356 2.513 2.670 2.827 2.985 3.142


Metode Runge-Kuta-Gill (Contoh 11) :

Solusi Eksak

y1

y2

y3

y

0.000000 0.157080 0.314159 0.471239 0.628319 0.785398 0.942478 1.099557 1.256637 1.413717 1.570796 1.727876 1.884956 2.042035 2.199115 2.356194 2.513274 2.670354 2.827433 2.984513 3.141593

0.000000 0.156434 0.309017 0.453990 0.587785 0.707107 0.809017 0.891007 0.951057 0.987688 1.000000 0.987688 0.951057 0.891007 0.809017 0.707107 0.587785 0.453990 0.309017 0.156434 -0.000000

1.000000 0.987688 0.951057 0.891007 0.809017 0.707107 0.587785 0.453990 0.309017 0.156434 -0.000000 -0.156434 -0.309017 -0.453990 -0.587785 -0.707107 -0.809017 -0.891007 -0.951057 -0.987688 -1.000000

0.000000 0.156434 0.309017 0.453990 0.587785 0.707107 0.809017 0.891007 0.951057 0.987688 1.000000 0.987688 0.951057 0.891007 0.809017 0.707107 0.587785 0.453990 0.309017 0.156434 0.000000

Halaman (28)


Contoh 11 : {Program Solusi Sistem PDB 'turunan kedua' (order 2) dengan 'konversi' atau 'substitusi variabel' menggunakan Metode RUNGE-KUTA-GILL PDB order 2 : y" + y' = - sin(x) + cos(x) y1(0) = 0 y2(0) = 0 dan y3(0) = 1} Type Real = Extended; Real01 = Array [0..1] of Real; Real50 = Array [1..50] of Real; Procedure F(N : Integer; x : Real; Y : Real50; Var DY : Real50); Begin DY[1] := 1; DY[2] := Y[3]; DY[3] := -Sin(Y[1]) + Cos(Y[1]) - Y[3]; End; {$I DRKGIL} Var I,N,NP : Integer; XV : Real01; Yi,Yf : Real50; Eps,Pi : Real; Begin Pi := 4*ArcTan(1); Eps := 1.0E-6; N := 3; NP := 20; xv[0] := 0; xv[1] := Pi/20; Yi[1] := 0.0; Yi[2] := 0.0; Yi[3] := 1.0; Writeln(xv[0]:0:3,' ',Yi[1]:0:6,' ',Yi[2]:0:6, ' ',Yi[3]:0:6,' ',Sin(xv[0]):0:6); For I := 1 to 20 do Begin DRKGIL(N,xv,Yi,Yf,NP,Eps); Writeln(xv[1]:0:3,' ',Yf[1]:0:6,' ',Yf[2]:0:6, ' ',Yf[3]:0:6,' ',Sin(xv[1]):0:6); xv[0] := xv[1]; xv[1] := xv[1] + Pi/20; Yi := Yf; End; Readln; End.


Halaman (29)

4 SOLUSI PERSAMAAN DIFERENSIAL BIASA dengan

Recommend Documents