PRINSIP DASAR PEMODELAN dan MODEL MATEMATIS

Seri Mata Kuliah : PEMODELAN dan MATEMATIKA TERAPAN

PRINSIP DASAR PEMODELAN

1

dan MODEL MATEMATIS

1.1. Prinsip Dasar Pemodelan Secara fundamental, pemodelan di dalam kajian-kajian proses teknik kimia dan proses adalah : • Penggambaran kinerja suatu aktivitas, sistem atau proses • Membangun persamaan matematis yang dapat menggambar-kan kinerja suatu proses (secara fisik)

1.2. Persamaan Matematis Didalam aktivitas pemodelan di dalam permasalahan teknik kimia dan proses, umumnya dihasilkan suatu bentuk atau sistem persamaan matematis. Secara garis besar, bentuk-bentuk persamaan yang mungkin terbentuk adalah : 1. Persamaan Aljabar : manakala proses berlangsung secara tunak atau penggambaran kinerja proses statik. – Hubungan antar variabel : linier atau non-linier (PAL atau PANL) – Jumlah persamaan (variabel anu) : tunggal atau jamak/serempak (PA atau SPA) – Pengungkapan : eksplisit atau implisit. 2. Persamaan Diferensial : bila proses yang digambar-kan berlangsung secara dinamis (unsteady state process, time dependent proses) : – Hubungan antar variabel : linier atau non-linier – Jumlah persamaan (jumlah variabel terikat yang dideferensialkan) : tunggal atau jamak – Dimensi perubahan (dinamisasi variabel) : biasa (PDB) atau parsial (PDP) Property of Setijo Bismo

Halaman (1)


– Pengungkapan : eksplisit atau implisit.

Persamaan Aljabar • K =yx

• q =mc p ΔT +mλ ⎧ y1 = x1 +2x 2 2 ⎩ y2 = x1 ⋅x 2 −x 2

• ⎨

Persamaan Diferensial • dx A dt = k ⋅ C A ,0 ⋅ ( 1 − x A ) d (N X C )= FR X CR + R−V YC dt ⎧ dy1 dθ = k1 y1 ⋅ y2 −k2 y22 • ⎨ 2 ⎩dy2 dθ = k3 y1

•

1.3. Strategi Pemahaman Suatu Model Matematis Secara konseptual, diperlukan pemahaman yang mendasar tentang persamaan-persamaan model yang terbentuk, sebagai berikut : A. Tidak mungkin semua variabel tidak diketahui (harganya) dan tidak mungkin semua variabel/ besaran diketahui (dull equation). B. Persamaan Tunggal : hanya satu variabel yang harus dihitung, simbol atau besaran lainnya disebut konstanta atau paramater (yang diketahui harganya). C. Persamaan Jamak (n buah) : hanya n buah variabel yang harus dihitung, teliti dan pelajari parameter/ besaran lain yang berperan sebagai konstanta. Penyelesaian model ini umumnya memerlukan harga awal atau tebakan, dilakukan secara iteratif dan serempak. D. Jika pembentukan model telah sesuai dengan kaidah dan sistematika yang benar : solusi dapat terarah (konvergen).

1.4. Penalaran atau Aliran Logika Pemodelan Model matematis yang terbentuk disyaratkan harus memenuhi aliran logika, baik yang bersifat fisika maupun matematis. Pada Property of Setijo Bismo

Halaman (2)


dasarnya, model matematis tersebut harus mampu menggambarkan diagram aliran informasi dasar dari permasalahan yang dimaksudkan. Di bawah ini, secara sederhana diberikan suatu bentuk persamaan atau model matematis sebagai berikut :

d {N X C } = FR X CR + R − V YC dt Sebagai Persamaan Diferensial Biasa Eksplisit, persamaan di atas memiliki 2 kelompok posisi variabel/parameter/ konstanta, yaitu satu fihak di RUAS KIRI (N dan XC) dan di fihak lain berada pada RUAS KANAN (FR, XCR, R, V dan YC). Bila diinginkan mencari atau menyelesaikan persamaan di atas, maka perlu difahami diagram aliran logika berikut : R FR, XCR V, YC

