PENYELESAIAN PERSAMAAN DIFERENSIAL BIASA KASUS: INITIAL KASUS: INITIAL VALUE PROBLEM (IVP) Materi Kuliah: Pengantar; Metode Euler; Perbaikan Metode Metode Euler; Metode Metode Runge-Kutta; Penyelesaian Sistem Sistem Persamaan Diferensial Diferensial Biasa secara secara Simultan
PENGANTAR Persamaan diferensial: Persamaan yang melibatkan turunan atau derivatif fungsi-fungsi Persamaan diferensial berorde n: Persamaan diferensial yang memuat turunan fungsi tertinggi berorde n Persamaan diferensial berorde satu (first order): Persamaan diferensial yang turunan fungsi tertingginya berorde 1 dy Contoh: = 2 x + 5 dx dC A k 1 C A
−
dt
=
1 + k 2 C A
y' = 10 x y 2 e 2 x Persamaan diferensial berorde dua (second order): Persamaan diferensial yang turunan fungsi tertingginya berorde 2
Contoh:
d 2 y dx 2
−
dy dx
− 6 y = e x
Beberapa penggolongan persamaan diferensial:
1. Berdasarkan banyaknya perubah bebas : a. Persamaan diferensial biasa (PDB) atau ordinary differential equation (ODE) equation (ODE) Yakni persamaan diferensial dengan perubah bebas tunggal. dC 2 Misal: − A = k C A dt b. Persamaan diferensial parsial (PDP) atau partial atau partial differential equation (PDE) equation (PDE) Yakni persamaan diferensial dengan jumlah perubah bebas lebih dari satu. Misal:
∂ T ∂ 2T = k 2 ρ Cp ∂ t ∂ z
2. Berdasarkan persoalan syarat atau nilainya : a. Persamaan diferensial dengan persoalan syarat/nilai awal (intial value problem, IVP). Yakni jika semua syarat diberikan pada satu nilai perubah bebas (yakni pada nol atau x 0) Misal:
d 2 y dx
2
= − y
dengan: y(0) = 2 dan y’(0) = -1
b. Persamaan diferensial dengan persoalan syarat/nilai batas (boundary value problem, BVP). Yakni jika syarat-syarat diberikan pada lebih dari satu nilai perubah bebas. Misal:
d 2 y dx
2
= − y
dengan: y(0) = 2 dan y’(3π/2) = 1
Yang akan dipelajari dalam materi kuliah ini: Penyelesaian atau integrasi numerik persamaan diferensial biasa (berorde satu) dengan persoalan nilai awal dy/analisis_numerik/penyelesaian_ dy/analisis_numerik/penyelesaian_pdb_kasus_ivp/ pdb_kasus_ivp/semester semester pendek_0607/july2007/halaman pendek_0607/july2007/halaman 1 dari 12
Menyelesaikan atau mengintegrasi persamaan diferensial:
dy dx
= f ( x , y )
... (1)
dengan syarat awal: y(x0) = y0 secara numerik berarti: #
menentukan atau menghitung nilai-nilai pendekatan y1, y2, y3, dst. dari penyelesaian eksak y1 , # # y2 , y3 , dst. pada x = x1, x = x2, x = x3, dst. #
#
#
(y1 , y2 , y3 , dst. sendiri biasanya justru tidak diketahui nilainya) Titik
(x0, y0) digunakan sebagai titik tolak pengintegrasian
Sebuah PDB disebut stabil , jika dalam arah integrasi, penyelesaiannya bersifat konvergen . Dan sebaliknya, PDB disebut tidak stabil, jika dalam arah integrasi, penyelesaiannya bersifat divergen . Dua jenis metode penyelesaian numerik persamaan diferensial ini:
1. Metode satu langkah (one-step methods) yang akan dipelajari dalam materi kuliah ini 2. Metode banyak langkah (multi-steps methods) One-step method yang akan dipelajari di sini: 1. Metode Euler (eksplisit) 2. Penyempurnaan atau perbaikan metode Euler (Metode Heun, Metode Titik Tengah) 3. Metode Runge-Kutta
METODE EULER Merupakan metode yang paling sederhana untuk mengintegrasikan PDB orde satu secara numerik. Kondisi atau syarat atau nilai awal (x0, y0) digunakan untuk menghitung besarnya slope (atau tangen arah) y(x) pada x = x0: dy dx x = x
= f ( x0 , y0 )
... (2)
0
Dengan menganggap bahwa slope (dy/dx) pada interval diperkirakan sebesar: y( x0
Δx
bernilai tetap, maka nilai y(x0+Δx) dapat
+ Δ x ) = y( x0 ) + Δ x f ( x0 , y0 )
Selanjutnya, nilai-nilai x dan y ini (yakni x = x0+Δx dan y = y(x0+Δx)) digunakan untuk memperkirakan besarnya slope pada titik yang baru. Atau, nilai y(x0+2Δx) dapat dihitung sbb:
+ 2Δ x ) = y( x0 + Δ x ) + Δ x f ( x0 + Δ x , y( x0 + Δ x ))
y( x0
Demikian seterusnya. Pola perhitungan yang beruntun ini digambarkan sebagai metode Euler: y( xi
+ Δ x ) = y( xi ) + Δ x . f ( xi , y( xi ))
atau:
yi + 1 = yi
atau:
yi + 1
dengan:
+ Δ x . f ( xi , yi )
... (3)
= yi + h . f ( xi , yi )
Δx = h
... (3)
menyatakan lebar langkah (step size)
f (xi,yi) merupakan bentuk persamaan diferensial seperti pada persamaan (1), sehingga: yi +1
= yi + Δ x
dy dx x , y i
atau:
yi + 1 = yi
i
+h
dy dx x , y i
... (3) i
Persamaan (3) merupakan formula metode Euler. Perhatikan bahwa formula metode Euler ini juga dapat dijabarkan dari ekspansi deret Taylor untuk yi+1 di sekitar yi: dy/analisis_numerik/penyelesaian_pdb_kasus_ivp/semester pendek_0607/july2007/halaman 2 dari 12
yi + 1 = yi
dy
+
dx x
Δ x + i
1 d 2 y 2
2 dx x i
Δ x
2
+
1 d 3 y 3
6 dx x i
Δ x 3 + ...
diabaikan dengan mengabaikan suku-suku berorde Δx (=h ) dan yang lebih tinggi. 2
2
Dengan kata lain, metode Euler ini mempunyai tingkat ketelitian yang dinyatakan dengan local truncation error sebesar: 2
ei = Ο ( Δ x )
2
atau:
ei = Ο (h )
Metode ini mempunyai global truncation error sebesar: E i = nilai eksak yi – nilai pendekatan numerik yi
CONTOH SOAL #:
Gunakan metode Euler untuk menghitung nilai y pada x = 1 jika:
dy dx
= x 2 y
dengan nilai awal: y = 1 pada x = 0
Penyelesaian: Formula metode Euler untuk kasus ini dapat dituliskan sebagai: yi + 1 = yi
(
+ Δ x xi 2 yi
)
Jika diambil step size Δx = 0,1, maka: pada x0 = 0 dan y0 = 1 dapat dihitung: 2
y1 = 1 + (0,1) (0) (1) = 1 Selanjutnya, pada x1 = x0 + Δx = 0 + 0,1 = 0,1 dan y1 = 1 dapat dihitung: 2
y2 = 1 + (0,1) (0,1) (1) = 1,001 Selanjutnya, pada x2 = x1 + Δx = 0,1 + 0,1 = 0,2 dan y2 = 1,001 dapat dihitung: 2
y3 = 1,001 + (0,1) (0,2) (1,001) = 1,005 Demikian seterusnya, hingga diperoleh y pada x = 1. Sebagai perbandingan, dapat diambil nilai step size yang lain, misalnya: Δx = 0,05, Δx = 0,02, dan Δx = 0,2. Dengan cara yang sama, maka dapat diperoleh hasil-hasil perhitungan sbb.: Nilai y
Δx
x
0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1
0,1 1,0000 1,0000 1,0010 1,0050 1,0140 1,0303 1,0560 1,0940 1,1476 1,2211 1,3200
0,05 1,0000 1,0001 1,0018 1,0069 1,0176 1,0361 1,0650 1,1070 1,1661 1,2468 1,3559
0,02 0,2 1,0000 1,0000 1,0002 1,0023 1,0000 1,0081 1,0199 1,0080 1,0400 1,0707 1,0403 1,1154 1,1718 1,1152 1,2635 1,3792 1,2579
1.5
Analitik
1,0000 1,0003 1,0027 1,0090 1,0216 1,0425 1,0747 1,1211 1,1861 1,2751 1,3956
1.5
Analitik
1.4
Δx = 0,2
1.4 1.3
Δx = 0,1
y1.3
Δx = 0,05
1.2 1.2
Δx = 0,02
1.1 1.1 1.0 0
0 .1 0.2 0 .3 0.4 0 .5 0 .6 0. 7 0.8 0 .9
1
x
(Bandingkan nilai y pada x = 1 yang diperoleh melalui perhitungan secara numerik dengan nilai Δ x yang berbeda-beda. Bandingkan juga dengan nilai y eksak (secara analitik).) Kesimpulannya: ......................................?
dy/analisis_numerik/penyelesaian_pdb_kasus_ivp/semester pendek_0607/july2007/halaman 3 dari 12
PERBAIKAN METODE EULER Sumber mendasar (pokok) terjadinya penyimpangan yang relatif besar pada penerapan metode Euler adalah karena: turunan fungsi pada awal interval yang diasumsikan tetap di sepanjang interval x. Untuk itu, dilakukan modifikasi terhadap metode Euler, antara lain: (1) Metode Heun, dan (2) Metode Titik Tengah (midpoint method).
Metode Heun Metode ini menyempurnakan metode Euler melalui penentuan dua nilai turunan fungsi sepanjang interval Δx, yakni: (a) di awal interval Δx, dan (b) di akhir interval Δx. Kedua nilai turunan ini selanjutnya dirata-ratakan untuk menghasilkan perkiraan nilai slope pada keseluruhan interval Δx. Tinjau kembali metode Euler di atas; nilai slope pada awal interval Δx: dy = y' i = f ( xi , yi ) dx x , y i
... (4)
i
digunakan untuk mengekstrapolasi linier nilai yi+1: yio+ 1 atau: yio+ 1
= yi + Δ x . f ( xi , yi )
... (5)
= yi + h . f ( xi , yi )
... (5)
Berbeda dengan metode Euler yang menjadikan bentuk pada persamaan (3) atau (5) sebagai jawaban akhir, metode Heun menjadikannya sebagai prediksi antara (intermediate prediction). (Dalam hal ini, superscript o digunakan untuk membedakannya). Persamaan (5) biasa disebut sebagai persamaan prediksi (predictor equation). Persamaan ini selanjutnya digunakan untuk memperkirakan besarnya slope pada akhir interval Δx yang ditinjau, yakni: dy dx x , y o i +1 i +1
= y' i +1 = f ( xi +1 , yio+1 )
... (6)
Slope rata-rata yang dihitung berdasarkan persamaan (4) dan (6) adalah: dy dx
= y' =
f ( xi , yi ) + f ( xi +1 , yio+1 )
... (7)
2
Slope rata-rata pada persamaan (7) ini selanjutnya digunakan untuk mengekstrapolasi linier dari yi ke yi+1 menggunakan metode Euler: yi +1 = yi atau: yi +1
+
= yi +
f ( xi , yi ) + f ( xi +1 , yio+1 ) 2 f ( xi , yi ) + f ( xi +1 , yio+1 ) 2
h
... (8)
Δ x
... (8)
Persamaan (8) biasa disebut sebagai persamaan koreksi (corrector equation). Metode Heun menggunakan pendekatan predictor-corrector , yang secara iteratif dapat dinyatakan sebagai berikut: Predictor :
yio+ 1
= yim + Δ x . f ( xi , yi )
atau: yio+1 Corrector :
j yi + 1
= yim
j atau: yi +1
... (9)
= yim + h . f ( xi , yi ) +
... (9)
j − 1 m f ( xi , yi ) + f ( xi +1 , yi + 1 )
= yim
2
+
h
j − 1 f ( xi , yim ) + f ( xi + 1 , yi + 1 )
2
... (10)
Δ x
... (10)
(untuk j = 1, 2, ..., m) (j menyatakan nomor langkah iterasi) dy/analisis_numerik/penyelesaian_pdb_kasus_ivp/semester pendek_0607/july2007/halaman 4 dari 12
(Perhatikanlah bahwa pada persamaan (10), yi+1 ada pada kedua ruas persamaan, sehingga perhitungan yi+1 dapat dilakukan secara berulang (atau iteratif) agar diperoleh hasil perkiraan yi+1 yang lebih baik.)
Metode Titik Tengah (Midpoint) Metode ini menggunakan metode Euler untuk memperkirakan sebuah nilai y pada titik tengah interval Δx yang ditinjau, yakni sebesar: yi 1 + 2
= yi + f ( xi , yi )
atau: yi 1 + 2
= yi + f ( xi , yi )
h 2 Δ x
... (11)
2
Persamaan (11) selanjutnya digunakan untuk memperkirakan nilai slope pada titik tengah interval: dy dx x , y i+ 1 i+ 1 2
= y' i + 1 = f ( xi + 1 2
, y 2
i+ 1
)
... (12)
2
2
yang dianggap dapat mewakili slope rata-rata pada keseluruhan interval Δx. Nilai slope pada persamaan (12) ini selanjutnya digunakan untuk mengekstrapolasi linier dari yi ke yi+1: yi +1 atau: yi + 1
= yi + f ( xi + 1
... (13)
, yi + 1 ) . h
2
2
= yi + f ( xi + 1 , yi + 1 ) . Δ x 2 2
... (13)
METODE RUNGE-KUTTA Merupakan metode yang paling banyak diterapkan untuk integrasi numerik persamaan diferensial biasa dengan initial value problem, karena menghasilkan pendekatan yang cukup baik. Metode ini menggunakan pendekatan deret Taylor yang cukup akurat, tanpa membutuhkan perhitungan turunan yang lebih tinggi. Bentuk umum metode-metode Runge Kutta:
yi +1
= yi + φ h
... (14)
φ biasa
disebut sebagai fungsi inkremen , yang dapat dianggap sebagai nilai slope pada keseluruhan interval h atau Δx yang ditinjau. Fungsi inkremen (φ) mempunyai bentuk umum: φ = a1 k 1 + a2 k 2
+ ... + an
k n
... (15)
a merupakan konstanta, dan k dapat dinyatakan sebagai: k 1 = f ( xi , yi ) k 2
= f ( xi + p1 h , yi + q11 k 1 h )
k 3
= f ( xi + p2
h , yi
+ q21 k 1 h + q22
k 2 h )
k n
= f ( xi + pn −1 h , yi + qn −1 ,1 k 1 h + qn −1 ,2
k 2 h + ... + qn −1 ,n −1 k n −1 h )
... (16)
p dan q merupakan konstanta. Parameter-parameter a, p, dan q dipilih sedemikian sehingga 2 perumusannya sesuai dengan ekspansi deret Taylor sampai dengan suku yang melibatkan faktor h 2 (atau (Δx) ). Perhatikan bahwa: Metode Euler merupakan salah satu jenis metode Runge-Kutta yang berorde satu (atau n = 1).
dy/analisis_numerik/penyelesaian_pdb_kasus_ivp/semester pendek_0607/july2007/halaman 5 dari 12
Metode Runge-Kutta Orde 4 Klasik Metode Runge-Kutta yang paling umum digunakan adalah metode Runge-Kutta berorde 4. Bentuk umumnya: yi +1
= yi +
1 6
( k 1 + 2 k 2
+ 2 k 3 + k 4 ) h
... (17)
dengan: k 1 = f ( xi , yi ) k 2
= f ( xi +
1
... (18)
+
1
k 1 h )
... (19)
h , yi + k 2 h ) 2 2 k 4 = f ( xi + h , yi + k 3 h )
... (20)
k 3
= f ( xi +
2 1
h , yi
2 1
... (21)
Perhatikan bahwa: Jika dy/dx atau f hanya merupakan fungsi x saja , maka metode Runge-Kutta orde 4 ini sama dengan integrasi numerik dengan metode Simpson 1/3. CONTOH SOAL #:
Lihat kembali contoh soal sebelumnya (pada Metode Euler). Gunakan metode Runge-Kutta orde 4 untuk menghitung nilai y pada x = 1 jika:
dy
= x 2 y
dx dengan nilai awal: y = 1 pada x = 0. Bandingkan hasilnya dengan perhitungan menggunakan metode Euler.
Penyelesaian: Formula metode Runge-Kutta untuk kasus ini: Δ x yi + 1 = yi + (k 1 + 2 k 2 + 2 k 3 + k 4 ) 6 dengan:
k 1 ,i
= xi 2 yi 2
Δ x ⎞ ⎛ Δ x ⎛ ⎞ k 2 ,i = ⎜ xi + ⎟ ⎜ yi + k 1 ,i ⎟ 2 ⎠ ⎝ 2 ⎝ ⎠ 2
Δ x ⎞ ⎛ Δ x ⎛ ⎞ k 3 ,i = ⎜ xi + ⎟ ⎜ yi + k 2 ,i ⎟ 2 ⎠ ⎝ 2 ⎝ ⎠ k 4 ,i
= ( xi + Δ x )2 ( yi + Δ x k 3 ,i )
Jika diambil step size Δx = 0,1, maka pada x0 = 0 dan y0 = 1 dapat dihitung: k 1 ,0
= x0 2 y0 = 0 2 1 = 0
k 2 ,0
0 ,1 ⎞ ⎛ 0 ,1 ⎞ = ⎛ ⎜ x0 + ⎟ ⎜ y0 + k 1 ,0 ⎟ = 0 ,05 2 1 = 0 ,0025 2 ⎠ ⎝ 2 ⎝ ⎠
2
2
0 ,1 ⎞ ⎛ 0 ,1 ⎛ ⎞ k 3 ,0 = ⎜ x0 + ⎟ ⎜ y0 + k 2 ,0 ⎟ = 0 ,05 2 1 ,000125 = 0 ,0025 2 ⎠ ⎝ 2 ⎝ ⎠ k 4 ,0 = ( x0 + 0 ,1)2 ( y0 + 0 ,1 k 3 ,0 ) = 0 ,12 1 ,000250 = 0 ,01000 sehingga: y1 = y( 0 ,1 ) = 1 +
0 ,1
(0 + 2 . 0 ,0025 + 2 . 0 ,0025 + 0 ,01) 6 y1 = y(0,1) = 1,000333
Demikian seterusnya, hingga diperoleh y pada x = 1. Sebagai perbandingan, dapat diambil nilai step size yang lain, misalnya: Δx = 0,2 dan Δx = 0,05. Dengan cara yang sama, maka dapat diperoleh hasil-hasil perhitungan sbb.: dy/analisis_numerik/penyelesaian_pdb_kasus_ivp/semester pendek_0607/july2007/halaman 6 dari 12
x x = 0,2 Euler RK-4 1,0000 1,000000
0 0,1 0,2 1,0000 0,3 0,4 1,0080 0,5 0,6 1,0403 0,7 0,8 1,1152 0,9 1 1,2579
1,002670 1,021562 1,074655 1,186094 1,395608
Nilai y x = 0,1 Euler RK-4 1,0000 1,000000 1,0000 1,000333 1,0010 1,002670 1,0050 1,009041 1,0140 1,021563 1,0303 1,042547 1,0560 1,074655 1,0940 1,121126 1,1476 1,186095 1,2211 1,275069 1,3200 1,395613
x = 0,05 Euler RK-4 1,0000 1,000000 1,0001 1,000333 1,0018 1,002670 1,0069 1,009041 1,0176 1,021563 1,0361 1,042547 1,0650 1,074655 1,1070 1,121126 1,1661 1,186095 1,2468 1,275069 1,3559 1,395612
Analitik
1,000000 1,000333 1,002670 1,009041 1,021563 1,042547 1,074655 1,121126 1,186095 1,275069 1,395612
Kesimpulannya: ......................................?
PENYELESAIAN SISTEM PERSAMAAN SECARA SIMULTAN Dalam prakteknya, persoalan-persoalan sains dan rekayasa (engineering) seringkali melibatkan penyelesaian sistem (atau sekumpulan) persamaan diferensial biasa secara simultan. Bentuk umum sistem persamaan sejumlah n dapat dinyatakan sebagai:
dy1 dx dy2 dx
= f 1 ( x , y1 , y2 ,..., yn ) = f 2 ( x , y1 , y2 ,..., yn )
dyn dx
= f n ( x , y1 , y2 ,..., yn )
... (22)
dengan sejumlah n nilai awal: x = x0; y1(x0) = y10; y2(x0) = y20; ..., yn(x0) = yn0 Persamaan (22) dapat diselesaikan dengan metode-metode yang sudah dipelajari sebelumnya (misal: metode Euler dan metode Runge-Kutta (RK) orde 4 ) melalui langkah perhitungan yang sangat identik dengan perhitungan untuk kasus persamaan diferensial tunggal.
Contoh Ilustratif: Penyelesaian sistem 2 buah persamaan diferensial biasa orde satu secara simultan dengan metode Runge-Kutta orde 4 Bentuk persamaan diferensial: dy dx
= f 1 ( x , y , z )
dan
dz dx
= f 2 ( x , y , z )
dengan 2 nilai awal: x = x0; y = y0; z = z0 Formula Runge-Kutta Orde 4 untuk menentukan xi+1, yi+1, dan zi+1 berdasarkan xi, yi, dan zi: xi +1 = xi
+ Δ x
yi + 1 = yi
+ Δ6 x (k 1 + 2.k 2 + 2.k 3 + k 4 )
z i + 1 = z i
+ Δ6 x (l 1 + 2.l 2 + 2.l 3 + l 4 )
dengan:
k 1 = f 1 ( xi , yi , z i ) l 1 = f 2 ( xi , yi , z i )
dy/analisis_numerik/penyelesaian_pdb_kasus_ivp/semester pendek_0607/july2007/halaman 7 dari 12
Δ x k l ⎞ = f 1 ⎛ ⎜ xi + , yi + 1 Δ x , z i + 1 Δ x ⎟ 2 2 2 ⎠ ⎝ k l ⎞ ⎛ Δ x l 2 = f 2 ⎜ xi + , yi + 1 Δ x , z i + 1 Δ x ⎟ 2 2 2 ⎠ ⎝ k 2
k l Δ x ⎞ = f 1 ⎛ ⎜ xi + , yi + 2 Δ x , z i + 2 Δ x ⎟ 2 2 2 ⎠ ⎝ k l ⎛ Δ x ⎞ l 3 = f 2 ⎜ xi + , yi + 2 Δ x , z i + 2 Δ x ⎟ 2 2 2 ⎠ ⎝ k 3
k 4 l 4
= f 1 ( xi + Δ x , yi + k 3 Δ x , z i + l 3 Δ x ) = f 2 ( xi + Δ x , yi + k 3 Δ x , z i + l 3 Δ x )
Penyelesaian Persamaan Diferensial Biasa Berorde Tinggi (n) Sebagai contoh ilustratif, tinjaulah sebuah persamaan diferensial biasa berorde 2: d 2 y dx 2
+ A( x )
dengan 2 nilai awal: Ambil pemisalan: sehingga:
dy dz
+ B( x ) y + C ( x ) = 0
y(x0) = a z = dz dx
dy
dan
dx x
... (23)
=b 0
dy
... (24)
dx
=
d ⎛ dy ⎞
d 2 y
dx ⎝ dx ⎠
2
⎜
⎟=
dx
... (25)
Substitusikan (24) dan (25) ke persamaan (23): dz
atau:
dx dy
Persamaan (24) dapat dituliskan sebagai:
dx
dz dx
+ A( x ) z + B( x ) y + C ( x ) = 0
= − A( x ) z − B( x ) y − C ( x )
... (26)
= z
... (27)
Berdasarkan 2 persamaan terakhir (yakni (26) dan (27)), terlihat bahwa persamaan diferensial biasa berorde 2 pada persamaan (23) dapat diubah menjadi 2 buah persamaan diferensial biasa berorde 1 yang dapat diselesaikan secara simultan. Dua nilai awal sistem PD ini sekarang berubah menjadi: y(x0) = a dan z(x0) = b Hal yang sama/identik dapat diterapkan untuk menyelesaikan persamaan diferensial biasa berorde lebih tinggi. Secara umum: PDB berorde n dapat diubah menjadi n buah PDB berorde 1, yang selanjutnya dapat diselesaikan secara simultan.
Atau: Untuk sebuah PDB berorde n:
dengan n buah nilai awal :
⎛ d n y d n −1 y d n − 2 y ⎞ dy ⎜ ⎟ F ⎜ dx n , dx n −1 , dx n − 2 , ... , dx , y , x ⎟ = 0 ⎝ ⎠
d n −1 y dx n −1
= an −1
d n − 2 y dx n − 2
= an − 2
Pada x = x0
dy
= a1 dx y = a0 dy/analisis_numerik/penyelesaian_pdb_kasus_ivp/semester pendek_0607/july2007/halaman 8 dari 12
Lakukan beberapa substitusi berikut:
z 1 = z 2
=
dy dx dz 1 dx
=
d 2 y dx 2
z n −1
=
dz n − 2 dx
=
d n −1 y dx n −1
Dan setelah dilakukan proses substitusi ke dalam persamaan diferensial, maka akan diperoleh: dz n −1 dx dz n − 2 dx
= G ( z n −1 , z n − 2 , ... , z 1 , y , x ) = z n −1 Merupakan n buah PDB orde 1 simultan
dz 1 dx dy dx
= z 2 = z 1
dengan n buah nilai awal: z n −1 = an −1 z n − 2 = an − 2
Pada x = x0
= a1 y = a0
z 1
CONTOH APLIKASI PROFIL KONSENTRASI DAN SUHU SEPANJANG WAKTU PADA SISTEM REAKTOR BATCH NON-ISOTERMAL (ADIABATIK)
Sebuah reaktor sistem batch non-isotermal beroperasi secara adiabatik. Di dalam reaktor berlangsung reaksi homogen fase cair: A P dengan kecepatan reaksi sebesar: r = k CA
⎛ E ⎞ ⎟⎟ R T ⎝ ⎠
dan: k = k 0 exp ⎜⎜ −
CA ≡ konsentrasi A; E ≡ energi aktivasi reaksi; R ≡ konstanta gas; dan T ≡ suhu mutlak. insulation
Karena sistem reaksi dianggap teraduk sempurna, maka neraca mol A pada unsteady state dapat dituliskan sebagai: dn A dt
= V R ( −r ) = −V R
⎛ ⎝
k 0 exp ⎜ −
E ⎞
⎟ C A
RT ⎠
Karena volume reaktor, VR , konstan, dan: C A maka:
dC A dt
= −k 0
⎛ ⎝
C A exp ⎜ −
E ⎞
⎟
RT ⎠
=
n A V R
,
... (*)
Neraca panas pada unsteady state dapat dituliskan sebagai: dT = −Δ H R r V R ρ V R C p dt ρ ≡ densitas campuran reaksi; C p ≡ kapasitas panas campuran reaksi rata-rata; dan ΔHR ≡ panas reaksi. Maka: dy/analisis_numerik/penyelesaian_pdb_kasus_ivp/semester pendek_0607/july2007/halaman 9 dari 12
dT dt
=
− Δ H R
k 0 C A V R
ρ V R C p
⎛ ⎝
exp ⎜ −
E ⎞
⎟
... (**)
RT ⎠
Jika ρ, Cp, dan ΔHR dianggap tidak terlalu dipengaruhi oleh suhu, maka persamaan (*) dan (**) dapat ditulis secara ringkas sbb.: E ⎞ = k 1 C A exp ⎛ ⎜− ⎟ dt ⎝ RT ⎠ dT E ⎞ = k 2 C A exp ⎛ ⎜− ⎟ dt ⎝ RT ⎠
dC A
Sistem 2 buah persamaan diferensial biasa berorde 1 bernilai awal, simultan
dengan nilai awal: T = T0 dan CA = CA0 pada t = 0.
SOAL-SOAL LATIHAN 1. Tinjaulah persamaan diferensial:
dy dx
= y e 3 x
dengan: y (0) = 1,0 Dengan menggunakan step size h = 0,1, tentukan nilai y (0,3) menggunakan: a. Metode Euler b. Metode Runge-Kutta orde 4 Tunjukkan semua langkah perhitungan yang Anda lakukan. Bandingkan hasilnya dengan hasil perhitungan secara analitik. 2. Reaksi fase gas homogen: A
2 P berlangsung dalam sebuah reaktor batch isotermal pada
2
tekanan tetap, dengan: r = 0,1 CA [=] gmol/liter.detik. Mula-mula reaktor berisi 0,01 gmol A dan 0,01 gmol gas inert dengan volume 0,5 liter. Tentukan volume reaktor setelah reaksi berlangsung 25 detik. Neraca mol A pada unsteady state dinyatakan dn A
0 ,1 n A
2
= V ( −r ) = − dt V Gas dianggap sebagai gas ideal, sehingga: sebagai:
V = V 0
⎛ nt ⎞ 0 ,01 + n A + 2 ( 0 ,01 − n A ) ⎞ ⎜⎜ ⎟⎟ = 0 ,5 ⎛ ⎜ ⎟ = 0 ,75 − 25 n A [=] liter n 0 , 02 ⎝ ⎠ ⎝ t 0 ⎠
Dengan demikian:
dn A dt
=−
0 ,1 n A2 0 ,75 − 25 n A
... (*)
dengan syarat awal: nA = 0,01 pada t = 0 Petunjuk: Integrasikan persamaan (*) secara numerik untuk menentukan nA pada t = 25 detik. Selanjutnya gunakan hasil yang diperoleh untuk menghitung volume reaktor.
3. Diketahui sistem persamaan diferensial: ⎛ y ⎞ 10 y1 y2 exp ⎜ 1 ⎟ dy1 ⎝ 100 ⎠ dan = 2 dx 1 + x
y ⎞ − y1 y2 exp ⎛ ⎜ 1 ⎟ dy2 ⎝ 100 ⎠ = dx 1+ 2
dengan: y1 = y2 = 1,0 pada x = 0 Tentukan y1 dan y2 pada x = 1,0 menggunakan metode Runge-Kutta orde 4. Pilihlah 2 step size yang berbeda. Bandingkan hasil yang diperoleh. 4. Selesaikan PD berikut dari t = 0 hingga t = 2, dengan: y (0) = 1: dy = y t 3 − 1 ,5 y dt a. Secara analitik b. Menggunakan metode Euler, dengan h = 0,5 dan 0,25 dy/analisis_numerik/penyelesaian_pdb_kasus_ivp/semester pendek_0607/july2007/halaman 10 dari 12
c. Menggunakan metode titik tengah, dengan h = 0,5 d. Menggunakan metode Runge-Kutta orde 4, dengan h = 0,5 Tunjukkan semua langkah perhitungan yang Anda lakukan. 5. Selesaikan sistem PD simultan berikut: dy dt
= −2 y + 5 z e − t
dz
dan
dt
=−
y z 2 2
dengan nilai awal: y(0) = 2 dan z(0) = 4 Lakukan perhitungan dari t = 0 hingga t = 0,4, dengan step size h = 0,1, menggunakan: a. metode Euler b. metode Runge-Kutta orde 4 Plotkan hasil perhitungan Anda dalam bentuk grafik. 6. Persamaan van der Pol yang merupakan salah satu model rangkaian listrik vacuum tubes d 2 y
dinyatakan sebagai:
dx
2
− ( 1 − y 2 )
dy dx
+ y = 0
Dengan kondisi awal: y(0) = y’(0) = 1, selesaikan persamaan ini dari x = 0 hingga x = 10 menggunakan metode Euler, dengan step size sebesar: (a) 0,2, dan (b) 0,1. Plotkan hasil perhitungan yang Anda peroleh dalam sebuah grafik. d 2 y
+ 9 y = 0 dx 2 dengan step size sebesar 0,1, dari x = 0 hingga x = 4, menggunakan:
7. Selesaikan persamaan diferensial:
a. Metode Euler b. Metode Runge-Kutta orde 4 Plotkan hasil perhitungan yang Anda peroleh dalam sebuah grafik. Bandingkan juga dengan penyelesaian eksak PD ini. 8. Reaktor Semi-Batch Tinjaulah sebuah reaktor semi batch berikut ini: Reaksi fase cair yang terjadi adalah: A P
Q0, CA0
dengan:
2
r = k C A
Mula-mula reaktor diisi dengan cairan inert dengan volume V0. Pada t = 0, cairan yang mengandung A dengan konsentrasi CA0 diumpankan ke dalam reaktor dengan laju alir volumetrik Q0. VR (t)
Neraca mol A pada unsteady state: Q0 C A0
Karena: C A
=
n A V R
, maka:
dn A dt
= Q0 C A0 −
− k C A2 V R =
k n A
dn A dt
2
... (*)
V R
Cairan ditambahkan ke dalam reaktor, sehingga volume reaktor (VR ) akan bertambah sepanjang d waktu. Neraca massa keseluruhan di dalam reaktor: ( ρ V R ) = Q0 ρ dt dV R Jika ρ dianggap tetap, maka: = Q0 dt dan diintegralkan menjadi: V R = Q0 t + V 0 ... (**) Substitusikan (**) ke (*), sehingga diperoleh:
dn A dt
= Q0 C A0 −
k n A 2 Q0 t + V 0
dengan: nA = 0 pada t = 0. dy/analisis_numerik/penyelesaian_pdb_kasus_ivp/semester pendek_0607/july2007/halaman 11 dari 12
Pertanyaan:
Gunakan integrasi numerik untuk mengetahui perilaku reaktor ini hingga t = 100 detik. Diketahui: CA0 = 1,0 gmol/liter; k = 0,1 liter/gmol.detik; Q0 = 10 liter/detik; dan V0 = 50 liter 9. Aliran Cairan Antara Dua Tangki: Dua tangki silinder tegak terbuka A dan B yang masingmasing berdiameter D dan tinggi H, diletakkan sama tinggi. Bagian dasar kedua tangki dihubungkan dengan pipa horizontal berdiameter Dp yang dilengkapi dengan kran. Volume pipa dapat diabaikan terhadap volume tangki. Kran mula-mula ditutup, tangki A berisi penuh cairan, sedangkan tangki B kosong. Mulai suatu saat kran dibuka, sehingga cairan mengalir dari tangki A ke B. Kecepatan aliran cairan (υ, m/s) tergantung beda tekanan pada ujung-ujung pipa (ΔP), sesuai persamaan: υ = k Δ P dengan: k ≡ tetapan. Bagaimanakah profil tinggi permukaan cairan pada tangki A (x) dan pada tangki B (y) pada berbagai waktu (t)...? Penggambaran proses:
Tangki A
Tangki B
x
h
y M Q
N
Beda tekanan pada ujung-ujung pipa: Δ P = P M − P N = P ud + ρ g x
− P ud + ρ g y = ρ g x − y
Kecepatan aliran cairan: υ = k Δ P = k ρ g x − y 2
Debit aliran:
π k D p π Q = D p2 υ = 4 4
(
ρ g x − y
)
Neraca massa cairan di tangki A:
⎛ Rate of ⎞ ⎛ Rate of ⎞ ⎛ Rate of ⎞ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ Input ⎟ − ⎜ Output ⎟ = ⎜ Accumulation ⎟ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ 0− dx dt Neraca massa cairan total:
π k D p2
=−
(
)
ρ g x − y . ρ =
4 k D p2 D 2
(
ρ g x − y
π 2 dx D ρ 4 dt
)
…. (*)
π 2 π π D h ρ = D 2 x ρ + D 2 y ρ 4 4 4
Tinggi permukaan cairan pada tangki B:
y = h − x
Dengan demikian, persamaan (*) dapat diubah menjadi: dx dt
=−
k D p2 D
2
(
ρ g 2 x − h
)
Keadaan batas: t = 0; x = h (Besarnya h dapat Anda simulasi sendiri...!) Misal, diambil: D = 2 m; Dp = 0,02 m;
ρ = 1000 kg/m ; 3
2
g = 10 m/s ; k = 0,4
m3 kg
--- Selamat Belajar dan Berlatih Mengerjakan Soal...!!! --dy/analisis_numerik/penyelesaian_pdb_kasus_ivp/semester pendek_0607/july2007/halaman 12 dari 12
