...

persamaan diferensial biasa

by user

on
Category: Documents
0

views

Report

Comments

Transcript

persamaan diferensial biasa
http://istiarto.staff.ugm.ac.id
PERSAMAAN DIFERENSIAL BIASA
Ordinary Differential Equations – ODE
Persamaan Diferensial Biasa
2
http://istiarto.staff.ugm.ac.id
q 
Acuan
q 
Chapra, S.C., Canale R.P., 1990, Numerical Methods for Engineers, 2nd
Ed., McGraw-Hill Book Co., New York.
n 
Chapter 19 dan 20, hlm. 576-640.
Persamaan Diferensial
3
http://istiarto.staff.ugm.ac.id
FU = −c v
q 
Benda bermassa m jatuh bebas dengan
kecepatan v
F FD + FU
=
m
m
dv
c
=g− v
dt
m
a=
Hukum Newton II
c = drag coefficient (kg/s)
g = gravitational acceleration (m/s2)
persamaan diferensial
FD = m g
dv
suku diferensial laju perubahan (rate of change)
dt
Persamaan Diferensial
4
http://istiarto.staff.ugm.ac.id
FU = −c v
q 
Sebuah benda jatuh bebas
q 
Jika pada saat awal benda dalam keadaan diam:
v (t = 0) = 0
syarat awal (initial condition)
dv
c
=g− v
dt
m
FD = m g
v (t ) =
gm
1 − e −(c m ) t
c
[
]
Persamaan Diferensial
5
http://istiarto.staff.ugm.ac.id
dv
c
=g− v
dt
m
q 
Persamaan diferensial
q 
q 
dv
c
=g− v
dt
m
q 
Persamaan diferensial biasa
(ordinary differential equations, ODE)
q 
∂C
∂ 2C
−D 2 = 0
∂t
∂x
q 
v variabel tak bebas (dependent variable)
t variabel bebas (independent variable)
hanya terdiri dari satu variabel bebas
Persamaan diferensial parsial
(partial differential equations, PDE)
q 
terdiri dari dua atau lebih variabel bebas
Persamaan Diferensial
6
http://istiarto.staff.ugm.ac.id
q 
Persamaan diferensial
q 
dv
c
=g− v
dt
m
q 
d2 x
dx
m 2 +c
+ kx = 0
dt
dt
q 
tingkat (order) tertinggi suku derivatif
Persamaan diferensial tingkat-1
q 
suku derivatif bertingkat-1
Persamaan diferensial tingkat-2
q 
suku derivatif bertingkat-2
Persamaan Diferensial Biasa
7
http://istiarto.staff.ugm.ac.id
q 
Beberapa contoh ODE di bidang engineering
q 
Hukum Newton II ttg gerak
dv F
=
dt m
q 
Hukum Fourier ttg panas
Heat flux = −k
q 
Hukum Fick ttg difusi
Mass flux = −D
dT
dx
dC
dx
Persamaan Diferensial Biasa
8
http://istiarto.staff.ugm.ac.id
diketahui:
fungsi polinomial tingkat 4
y = −0.5x 4 + 4 x 3 − 10 x 2 + 8.5x + 1
diperoleh: ODE
di-diferensial-kan
dy
= −2x 3 + 12 x 2 − 20 x + 8.5
dx
Persamaan Diferensial Biasa
9
http://istiarto.staff.ugm.ac.id
5
4
3
8
2
y = −0.5x + 4x −10x + 8.5x +1
4
Y
dy
= −2x3 +12x2 − 20x + 8.5
dx
4
3
dy
0
dx
2
-4
1
-8
0
0
1
2
X
3
4
0
1
2
X
3
4
Persamaan Diferensial Biasa
10
http://istiarto.staff.ugm.ac.id
diketahui: ODE
dy
= −2x 3 + 12 x 2 − 20 x + 8.5
dx
fungsi asal
di-integral-kan
y = ∫ (− 2 x 3 + 12 x 2 − 20 x + 8.5)d x
y = −0.5x 4 + 4 x 3 − 10 x 2 + 8.5x + C
C disebut konstanta integrasi
Persamaan Diferensial Biasa
11
http://istiarto.staff.ugm.ac.id
8
dy
= −2x 3 + 12 x 2 − 20 x + 8.5
dx
4
dy
dx
8
y = −0.5x 4 + 4 x 3 − 10 x 2 + 8.5x + C
C=3
4
2
1
0
−1
−2
Y
0
0
-4
-4
-8
0
1
2
X
3
4
0
1
2
X
3
4
Persamaan Diferensial Biasa
12
http://istiarto.staff.ugm.ac.id
y = −0.5x 4 + 4 x 3 − 10 x 2 + 8.5x + C
q 
q 
q 
Hasil dari integrasi adalah sejumlah tak berhingga polinomial.
Penyelesaian yang unique (tunggal, satu-satunya) diperoleh dengan
menerapkan suatu syarat, yaitu pada titik awal x = 0, y = 1 à ini disebut
dengan istilah syarat awal (initial condition).
Syarat awal tersebut menghasilkan C = 1.
y = −0.5x 4 + 4 x 3 − 10 x 2 + 8.5x + 1
Persamaan Diferensial Biasa
13
http://istiarto.staff.ugm.ac.id
q 
Syarat awal (initial condition)
q 
q 
q 
mencerminkan keadaan sebenarnya, memiliki arti fisik
pada persamaan diferensial tingkat n, maka dibutuhkan sejumlah n
syarat awal
Syarat batas (boundary conditions)
q 
syarat yang harus dipenuhi tidak hanya di satu titik di awal saja, namun
juga di titik-titik lain atau di beberapa nilai variabel bebas yang lain
Persamaan Diferensial Biasa
14
http://istiarto.staff.ugm.ac.id
q 
Metode penyelesaian ODE
q 
q 
q 
q 
Metode Euler
Metode Heun
Metode Euler Modifikasi (Metode Poligon)
Metode Runge-Kutta
http://istiarto.staff.ugm.ac.id
15
Penyelesaian ODE
Metode Euler
http://istiarto.staff.ugm.ac.id
Metode Satu Langkah
16
16
http://istiarto.staff.ugm.ac.id
q 
y
q 
Dikenal pula sebagai metode satu
langkah (one-step method)
Persamaan:
q 
yi+1 = yi + φ h
q 
Dalam bahasa matematika:
yi+1 = yi + φ h
slope = φ
q 
xi
xi+1
step size = h
new value = old value + slope x step size
x
jadi, slope atau gradien φ dipakai untuk
meng-ekstrapolasi-kan nilai lama yi ke
nilai baru yi+1 dalam selang h
Metode Satu Langkah
17
http://istiarto.staff.ugm.ac.id
yi+1 = yi + φ h
y
q 
q 
yi+1 = yi + φ h
slope = φ
xi
xi+1
step size = h
x
q 
Semua metode satu langkah dapat
dinyatakan dalam persamaan tsb.
Perbedaan antara satu metode dengan
metode yang lain dalam metode satu
langkah ini adalah perbedaan dalam
menetapkan atau memperkirakan slope φ.
Salah satu metode satu langkah adalah
Metode Euler.
Metode Euler
18
http://istiarto.staff.ugm.ac.id
q 
q 
q 
Dalam Metode Euler, slope di xi diperkirakan
dengan derivatif pertama di titik (xi,yi).
Metode Euler dikenal pula dengan nama
Metode Euler-Cauchy.
Jadi nilai y baru diperkirakan berdasarkan
slope, sama dengan derivatif pertama di titik
x, untuk mengekstrapolasikan nilai y lama
secara linear dalam selang h ke nilai y baru.
φ = f (x i , y i ) =
dy
d x x i ,y i
yi+1 = yi + f (xi , yi )h
Metode Euler
19
http://istiarto.staff.ugm.ac.id
q 
Pakailah Metode Euler untuk mengintegralkan ODE di bawah ini, dari x = 0
s.d. x = 4 dengan selang langkah h = 0.5:
f (x , y ) =
q 
q 
dy
= −2 x 3 + 12 x 2 − 20 x + 8.5
dx
Syarat awal yang diterapkan pada ODE tsb adalah bahwa di titik x= 0,
y=1
Ingat, penyelesaian eksak ODE di atas adalah:
y = −0.5x 4 + 4 x 3 − 10 x 2 + 8.5x + 1
Metode Euler
20
http://istiarto.staff.ugm.ac.id
q 
Selang ke-1, dari x0 = 0 s.d. x1 = x0 + h = 0.5:
f (x0 , y0 ) = −2(0 )3 + 12(0 )2 − 20(0 ) + 8.5 = 8.5
y1 = y0 + f (x0 , y0 )h
= 1 + 8.5 (0.5) = 5.25
q 
Nilai y1 sesungguhnya dari penyelesaian eksak:
y1 = −0.5(0.5)4 + 4(0.5)3 − 10(0.5)2 + 8.5(0.5) + 1
= 3.21875
q 
Error, yaitu selisih antara nilai y1 sesungguhnya dan estimasi:
Et = 3.21875 − 5.25 = −2.03125 atau εt = − 2.03125 3.21875 = 63%
Metode Euler
21
http://istiarto.staff.ugm.ac.id
i
xi
yi(eksak)
yi(Euler)
0
0
1
1
1
0.5
3.21875
2
1
3
εt
8
5.25
-63%
6
3
5.875
-96%
1.5
2.21875
5.125
-131%
4
2
2
4.5
-125%
5
2.5
2.71875
4.75
-75%
6
3
4
5.875
-47%
7
3.5
4.71875
7.125
-51%
8
4
3
7
-133%
Y
Euler
4
Eksak
2
X
0
0
1
2
3
4
Metode Euler
22
http://istiarto.staff.ugm.ac.id
q 
Error atau kesalahan terdiri dari dua aspek
q 
Truncation or discretization errors (kesalahan pemotongan) yang
disebabkan oleh teknik penyelesaian dalam mengestimasikan nilai y.
n 
n 
q 
local truncation error, yaitu kesalahan pada satu langkah
propagated truncation error, yaitu kesalahan-kesalahan pada langkahlangkah terdahulu
Round-off errors yang disebabkan oleh keterbatasan jumlah digit dalam
hitungan atau jumlah digit dalam alat hitung (kalkulator, komputer).
Metode Euler
23
http://istiarto.staff.ugm.ac.id
yʹ′ = f (x , y ) =
q 
dy
dx
Deret Taylor
yʹ′ʹ′ 2
y (n ) n
yi +1 = yi + yʹ′ h + h + ... +
h + Rn
2
n!
y (n+1) (ξ) n+1
h = xi +1 − xi , Rn =
h
(n + 1)!
ξ adalah sembarang titik di antara xi dan xi+1.
Deret Taylor dapat pula dituliskan dalam bentuk lain sbb.
f ʹ′(xi , yi ) 2
f (n−1 ) (xi , yi ) n
yi +1 = yi + f (xi , yi )h +
h + ... +
h + O(hn+1 )
2
n!
O(hn+1) menyatakan bahwa local truncation error adalah proporsional terhadap selang
jarak dipangkatkan (n+1).
Metode Euler
24
http://istiarto.staff.ugm.ac.id
f ʹ′(xi , yi ) 2
f (n−1 ) (xi , yi ) n
yi +1 = yi + f (xi , yi )h +
h + ... +
h + O(hn+1 )
2
n!
Euler
Error, Et
q 
q 
q 
true local truncation error of the Euler Method (Et)
untuk selang h kecil, error mengecil seiring dengan
peningkatan tingkat
error Et dapat didekati dengan Ea
Et ≈ E a =
f ʹ′(xi , yi ) 2
h
2
atau
Ea = O(h2 )
Metode Euler
25
http://istiarto.staff.ugm.ac.id
f (x , y ) = −2x 3 + 12 x 2 − 20 x + 8.5
q 
q 
q 
Hitunglah error yang terjadi (Et) pada penyelesaian ODE tersebut;
hitunglah komponen error setiap suku pada persamaan Et.
Selesaikan ODE tersebut dengan memakai h = 0.25;
bandingkan dengan penyelesaian sebelumnya;
bandingkan juga error yang terjadi.
Baca buku acuan pada hlm 580-584 untuk membantu Sdr dalam membuat
diskusi hasil hitungan Sdr.
Metode Euler
26
http://istiarto.staff.ugm.ac.id
q 
Error pada Metode Euler dapat dihitung dengan memanfaatkan Deret Taylor
q 
Keterbatasan
q 
q 
q 
Deret Taylor hanya memberikan perkiraan/estimasi local truncation error, yaitu error yang
timbul pada satu langkah hitungan Metode Euler, bukan propagated truncation error.
Hanya mudah dipakai apabila ODE berupa fungsi polinomial sederhana yang mudah
untuk di-diferensial-kan, fi(xi,yi) mudah dicari.
Perbaikan Metode Euler, memperkecil error
q 
Pakailah selang h kecil.
q 
Metode Euler tidak memiliki error apabila ODE berupa fungsi linear.
http://istiarto.staff.ugm.ac.id
27
Penyelesaian ODE
Metode Heun
http://istiarto.staff.ugm.ac.id
Metode Heun
28
http://istiarto.staff.ugm.ac.id
Slope di selang antara xi dan di xi+1 ditetapkan
sebagai nilai rata-rata slope di awal dan di akhir
selang, yaitu di xi dan di xi+1:
y
yiʹ′ = f (xi , yi )
φi + φ0i +1
slope =
2
yi0+1 = yi + f (xi , yi )h
Euler à sebagai prediktor
yiʹ′+1 = f (xi +1 , yi0+1 )
slope = φi
slope = φ0i+1 xi
xi+1
step size = h
x
yiʹ′ + yiʹ′+1 f (xi , yi ) + f (xi +1 , yi0+1 )
yʹ′ =
=
2
2
f (xi , yi ) + f (xi +1 , yi0+1 )
yi +1 = yi +
h à sebagai korektor
2
Metode Heun
29
http://istiarto.staff.ugm.ac.id
q 
Pakailah Metode Heun untuk mengintegralkan ODE di bawah ini, dari x = 0
s.d. x = 4 dengan selang langkah h = 1
yʹ′ =
q 
q 
dy
= 4e 0.8 x − 0.5y
dx
Syarat awal yang diterapkan pada ODE tsb adalah bahwa pada x = 0,
y=2
Penyelesaian eksak ODE tsb yang diperoleh dari kalkulus adalah:
y=
4 0.8 x −0.5 x
(e − e )+ 2e−0.5x
1.3
Metode Heun
30
http://istiarto.staff.ugm.ac.id
q 
Selang ke-1, dari x0 = 0 s.d. x1 = x0 + h = 1:
[
]
f (x0 , y0 ) = f (0,2) = 4 e 0.8 (0 ) − 0.5(2) = 3
slope di titik ujung awal, (x0, y0)
y10 = y0 + f (x0 , y0 )h = 2 + 3 (1) = 5
prediktor y1
[
]
f (x1 , y10 ) = f (1,5) = 4 e0.8(1) − 0.5(5) = 6.4021637
yʹ′ =
3 + 6.4021637
= 4.70108185
2
y1 = y0 + y ʹ′ h = 2 + 4.70108185(1) = 6.7010819
slope di titik ujung akhir, (x1, y1)
slope rata-rata selang ke-1
korektor y1
Metode Heun
31
http://istiarto.staff.ugm.ac.id
yi (eksak)
f(xi,yi)awal
yi (prediktor)
f(xi,yi)akhir
f(xi,yi)rerata
---
---
---
yi (korektor)
εt
2
---
i
xi
0
0
2
3
1
1
6.1946
5.5516
5.0000
6.4022
4.7011
6.7011
-8%
2
2
14.8439
11.6522
12.2527
13.6858
9.6187
16.3198
-10%
3
3
33.6772
25.4931
27.9720
30.1067
20.8795
37.1992
-10%
4
4
75.3390
---
62.6923
66.7840
46.1385
83.3378
-11%
Metode Heun
32
http://istiarto.staff.ugm.ac.id
90
Heun
60
Eksak
Y
30
0
0
1
2
X
3
4
Metode Heun
33
http://istiarto.staff.ugm.ac.id
q 
Metode Heun dapat diterapkan secara iteratif pada saat
menghitung slope di ujung akhir selang dan nilai yi+1
korektor
q 
q 
nilai yi+1 korektor pertama dihitung berdasarkan nilai yi+1
prediktor
q 
nilai yi+1 korektor tersebut dipakai sebaga nilai yi+1 prediktor
q 
hitung kembali nilai yi+1 korektor yang baru
q 
ulangi kedua langkah terakhir tersebut beberapa kali
Perlu dicatat bahwa
q 
error belum tentu selalu berkurang pada setiap langkah iterasi
q 
iterasi tidak selalu konvergen
f (xi , yi ) + f (xi +1 , yi0+1 )
yi +1 = yi +
h
2
f (xi , yi ) + f (xi +1 , yi0+1 )
yi +1 = yi +
h
2
Metode Heun
34
http://istiarto.staff.ugm.ac.id
q 
Iterasi kedua pada selang ke-1, dari x0 = 0 s.d. x1 = x0 + h = 1:
y10 = y1(old ) = 6.7010819
prediktor y1 = korektor y1(lama)
f (x1 , y10 ) = f (1,6.7010819)
slope di titik ujung akhir, (x1, y1)
[
]
= 4 e 0.8 (1 ) − 0.5(6.7010819)
= 5.5516228
yʹ′ =
3 + 5.5516228
= 4.2758114
2
y1 = y0 + y ʹ′ h = 2 + 4.27581145(1) = 6.2758114
q 
Iterasi di atas dapat dilakukan beberapa kali
slope rata-rata selang ke-1
korektor y1
http://istiarto.staff.ugm.ac.id
35
Penyelesaian ODE
Metode Poligon
(Modified Euler Method)
http://istiarto.staff.ugm.ac.id
Metode Poligon
36
http://istiarto.staff.ugm.ac.id
Slope di selang antara xi dan di xi+1 ditetapkan
sebagai nilai slope di titik tengah selang, yaitu di
xi+½:
y
yiʹ′ = f (xi , yi )
slope = φi+½
slope di titik awal
yi + 1 = yi + f (xi , yi )
2
slope = φi
slope = φi+½
(
yiʹ′+ 1 = f xi+ 1 , yi+ 1
2
xi
xi+½ xi+1
step size = h
x
2
2
(
h
2
ekstrapolasi ke titik tengah
)
slope di titik tengah
)
yi+1 = yi + f xi+ 1 , yi+ 1 h
2
2
ekstrapolasi ke titik akhir
Metode Poligon
37
http://istiarto.staff.ugm.ac.id
q 
Pakailah Metode Poligon untuk mengintegralkan ODE di bawah ini, dari x = 0
s.d. x = 4 dengan selang langkah h = 1
yʹ′ =
q 
q 
dy
= 4e 0.8 x − 0.5y
dx
Syarat awal yang diterapkan pada ODE tsb adalah bahwa pada x = 0,
y=2
Ingat penyelesaian eksak ODE tsb yang diperoleh dari kalkulus adalah:
y=
4 0.8 x −0.5 x
(e − e )+ 2e−0.5x
1.3
Metode Poligon
38
http://istiarto.staff.ugm.ac.id
q 
Selang ke-1, dari x0 = 0 sd x1 = x0 + h = 1:
[
]
f (x0 , y0 ) = f (0,2) = 4 e 0.8 (0 ) − 0.5(2) = 3
slope di titik ujung awal, (x0, y0)
h
y 1 = y0 + f (x0 , y0 ) = 2 + 3 (1 2) = 3.5
2
2
titik tengah y½
( )
[
]
f x 1 , y 1 = f (0.5,3.5) = 4 e0.8(0.5) − 0.5(3.5) = 4.2173 slope di titik tengah, (x½, y½)
2
2
( )
y1 = y0 + f x 1 , y 1 h = 2 + 4.2173(1) = 6.2173
2
2
ekstrapolasikan y0 ke y1
Metode Poligon
39
http://istiarto.staff.ugm.ac.id
yi (eksak)
f(xi,yi)
xi
0
0
2
3
1
1
6.1946
5.7935
0.5
3.5000
4.2173
6.2173
-0.4%
2
2
14.8439
12.3418
1.5
9.1141
8.7234 14.9407
-0.7%
3
3
33.6772
27.1221
2.5 21.1116
19.0004 33.9412
-0.8%
4
4
75.3390
3.5 47.5022
42.0275 75.9686
-0.8%
---
xi +½
yi +½
f(xi +½,yi +½)
---
---
---
εt
i
yi
2
---
Metode Poligon
40
http://istiarto.staff.ugm.ac.id
90
Poligon
60
Eksak
Y
30
0
0
1
2
X
3
4
http://istiarto.staff.ugm.ac.id
41
Penyelesaian ODE
Metode Runge-Kutta
http://istiarto.staff.ugm.ac.id
Metode Runge-Kutta
42
http://istiarto.staff.ugm.ac.id
q 
Metode Euler
q 
q 
q 
kurang teliti
ketelitian lebih baik diperoleh dengan cara memakai pias kecil atau
memakai suku-suku derivatif berorde lebih tinggi pada Deret Taylor
Metode Runge-Kutta
q 
q 
lebih teliti daripada Metode Euler
tanpa memerlukan suku derivatif
Metode Runge-Kutta
43
http://istiarto.staff.ugm.ac.id
q 
Bentuk umum penyelesaian ODE dengan Metode Runge-Kutta adalah:
yi+1 = yi + φ(xi , yi , h) h
q 
φ(xi , yi , h) adalah increment function yang dapat
diinterpretasikan sebagai slope atau gradien
fungsi y pada selang antara xi s.d. xi+1
Fungsi φ dapat dituliskan dalam bentuk umum sbb:
φ = a1k1 + a2k2 + ... + ankn
a adalah konstanta dan k adalah: k1 = f (xi , yi )
setiap k saling terhubung dengan
k yang lain à k1 muncul pada
pers k2 dan k2 muncul pada pers
k3 dst.
k2 = f (xi + p1h, yi + q11k1h)
k3 = f (xi + p2h, yi + q21k1h + q22k2h)
kn = f (xi + pn−1h, yi + qn−1,1k1h + qn−1,2k2h + ... + qn−1,n−1kn−1h)
Metode Runge-Kutta
44
http://istiarto.staff.ugm.ac.id
q 
q 
Terdapat beberapa jenis Metode Runge-Kutta yang dibedakan dari jumlah
suku pada persamaan untuk menghitung k:
q  RK tingkat-1 (first-order RK):
n=1
yi +1 = yi + φ(xi , yi , h) h
q  RK tingkat-2 (second-order RK): n = 2
φ(xi , yi , h) = a1k1 + a2k2 + ... + ankn
q  RK tingkat-3 (third-order RK):
n=3
q  RK tingkat-4 (fourth-order RK): n = 4
Order of magnitude kesalahan penyelesaian Metode RK tingkat n:
q  local truncation error = O(hn+1)
q  global truncation error = O(hn)
Second-order Runge-Kutta Method
45
http://istiarto.staff.ugm.ac.id
q 
Bentuk umum persamaan penyelesaian ODE dengan 2nd-order RK
yi+1 = yi + (a1k1 + a2k2 ) h
k1 = f (xi , yi )
k2 = f (xi + p1h, yi + q11k1h )
a1, a2, p1, q11 unknowns à perlu 4 persamaan
q 
Deret Taylor
h2
yi +1 = yi + f (xi , yi )h + f ʹ′(xi , yi )
2
⎛ ∂f ∂x d y ⎞ h2
yi +1 = yi + f (xi , yi )h + ⎜ +
⎟
⎝ ∂x ∂y d x ⎠ 2
f ʹ′(xi , yi ) =
∂f ∂f d y
+
∂x ∂y d x
Second-order Runge-Kutta Method
46
http://istiarto.staff.ugm.ac.id
q 
Ingat, Deret Taylor untuk fungsi yang memiliki 2 variabel
g(x + r , y + s ) = g(x , y ) + r
q 
∂g
∂g
+ s + ...
∂x
∂y
Bentuk di atas diterapkan pada persamaan k2
∂f
∂f
+ q11k1h + O(h2 )
∂x
∂y
Bentuk umum penyelesaian ODE metode 2nd-order RK menjadi:
∂f
∂f
yi+1 = yi + a1hf xi ,yi + a2hf xi ,yi + a2 p1h2
+ a2q11h2 f xi ,yi
+ O h3
∂x
∂x
⎡
∂f
∂f ⎤ 2
3
= yi + ⎡⎣a1 f xi ,yi + a2 f xi ,yi ⎤⎦ h + ⎢a2 p1
+ a2q11 f xi ,yi
⎥ h +O h
∂x
∂x ⎦
⎣
!
k2 = f (xi + p1h, yi + q11k1h) = f (xi , yi ) + p1h
q 
(
)
(
)
(
)
(
)
(
(
)
)
( )
( )
Second-order Runge-Kutta Method
47
http://istiarto.staff.ugm.ac.id
q 
Bandingkan persamaan di atas dengan persamaan semula
∂f
∂f ⎤
⎡
yi +1 = yi + [a1 f (xi , yi ) + a2 f (xi , yi )] h + ⎢a2 p1 + a2q11 f (xi , yi ) ⎥ h2 + O(h3 )
∂x
∂x ⎦
⎣
⎛ ∂f ∂x d y ⎞ h2
yi +1 = yi + f (xi , yi )h + ⎜ +
⎟
∂
x
∂
y
d
x
⎝
⎠ 2
q 
Agar kedua persamaan di atas ekuivalen, maka:
a1 + a2 = 1
a2 p1 =
1
q 
2
a2q11 = 1 2
q 
Karena hanya ada 3 persamaan untuk 4 unknowns, maka nilai
salah satu variabel harus ditetapkan.
Misalkan nilai a2 ditetapkan, maka a1, p1, dan q11 dapat dihitung.
Second-order Runge-Kutta Method
48
http://istiarto.staff.ugm.ac.id
q 
Jika a2 ditetapkan, maka:
a1 + a2 = 1
a2 p1 = 1 2
a2q11 = 1 2
a1 = 1 − a2
p1 = q11 =
q 
1
2a2
q 
Karena ada sejumlah tak berhingga nilai a2,
maka terdapat pula sejumlah tak berhingga
2nd-order RK methods.
Setiap versi 2nd-order RK akan memberikan
hasil yang persis sama jika fungsi
penyelesaian ODE yang dicari adalah fungsi
kuadrat, linear, atau konstanta.
Second-order Runge-Kutta Method
49
http://istiarto.staff.ugm.ac.id
q 
Metode Heun dengan korektor tunggal
a2 = 12 ⇒ a1 = 12 , p1 = q11 = 1
yi+1 = yi + (12 k1 + 12 k2 ) h k1 = f (xi , yi )
k2 = f (xi + h, yi + hk1 )
q 
Metode poligon yang diperbaiki (improved polygon method)
a2 = 1 ⇒ a1 = 0, p1 = q11 = 12
yi+1 = yi + k2 h
k1 = f (xi , yi )
k2 = f (xi + 12 h, yi + 12 hk1 )
q 
Metode Ralston
a2 = 2 3 ⇒ a1 = 1 3 , p1 = q11 = 3 4
yi+1 = yi + (13 k1 + 23 k2 ) h
k1 = f (xi , yi )
k2 = f (xi + 43 h, yi + 43 hk1 )
Second-order Runge-Kutta Method
50
http://istiarto.staff.ugm.ac.id
q 
Pakailah berbagai 2nd-order RK methods untuk mengintegralkan ODE di
bawah ini, dari x = 0 s.d. x = 4 dengan selang langkah h = 0.5:
f (x , y ) =
q 
q 
dy
= −2 x 3 + 12 x 2 − 20 x + 8.5
dx
Syarat awal yang diterapkan pada ODE tsb adalah bahwa di titik x= 0,
y=1
Bandingkan hasil-hasil penyelesaian dengan berbagai metode RK tsb.
Second-order Runge-Kutta Method
51
http://istiarto.staff.ugm.ac.id
i
xi
0
0
1
0.5
2
1
3
1.5
4
2
5
2.5
6
3
7
3.5
8
4
Single-corrector Heun
k1
k2
yi (eksak)
φ
εt
yi
1
8.5
1.25
4.875
1
0.0%
3.21875
1.25
-1.5
-0.125
3.4375
-6.8%
3
-1.5
-1.25
-1.375
3.375
-12.5%
2.21875
-1.25
0.5
-0.375
2.6875
-21.1%
2
0.5
2.25
1.375
2.5
-25.0%
2.71875
2.25
2.5
2.375
3.1875
-17.2%
4
2.5
-0.25
1.125
4.375
-9.4%
4.71875
-0.25
-7.5
-3.875
4.9375
-4.6%
3
---
---
---
3
0.0%
Second-order Runge-Kutta Method
52
http://istiarto.staff.ugm.ac.id
i
xi
0
0
1
0.5
2
1
3
1.5
4
2
5
2.5
6
3
7
3.5
8
4
Improved Polygon
k1
k2
yi (eksak)
φ
εt
yi
1
8.5
4.21875
4.21875
1
0.0%
3.21875
1.25
-0.59375
-0.59375
3.109375
3.4%
3
-1.5
-1.65625
-1.65625
2.8125
6.3%
2.21875
-1.25
-0.46875
-0.46875
1.984375
10.6%
2
0.5
1.46875
1.46875
1.75
12.5%
2.71875
2.25
2.65625
2.65625
2.484375
8.6%
4
2.5
1.59375
1.59375
3.8125
4.7%
4.71875
-0.25
-3.21875
-3.21875
4.609375
2.3%
3
---
---
---
3
0.0%
Second-order Runge-Kutta Method
53
http://istiarto.staff.ugm.ac.id
i
xi
0
0
1
0.5
2
1
3
1.5
4
2
5
2.5
6
3
7
3.5
8
4
Second-order Ralston Runge-Kutta
yi (eksak)
k1
k2
φ
εt
yi
1
8.5
2.582031
4.554688
1
0.0%
3.21875
1.25
-1.15234
-0.35156
3.277344
-1.8%
3
-1.5
-1.51172
-1.50781
3.101563
-3.4%
2.21875
-1.25
0.003906
-0.41406
2.347656
-5.8%
2
0.5
1.894531
1.429688
2.140625
-7.0%
2.71875
2.25
2.660156
2.523438
2.855469
-5.0%
4
2.5
0.800781
1.367188
4.117188
-2.9%
4.71875
-0.25
-5.18359
-3.53906
4.800781
-1.7%
3
---
---
---
3.03125
-1.0%
Second-order Runge-Kutta Method
54
http://istiarto.staff.ugm.ac.id
6
Y
4
Exact
Heun
Improved Polygon
2
Ralston
X
0
0
1
2
3
4
Third-order Runge-Kutta Method
55
http://istiarto.staff.ugm.ac.id
q 
Persamaan penyelesaian ODE 3rd-order RK methods:
yi+1 = yi + [16 (k1 + 4k2 + k3 )] h
k1 = f (xi , yi )
k2 = f (xi + 12 h, yi + 12 hk1 )
k3 = f (xi + h, yi − hk1 + 2hk2 )
q 
Catatan:
q 
Jika derivatif berupa fungsi x saja, maka 3rd-order RK sama dengan
persamaan Metode Simpson ⅓
Third-order Runge-Kutta Method
56
http://istiarto.staff.ugm.ac.id
i
xi
0
0
1
0.5
2
1
3
1.5
4
2
5
2.5
6
3
7
3.5
8
4
yi (eksak)
Third-order Runge-Kutta
k1
k2
k3
φ
εt
yi
1
8.5
4.219
1.25
4.438
1
0.0%
3.21875
1.25
-0.594
-1.5
-0.438
3.21875
0.0%
3
-1.5
-1.656
-1.25
-1.563
3
0.0%
2.21875
-1.25
-0.469
0.5
-0.438
2.21875
0.0%
2
0.5
1.469
2.25
1.438
2
0.0%
2.71875
2.25
2.656
2.5
2.563
2.71875
0.0%
4
2.5
1.594
-0.25
1.438
4
0.0%
4.71875
-0.25
-3.219
-7.5
-3.438
4.71875
0.0%
3
---
---
---
---
3
0.0%
Fourth-order Runge-Kutta Method
57
http://istiarto.staff.ugm.ac.id
q 
Persamaan penyelesaian ODE 4th-order RK methods:
yi+1 = yi + [16 (k1 + 2k2 + 2k3 + k4 )] h
k1 = f (xi , yi )
k2 = f (xi + 12 h, yi + 12 hk1 )
k3 = f (xi + 12 h, yi + 12 hk2 )
k4 = f (xi + h, yi + hk3 )
q 
Catatan:
q 
Jika derivatif berupa fungsi x saja, maka 4th-order RK sama dengan
persamaan Metode Simpson ⅓
Fourth-order Runge-Kutta Method
58
http://istiarto.staff.ugm.ac.id
k1
Fourth-order Runge-Kutta
k2
k3
k4
φ
εt
i
xi
yi (eksak)
yi
0
0
1
8.5
4.219
4.219
1.25
4.44
1
0.0%
1
0.5
3.21875
1.25
-0.594
-0.594
-1.5
-0.44
3.21875
0.0%
2
1
3
-1.5
-1.656
-1.656
-1.25
-1.56
3
0.0%
3
1.5
2.21875
-1.25
-0.469
-0.469
0.5
-0.44
2.21875
0.0%
4
2
2
0.5
1.469
1.469
2.25
1.44
2
0.0%
5
2.5
2.71875
2.25
2.656
2.656
2.5
2.56
2.71875
0.0%
6
3
4
2.5
1.594
1.594
-0.25
1.44
4
0.0%
7
3.5
4.71875
-0.25
-3.219
-3.219
-7.5
-3.44
4.71875
0.0%
8
4
3
---
---
---
---
---
3
0.0%
3rd- and 4th-order Runge-Kutta Methods
59
http://istiarto.staff.ugm.ac.id
q 
Pakailah 3rd-order dan 4th-order RK methods untuk mengintegralkan ODE di
bawah ini, dari x = 0 s.d. x = 4 dengan selang langkah h = 0.5:
f (x , y ) =
dy
= −2 x 3 + 12 x 2 − 20 x + 8.5
dx
§  Syarat awal yang diterapkan pada ODE tsb adalah bahwa di titik x = 0,
y=1
3rd- and 4th-order Runge-Kutta Methods
60
http://istiarto.staff.ugm.ac.id
Node
Exact Solution
i
xi
0
0
1
0.5
2
1
3
1.5
4
2
5
2.5
6
3
7
3.5
8
4
Third-order RK
yi
Fourth-order RK
εt
yi
εt
yi
1
1
0.0%
1
0.0%
3.21875
3.21875
0.0%
3.21875
0.0%
3
3
0.0%
3
0.0%
2.21875
2.21875
0.0%
2.21875
0.0%
2
2
0.0%
2
0.0%
2.71875
2.71875
0.0%
2.71875
0.0%
4
4
0.0%
4
0.0%
4.71875
4.71875
0.0%
4.71875
0.0%
3
3
0.0%
3
0.0%
61
http://istiarto.staff.ugm.ac.id
Fly UP