d {N XC } = FR XCR + R − V YC dt

XC

N

Property of Setijo Bismo

Halaman (3)


BEBERAPA CONTOH MODEL Kasus #1 :

F1 ZZ0

Z F2

Mempunyai DIAGRAM ALIR INFORMASI DASAR sbb : Masukan

Persamaan-persamaan Sistem

F1, F2

dZ A = F1 − F2 dt

A

dZ dt

Z =

Keluaran

⎧ dZ ⎫ ⎬ dt dt ⎭

∫ ⎨⎩

Z

Blok Diagram Informasi Dasar

Secara lebih ringkas, diagram alir kerja MODEL dari Kasus #1 ini adalah sbb : Masukan

Persamaan Sistem

F1, F2 A

A

dZ = F1 − F2 dt

Keluaran

Z

Model Kasus #1


Halaman (4)


Kasus # 2 : k

k

1 2 A ⎯⎯ ⎯ → B ⎯⎯ ⎯ → C

RA = k1 C A ,

RB = k1 C A − k2 C B

Secara sistematik, diagram kerja MODEL dari Kasus #2 ini adalah sbb :

Masukan CA,0, CB,0 k1, k2

Persamaan Sistem dC A = k1 C A dt

Keluaran

CA, CB

dC B = k1 C A − k2 C B dt

Model Kasus #2 Dalam mencari solusi (jawab) dari Kasus #2 ini sebagai fungsi dari waktu, baik secara analitis maupun numeris diperlukan sejumlah HARGA AWAL (initial values), sbb : 1. CA,0 2. CB,0 dan sejumlah TETAPAN (laju reaksi), sbb : 1. k1 2. k2

ÖÖÖ Property of Setijo Bismo

SOLUSI “PDB” dengan “HARGA AWAL”

Halaman (5)


2

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

2.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 ) = y0 x0 ≤ x ≤ x N 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 :

x i = x 0 + i h,

x N − x0 N

i = 0,1, 2,… , 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 Property of Setijo Bismo

Halaman (6)


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!

2.1. Metode EULER Metode yang paling sederhana dalam menterjemahkan deret Taylor diatas, yaitu dengan dengan cara “pemotongan” term kedua (atau di atasnya). Dan bila dinotasikan ui ≈ y( x i ) , maka : ui +1 = ui + h f ( x i , ui ), i = 0,1,… , N − 1

u0 = y0 Dari persamaan di atas, dapat diketahui bahwa untuk menghitung

ui +1 diperlukan informasi tentang harga-harga x i dan ui , yang disebut sebagai harga awal. Akurasi dari metode ini adalah dalam “order satu” (first-order approximation) : ei +1 = 0( h1 )


Halaman (7)


Gambar 1. Representasi METODE EULER.

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

2.3. Metode RUNGE-KUTTA Metode ini termasuk algoritma eksplisit yang melibatkan evaluasi fungsi f di antara x i dan x i +1 . Formula umum dari metode ini adalah sebagai berikut :

ui +1 = ui +

ν

∑ ωj K j

j =1

dengan :

K j = h f ( x i + c j h, ui + c1 = 0

j −1

∑ a jl K l )

l =1

Perlu dicatat, bila ν = 1 , ω = 1 , dan K1 = h f ( x i , ui ) 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, parameter-parameter ω , 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 :


Halaman (8)


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 :

∂f dy dfi ∂f = i + i = ( fx + f y f )i dx ∂x ∂y dx x = x i

Jika persamaan terakhir disubstitusikan kedalam persamaan di atasnya dengan pemotongan term di atasnya, maka akan diperoleh : ui +1

h2 = ui + h fi + ( fx + f y f )i 2!

Untuk ekspansi formula K j pada posisi ke- i , perlu dicatat bahwa harga K 1 merupakan bentuk paling sederhana, sebagai berikut :

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

K 2 = h f ( x i + c2 h, ui + a21 K1 ) Jika diketahui bahwa sembarang dua fungsi η dan φ yang lokasinya berturut-turut berdekatan dengan x i dan ui , maka akan diperoleh : f ( η, φ) ≈ f ( x i , ui ) + ( η − x i ) f ( x i , ui ) + ( φ − ui ) f y ( x i , ui )

Dengan menggunakan persamaan di atas untuk K 2 , maka akan didapatkan bentuk persamaan untuk menghitung K 2 : Property of Setijo Bismo

Halaman (9)


K 2 = h ( fi + c2 h fx + a21 K1 f y )

atau K 2 = hfi + h 2 ( c2 fx + a21 f y f )i

Bila disubstitusikan persamaan terakhir dan persamaan K1 di atas kedalam formula dasar Runge-Kutta ( ui +1 ) di atas, maka akan diperoleh : ui +1 = ui + ω1 h fi + ω2 h fi + ω2 h 2 c2 ( fx )i + a21 ω2 h 2 ( f y f )i

Jika dibandingkan persamaan ui +1 yang terakhir ini dengan persaamaan ui +1 sebelumnya, maka akan diperoleh : ω1 + ω2 = 1,0 ω2 c2 = 0,5 ω2 a21 = 0,5

Algoritma METODE RUNGE-KUTTA dilengkapi dengan cara memilih salah satu di antara paramater-parameter ω1 , ω2 , c2 atau a21 , 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 : ui +1 = ui + h f ( x i + 0,5 h, ui + 0,5 h fi ), i = 0,1,… , N − 1 u0 = y0 Property of Setijo Bismo

Halaman (10)


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

Gambar 2. Representasi METODE TITIK TENGAH.

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


Halaman (11)


ui +1 = ui +

h [fi + f ( x i + h, ui + h fi )], i = 0,1,… , N − 1 2

u0 = 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 2 pemotongan deret Taylor dilakukan setelah 0( h ) . 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 (12)


2.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 :

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


Halaman (13)


2.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 dituliskan seperti di bawah ini :

1 6

RKG

ini

dapat

1 3

ui +1 = ui + ( K1 + K 4 ) + (b K 2 + d K 3 )

K1 = h f ( x i , ui ) K 2 = h f ( x i + 12 h, ui + 12 K1 ) K 3 = h f ( x i + 12 h, ui + a K1 + b K 2 )

K 4 = h f ( x i + h, ui + 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 : u0 = y0


Halaman (14)


2. Metode RUNGE-KUTTA-MERSON : Metode Runge-Kutta-Merson (RKM) tergolong dalam keluarga metode Runge-Kutta 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 order5. Formulasi ringkas dari metode RKM ini dapat dituliskan seperti di bawah ini : 1 2

3 2

1 6

2 3

ui +1 = ui + K1 − K 3 + 2 K 4 1 6

ui +1 = ui + K1 + K 4 + K 5

K1 = h f ( x i , ui ) K 2 = h f ( x i + 13 h, ui + 13 K1 ) K 3 = h f ( x i + 13 h, ui + 16 K1 + 16 K 2 ) K 4 = h f ( x i + 12 h , ui + 18 K 1 + 38 K 3 )

K 5 = h f ( x i + h, ui + 12 K1 − 32 K 3 + 2 K 4 ) untuk : i = 0, 1, 2,…,N-1 dan harga kondisi awal : u0 = y0


Halaman (15)


3. Metode RUNGE-KUTTA-FEHLBERG : Sama halnya dengan metode RKM, metode RungeKutta-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 : K1 = h f ( x i , ui ) K 2 = h f ( x i + 14 h, ui + 14 K1 ) K 3 = h f ( x i + 38 h, ui + 323 K1 + 329 K 2 ) K 4 = h f ( x i + 12 h, ui + 1932 K1 − 7200 K 2 + 7296 K3 ) 13 2197 2197 2197 845 K 5 = h f ( x i + h, ui + 439 K1 − 8 K 2 + 3680 K 3 − 4104 K4 ) 216 513

K 6 = h f ( x i + 12 h, ui − 278 K 1 + 2 K 2 − 3544 K 3 + 1859 K 4 − 11 K5 ) 2565 4104 40

• Formula ‘update’ order-4 : ui +1 = ui +

25 216

K1 +

1408 2565

K3 +

2197 4104

1 5

K4 − K5

• Formula order-5 : uˆ i +1 = ui +

16 135

K1 +

6656 12825

K3 +

28561 56437

K4 −

9 50

K5 +

2 55

K6

• Galat ‘pembabatan’ order-4 : uˆ i +1 − ui +1 =

untuk :

1 360

K1 −

128 4275

K3 −

2197 75240

K4 +

1 50

K5 +

2 55

K6

i = 0, 1, 2,…,N-1 dan u0 = y0


Halaman (16)


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.

Hasil Perhitungan Numerik (ui+1)

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


Halaman (17)


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 (18)


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 (19)


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 (20)


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 (21)


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 (22)


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 (23)


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 (24)


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 (25)


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,' For I := 1 to 10 do Begin DRKMER(N,xv,Yi,Yf,s,Eps,MSG); Writeln(xv[1]:0:3,' ',Yf[1]:0:6,' ',MSG); xv[0] := xv[1]; xv[1] := xv[1] + 0.1; Yi := Yf; End; Readln; End.


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

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

Halaman (26)


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

Hasil Perhitungan Numerik untuk Sistem PDB : x

Contoh 6 : RK-Gill

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 (27)


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,' For I := 1 to 10 do Begin DRKF45(N,xv,Yi,Yf,s,Eps,MSG); 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 (28)


Solusi PDB Order-2 menggunakan Metode RKG dan RKM dengan ‘Teknik Shooting’ :

Persamaan Diferensial Biasa Order-2 non-linier berikut : d2 y 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 : ⎧ dy1 ⎪⎪ dx = y2 ⎨ ⎪ dy2 = − 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 y2 ,


Halaman (29)


• Harga awal (kondisi awal) dari y1 diketahui, namun untuk y2 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 . 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 (30)



Halaman (31)


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,' For I := 1 to 20 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 (32)


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 (33)


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) : ⎧ dy1 ⎪ dx = 1 ⎪ ⎪ dy2 = y3 ⎨ dx ⎪ ⎪ dy3 ⎪ dx = − sin( y1 ) + cos( y1 ) - y3 ⎩

dy1 ⎧ ⎪ y1 = x ⇒ dx = 1 ⎪ ⎨ y2 = y ⎪ dy ⎪ y3 = 2 dx ⎩

⎧ y1 (0) = 0 ⎪ ⎪⎪ ⎨ y2 (0) = 0 ⎪ ⎪ ⎪⎩ y3 (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 y2 , • Permisalan ketiga (yang terakhir) adalah untuk ‘turunan variabel terikat’ ( dy dx = dy2 dx ) sebagai y3 , • Sistem PDB baru yang diperoleh adalah dengan cara menyusun turunanturunan dari variabel-variabel permisalan di atas ( y1 dan y2 ) yang digabungkan dengan penyusunan ulang PDB tunggal di atas untuk dy3 dx (lihat kolom kiri)


Halaman (34)


• 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 y2 (0) ), dan turunan pertama variabel terikat pada saat x = 0 (dalam hal ini y3 (0) ).

• Solusi yang diinginkan adalah harga-harga y2 sebagai fungsi x .

Perbandingan harga-harga solusi numerik dengan metode RKG dari PDB di atas ( y2 ) dengan solusi eksak disajikan pada tabel di bawah ini : Tabel 7. Perbandingan solusi numerik sistem PDB order-2 dengan metode RKGill dengan solusi analitis.

x

Metode Runge-Kuta-Gill (Contoh 11) :

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


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 (35)


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 (36)

PRINSIP DASAR PEMODELAN dan MODEL MATEMATIS

Recommend Documents