5
MASALAH NILAI AWAL PERSAMAAN DIFERENSIAL BIASA
Pada bab ini dibahas tentang persamaan diferensial biasa, ordinary differential equations (ODE) yang diklasifikasikan kedalam masalah nilai awal (initial value) dan masalah nilai batas (boundary value), dimana kedua keadaan ini solusinya dispesifikasi pada waktu awal (initial time). Akan disajikan beberapa metode pendekatan komputasi numerik untuk menangani permasalahan yang berkaitan dengan persamaan diferensial biasa, dengan beberapa contoh kasus terapan.
A. SASARAN UMUM Sasaran umum dari perkuliahan ini adalah memberikan pemahaman kepada mahasiswa mengenai penyelesaian persamaan diferensial biasa ke dalam model komputasi numerik, yang banyak dijumpai dalam formulasi fenomena atau sistem fisika.
B. SASARAN KHUSUS Setelah perkuliahan selesai dilaksanakan, mahasiswa diharapkan mampu: 1. Menjelaskan konsep nilai awal dan nilai batas dalam persamaan diferensial biasa dalam pendekatan komputasi numerik. 2. Menyebutkan beberapa metode yang populer digunakan dalam menyelesaikan permasalahan persamaan diferensial biasa. 3. Menjelaskan strategi dasar dari beberapa metode komputasi numerik untuk persamaan diferensial biasa 4. Mengimplementasikan dalam bentuk program tentang persamaan diferensial biasa dalam beberapa kasus yang ditangani.
C. URAIAN MATERI Äfisika-komputasi ⊇
86
Banyak hukum-hukum fisika yang
‘sangat pas ’ diformulasikan dalam
bentuk persamaan diferensial. Lebih lanjut, tidak mengherankan bahwa solusi komputasi numerik dari persamaan-persamaan diferensial menjadi bagian yang umum dalam pemodelan sistem-sistem fisika.
Beberapa hukum mendasar
diantaranya sebagai berikut:
Hukum Hukum Newton II tentang gerak Hukum Panas Fourier
Formulasi Matematika dv F = dt m q = −k '
Hukum Faraday
dT dx
∆V L = L
di dt
Variabel & Parameter Kecepatan (v), gaya (F) dan massa (m) Flux panas(q), konduktivitas termal (k’) dan temperatur (T) Tegangan Drop ( ∆V L ), induktansi (L) dan arus (i)
Bentuk umum persamaan diferensial biasa adalah set M pasang persamaan orde satu, dy = f ( x, y) dx
(5.1)
dimana x adalah variabel bebas, dan y adalah sebuah set dari M variabel takbebas ( f adalah vektor komponen
M). Persamaan diferensial orde yang lebih tinggi bisa
dituliskan dalam bentuk persamaan diferensial orde satu. Gerak 1 dimensi sebuah partikel bermassa m, dibawah pengaruh gaya sebesar F(z), dituliskan dalam bentuk persamaan orde dua: m
d2z = F ( z) dt
(5.2)
Jika momentum didefinisikan sebagai perkalian massa dengan kecepatan, yang dituliskan sebagai p (t ) = m
dz dt
(5.3)
maka persamaan (5.2) menjadi dua pasang persamaan orde satu (Hamilton), dz p dp = dan = F (z ) dt m dt
(5.4)
Äfisika-komputasi ⊇
87
Akan dibahas beberapa metode untuk menyelesaikan persamaan diferensial biasa, dengan penekanan pada masalah nilai awal. Artinya, mencari y(x) yang diberikan oleh nilai y pada beberapa titik awal, y(x=0)=y0. Kasus atau masalah dalam tipe ini, diantaranya saat memberikan posisi awal dan momentum suatu partikel dan diharapkan memberikan atau menemukan gerak selanjutnya, menggunakan persamaan (5.3). Masalah nilai awal pada persamaan diferensial biasa bisa dituliskan dalam bentuk y ' (t ) = f ( y , t ) ; y (0 ) = y 0
(5.5)
dimana f(y,t) adalah fungsi terhadap y dan t, dan persamaan kedua adalah keadaan awal. Pada persamaan (5.5) turunan pertama terhadap y diberikan sebagai fungsi y dan t, dan akan dicari fungsi y yang tidak diketahui dengan melakukan integrasi f’(y,t). Banyak contoh untuk masalah nilai awal persamaan diferensial biasa, antara lain: (a) y' (t ) = 3 y + 5, y( 0) = 1 (b) y' (t ) = ty + 1, y(0 ) = 0 (c) y' (t ) = −
1 , y (0) = 1 1+ y2
(d) y' = z, z' = − y, y(0 ) = 1, z(0 ) = 0 Prinsip metode komputasi numerik untuk persamaan diferensial bia sa adalah menentukan solusi pada titik-titik, tn =t n-1+h, dimana h adalah lebar langkah ( atau interval waktu) . Ada beberapa pemilahan tipe metode, diantaranya metode Euler, metode Runge-Kutta, dan metode Predictor-Corrector. Berdasarkan kaku tidaknya persa maan (stiff equation), ada kategori nonstiff equation, meliputi metode Euler, metode Runge -Kutta dan metode Predictor-Corrector. Sedangkan kategori stiff equation meliputi metode implisit dan metode transformasi eksponensial. Dengan pertimbangan ragam pema haman, pada bab ini ditelaah 3 macam metode berdasarkan tingkatan kasus dan metode yang ditangani yaitu metode Sederhana, metode Implisit dan metode Runge-Kutta.
Äfisika-komputasi ⊇
88
5.1 METODE SEDERHANA Secara umum, solusi persamaan diferensial dari persamaan (5.4) didekati dengan keadaan awal, y(x=0)=y 0. Lebih spesifik, sebenarnya bisa dianalisis pada nilai y dengan sebarang x, misalnya 1. Strategi umum pada metode sederhana ini adalah membagi interval [0,1] menjadi sejumlah N,
sehingga lebar sub-interval sama h=1/N, kemudian
membangun sebuah formula rekursi terkait yn pada { yn-1 , yn-2 ,…} dimana yn adalah pendekatan pada y(xn =nh). Sebagaimana sebuah rekursi, pola integrasi selangkah demi langkah pada persamaan diferensial dari x=0 sampai x=1. Salah satu algoritma metode sederhana adalah metode Euler, dimana persamaan (5.1) ditinjau pada titik xn dengan pendekatan diferensiasi beda maju, seperti yang telah dibahas pada bab IV, sehingga y n +1 − y n + O( h ) = f ( x n , y n ) h
(5.6)
sehingga relasi rekursi yn+1 yang dinyatakan da lam bentuk yn adalah yn + 1 = yn + hf ( xn , yn ) + O( h 2 )
(5.7)
Formula ini memiliki kesalahan lokal ( yang terjadi karena langkah tunggal dari yn ke yn+1) sebesar O(h2). Kesalahan global yang terjadi dalam mencari y(1) dengan memberikan sejumlah N langkah integrasi dari x=0 sampai x=1 sebesar N kali O(h2)
O(h).
Contoh 5.1 Selesaikan set persamaan diferensial biasa berikut dengan metode dengan h=0,005ð dan h=0,0005 ð : y' = z , z ' = − y,
y(0) = 1 z( 0) = 0
solusi Kalkulasi untuk langkah pertama de ngan h=0,0005ð ditunjukkan berikut ini.
Äfisika-komputasi ⊇
89
t 0 = 0;
y0 = 1 z0 = 0
t1 = 0 ,0005 π :
y1 = y0 + hz0 = 1 + ( 0,0005 π )(0 ) = 1 .0 z1 = z 0 − hy0 = 0 − ( 0,0005 π )(1) = −0 ,00157
t2 = 0,001π :
y2 = y1 + hz1 = 1 + (0,0005 π )(−0 ,00157 ) = 0,99999 z 2 = z1 − hy1 = −0,00157 − (0,0005 π )(1) = −0,00314
pada tabel dibawah disajikan hasil perhitungan untuk nilai-nilai t yang dipilih dibandingkan dengan solusi eksak, y=cos(t) dan z=-sin(t). Eksak
t
h=0,005ð
h=0,0005 ð
y
z
y
z
y
Z
0.5 ð
0
-1
1,3E-4
-1,01241
2,62E-6
-1,00123
ð
-1
0
-1,02497
-2,67E-4
-1,00247
-5,25E-6
1,5 ð
0
1
-4,01E-4
1,03770
-7,88E-6
1,00371
2ð
1
0
1,05058
5,48E-4
1,00495
1,05E-5
3ð
-1
0
-1,07682
-8,43E-4
-1,00743
-1,58E-5
6ð
1
0
1,15954
1,82E-3
1,01491
3,19E-5
8ð
1
0
1,21819
2,54E-3
1,01994
4,27E-5
Dari hasil diatas menunjukkan bahwa kesalahan pada y bertambah dengan pertambahan t, dan sebanding dengan h ( lihat bahwa nilai y untuk t= ð, 2ð, 3ð, 6ð dan 8ð: nilai z untuk t tersebut tidak mengikuti trend yang sama, karena ketika z nol , kesalahan z secara signifikan dipengaruhi oleh pergeseran fase.
Contoh 5.2 Buatlah program untuk mencari nilai pendekatan komputasi numerik, dari persamaan diferensial dengan keadaan batas, dy = − xy; y(0 ) = 1 dx Integrasikan maju dari x=0 sampai x=3, dimana solusinya adalah e − 12 x 2 . solusi Program dalam BASIC Inputnya adalah lebar langkah (h) dan output yang ditampilkan adalah hasil integrasi dan kesalahan yang terjadi. Äfisika-komputasi ⊇
90
10 DEF FNF(X,Y)=-X*Y 20 INPUT “masukkan nilai h”;H 30 IF H <= 0 THEN STOP 40 N=3/H 50 Y=1 60 FOR I%=0 TO N%-1 70
X=I%*H
80
Y=Y+H*FNF(X,Y)
90
DIFF=EXP(-0,5*(X+H)^2)-Y
100
PRINT I%, X+H, Y, DIFF
110 NEXT I% 120 GOTO 20 Kesalahan jawaban atau hasil ditentukan dari y(1)=e-1/2 =0,606531 dan y(3)=e9/2
=0,011109. Tingkat kesalahan dalam integrasi dy/dx=-xy dengan y(0)=1, pada
variasi lebar langkah seperti terlihat pada tabel berikut. Euler dengan persamaan (5.7) H
y(1)
y(3)
0,500
-1,43469
0,011109
0,200
-0,46330
0,006519
0,100
-0,21625
0,03318
…
…
…
…
…
…
0,005
-0,01014
0,000167
0,002
-0,00405
0,000067
0,001
-0,00203
0,000038
Kesalahan akan berkurang secara linear dengan semakin kecilnya h. Bagaimanapun, kesalahan fraksi ( kesalahan dibagi oleh y) akan meningkat terhadap x ketika langkah semakin banyak diberikan di dalam integrasi dan y menjadi kecil. Meskipun metode Euler kelihatannya cukup baik, tetapi secara umum tidak memuaskan karena akurasinya rendah. Metode yang lebih tinggi satu klas dari Äfisika-komputasi ⊇
91
kesederhanaannya bisa diturunkan dari deret Taylor yang ekspansinya untuk yn+1 disekitar yn : yn +1 = y( xn + h ) = yn + hyn' +
1 2 '' h y n + O (h 2 ) 2
(5.8)
dari persamaan (5.1) kita memiliki yn ' = f ( xn , yn )
(5.9)
Jika disubtitusikan pada persamaan (5.8) hasilnya adalah: yn " =
df ∂f ∂f ∂y ∂ f ∂f ( xn , yn ) = + = + f dx ∂ x ∂y ∂x ∂x ∂y
(5.10)
subtitusikan ke (5.8) hasilnya adalah yn +1 = yn + hf +
1 2 ∂f ∂f h + f . + O(h 3 ) 2 ∂x ∂y
(5.11)
dimana f dan turunannya dievaluasi pada (xn ,yn). Relasi rekursi memiliki kesalahan lokal O(h3 ) dan kesalahan global O(h), satu tingkat lebih akurat daripada kesalahan pada metode Euler. Metode ini sangat bermanfaat ketika f diketahui secara analitik dan cukup sederhana untuk turunannya.
5.2 MULTI STEP & METODE IMPLISIT Metode lain yang akurasinya lebih baik, adalah menggunakan relasi rekursi yang menghubungkan yn+1 tidak hanya pada yn , tetapi juga pada titik-titik lebih lanjut, katakan yn-1 ,yn-2 ,…. Untuk menurunkan formulanya, kita integralkan satu langkah persamaan diferensial (5.1) y n +1 = y n +
xn +1
∫ f ( x, y) dx
(5.12)
xn
Masalahnya adalah kita tidak tahu f melampaui interval integrasi. Tapi kita bisa gunakan nilai y pada xn dan xn-1 untuk menentukan ekstrapolasi linear f, melalui interval f ≈
( x − xn −1 ) ( x − xn) fn − f n −1 + O( h 2 ) h h
(5.13)
dimana f i ≡ f ( xi , yi ) . Subtitusikan ke dalam (5.12) dan kerjakan integral terhadap x, lalu hasiln ya dalam metode 2 langkah Adam-Bashforth adalah: Äfisika-komputasi ⊇
92
1 3 yn +1 = yn + h f n − f n −1 + O(h 3 ) 2 2
(5.14)
Dengan cara yang sama, metode -metode lain untuk orde tinggi bisa diturunkan dengan ekstrapolasi terhadap polinomial orde tinggi. Sebagai contoh, jika f diekstrapolasi oleh polinomial kubik, dipetakan pada fn , fn-1 , fn-2 dan fn-3, metode empat langkah Adam-Bashforth menghasilkan: y n +1 = y n +
h (55 f n − 59 f n−1 + 37 f n−2 − 9 f n−3 ) + O( h 4 ) 24
(5.15)
Untuk metode implisit, persamaan diselesaikan dengan menentukan yn+1 , pada persamaan (5.1) pada sebuah titik xn + ≡ ( n + 12 )h di posisi tengah antara 2 1 2
titik-titik kisi: dy = f ( xn + , yn + ) dx xn + 1 2
1 2
1 2
(5.16)
Jika kita gunakan pendekatan simetri yang berbeda untuk turunan ( analog dengan f '≈
f 1− f −1 , h → 12 h ) dan menggantikan fn+1/2 dengan rata-rata nilai titik-titik kisi 2h
yang berdekatan ( kesalahan pada penggantian ini O(h2)), persamaan (5.16) dituliskan y n +1 − y n 1 + O( h 2 ) = ( f n + f n +1 ) + O(h 2 ) h 2
(5.17)
Bentuk (5.17) baik dan bagus, tetapi sifat yn+1 pada kedua sisi mengindikasikan harus ada penyelesaian non trivial ( mirip dengan Newton Raphson pada integrasi) pada setiap langkah integrasinya, yang tentunya cukup menghabiskan waktu. Penyederhanaannya adalah jika f adalah fungsi linear f(x,y)=g(x,y), sehingga pada kasus pada contoh 5.1, bisa diselesaikan 1 + 1 g ( x n )h y n +1 = 1 2 yn 1 − 2 g ( xn +1 ) h
(5.18)
Contoh 5.3 Selesaikan kasus contoh 5.2 dengan metode ekspansi deret Taylor , persamaan (5.10) dan metode Multiple & Implisit, persamaan (5.18). solusi Äfisika-komputasi ⊇
93
Dengan merubah fungsi yang dideklarasikan dengan program BASIC, pada line 80 list program contoh 5. 2, dengan fungsi rekursi Taylor (5.10) dan Implisit (5.18), maka program akan memberikan hasil output dalam besar kesalahan. Berikut adalah hasil setelah di running, Taylor dengan persamaan
Implisit dengan persamaan
(5.10)
(5.18)
H
0,500 0,200 0,100 … … 0,005 0,002 0,001
y(1)
y(3)
y(1)
y(3)
0,032312 0,005126 0,001273 … … 0.000003 0,000001 0,000000
-0,006660 -0,000712 -0,000149 … … 0,000000 0,000000 0,000000
-0,015691 -0,002525 -0,000631 … … -0,000001 0,000000 0,000000
0,001785 0,000255 0,000063 … … 0,000000 0,000000 0,000000
Dilihat dari hasil, dan dibandingkan dengan metode Euler pada contoh 5. 2, maka terlihat pendekatan Taylor dan Implisit memberikan hasil yang lebih baik pada lebar langkah yang sama. Dan pada kasus ini, metode Implisit lebih cepat konvergen dibanding Taylor maupun Euler.
5.3 METODE RUNGE-KUTTA Kekurangan utama dari metode Euler adalah tingkat akurasinya yang rendah. Keadaan ini membuat kerugian ganda. Untuk mencapai tingginya akurasi memerlukan h yang sangat kec il, disamping meningkatnya waktu komputasi dan mengakibatkan terjadinya kesalahan pembulatan (round off error). Pada metode Runge -Kutta, tingkat akurasinya meningkat oleh penggunaan titik-titik lanjutan pada tiap-tiap langkah interval. Akurasi yang lebih tinggi juga memberikan implikasi kesalahan berkurang lebih cepat dibanding metode akurasi tingkat rendah ketika h dikurangi. Tinjau persamaan diferensial biasa dengan nilai batas secara umum dituliskan y ' (t ) = f ( y , t ) ; y (0 ) = y 0
(5.19)
untuk menghitung yn +1 pada t n +1 = t n + h dengan nilai yn yang diketahui, kita integrasikan persamaan (5.19) pada interval [t n , t n +1 ] sebagai
Äfisika-komputasi ⊇
94
t n +1
yn +1 = yn +
∫ f ( y, t )dt
(5.20)
tn
dimana dari aturan trapesium pendekatan untuk integral: t n +1
∫
tn
1 f ( y, t )dt ≅ h [ f ( yn , tn ) + f ( yn +1 , tn +1 ) ] 2
(5.21)
pada persamaan (5.21) yn +1 tidak diketahui, sehingga bentuk kedua didekati oleh
f ( yn +1 , t n +1 ) , dimana y n +1 adalah perkiraan pertama untuk yn +1 yang dihitung dengan metode Euler. Skema yang didapatkan disini dinamakan metode Runge -Kutta dan diringkas sebagai berikut: yn +1 = yn + hf ( yn , tn )
[
h f ( yn ,t n ) + f ( yn +1 , tn +1 2 atau dalam bentuk yang lebih standar k1 = hf ( yn , t n ) yn +1 = yn −
]
k 2 = hf ( y n + k1 , t n +1 )
(5.22)
1 y n + 1 = y n + [k 1 + k 2 ] 2 Metode Runge -Kutta orde kedua ini identik dengan metode Euler Predictor Corrector, yaitu metode predictor corrector yang paling sederhana. Persamaan ini juga ekivalen dengan metode Euler modifikasi yang diaplikasikan hanya dengan satu langkah iterasi.
Contoh 5.4 Sebua h rangkaian ditunjukkan pada gambar 5.1 memiliki induktansi diri L=50H, resistensi R=20 ohm, dan sumber tegangan V=10 volt. Jika saklar ditutup pada t=0, arus I(t) memenuhi L
d I (t ) + RI (t ) = E , dt
I (0) = 0
Tentukan arus listrik yang mengalir untuk 0 t
(5.23) 10 detik menggunakan metode
Runge Kutta orde kedua dengan h=0,1 Solusi:
L
I(t)
E
Gambar 5.1 Rangkaian Listrik
Äfisika-komputasi ⊇ R
95
Persamaan (5.23) dituliskan kembali dalam bentuk d R E I (t ) + I (t ) + ≡ f ( I , t ) dt L L kemudian metode Runge Kutta orde dua menjadi: E R k 1 = h − I + L L E R k 2 = h − (In + k1) + L L 1 I n +1 = I n + (k1 + k 2) 2 Komputasi untuk dua langkah diperlihatkan sebagai berikut: n = 0(t = 0,1) : k1 = 0,1[( −0,4 )(0) + 0 ,2] = 0,02
k 2 = 0,1[( −0,4 )(0 + 0,02 ) + 0 ,02 ] = 0,0192 1 1 I1 = I 0 + ( k1 + k 2 ) = 0 = 0 + (0 ,02 + 0 ,0192 ) = 0,0196 2 2 n = 0(t = 0,1) : k1 = 0,1[( −0,4 )(0,0196 ) + 0,2 ] = 0,019216 k 2 = 0,1[( −0,4 )(0,0196 + 0,019216 ) + 0,02 ] = 0,018447 1 1 I 2 = I 1 + (k1 + k 2 ) = 0,0196 − ( 0,019216 + 0,018447 ) = 0,038431 2 2 Hasil akhir komputasi (sampai 10 langkah) sebagai berikut: t(detik) I(amp) t(detik) I(amp) 0 0 6 0,4546 1 0,1648 7 0,4695 2 0,2752 8 0,4796 3 0,3493 9 0,4863 4 0,3990 10 0,4908 5 0,4332 (∞ ) (0,500) Penggunaan metode Runge Kutta pada persamaan diferensial orde yang lebih tinggi, mudah. Sebagai ilustrasi, ditinjau persamaan diferensial orde kedua berikut: y" (t ) + ay' (t ) + by(t ) = q (t ), y(0 ) = 1, y' (0 ) = 0
(5.24)
dimana a dan b adalah koefisien dan q(t) adalah fungsi yang diketahui, dan dua keadaan awal diberikan. Dengan definisi z (t ) = y ' (t )
(5.25) Äfisika-komputasi ⊇
96
persamaan (5.24) bisa direduksi menjadi sepasang persamaan diferensial orde pertama: y ' = f ( y , z , t ) ≡ z,
y(0) = 1
z ' = g ( y, z, t ) ≡ − az − by + q ,
z( 0) = 0
(5.26)
sehingga metode Runge Kutta orde kedua untuk persamaan diatas dituliskan sebagai berikut: k1 = hf ( yn , z n , t n ) = hzn l1 = hg ( yn , z n , t n ) = h (− azn − bz = y n + q n ) k 2 = hf (( y n + k1 , zn + l1 , t n + 1) = h( zn + 1) l 2 = hg (( y n + k1 , z n + l1 , t n + 1) = h (− a ( z n +1 ) − b ( yn +1 ) + q n +1 ) 1 y n + 1 = y n + (k 1 + k 2 ) 2 1 z n +1 = z n + (l1 + l 2 ) 2
::: Studi Kasus Fisika ::: Sebuah kotak kubus bermassa M=0,5 kg terikat pada ujung bawah sebuah tali tak bermassa. Pada ujung tali atas diikatkan pada tiang diam. Kubus menerima resistansi R= -B dy/dtndari udara, dimana B adalah konstanta damping (lihat gambar 5.2). persamaan geraknya adalah: M
d2 d y + B y + ky = 0, 2 dt dt
y (0 ) = 1, y ' (0 ) = 0
(5.27)
dimana y adalah perpindahan dari posisi statis, k adalah konstanta pegas sebesar 100 kg/det2, dan B=10 kg/det. Hitunglah menggunakan metode Runge Kutta orde kedua: (a) y(t) untuk 0
Logam kubus
y
Gambar 5.2 Sistem massa pegas
Äfisika-komputasi ⊇
97
y' = z ≡ f ( y, z , t ),
y (0 ) = 1 B k z'= − z− y ≡ g ( y, z, t ), z (0) = 0 M M
dengan mengatur a =
B k = 20 , b = = 200 , , dan g=0, metode Runge Kutta orde M M
kedua untuk persamaan (5.27) menjadi bentuk persamaan (5.26) (a) Untuk n=1: t=0,025 k1 = hf ( y0 , z0 , t 0 ) = hz 0 = 0,025 (0 ) = 0 l1 = hg ( y 0 , z 0 , t 0 ) = h (−20 z 0 − 200 y0 ) = 0 ,025 (−20 (0 ) − 200 (1)) = −5 k 2 = hf (( y 0 + k1 , z 0 + l1 , t 0 ) = h ( z 0 + l1 ) = 0,025 (0 − 5 ) = −0,125 l 2 = hg (( y0 + k1 , z0 + l1 , t1 ) = h[( −20 ( z 0 + l1 ) − 200 ( y 0 ) + k1 )] = 0,025[ −20 (0 − 5) − 200 (1 + 0 )] = −2 ,5 1 y1 = y 0 + (0 − 0 ,125 ) = 0 .9375 2 1 z1 = z 0 + ( −5 − 2,5) = −3,75 2 untuk n=2: t=0,05 k 2 = hf ( y1 , z1 , t1 ) = hz1 = 0,025 (− 3,75 ) = −0 ,09375 l1 = hg ( y1 , z1 , t1 ) = h (−20 z1 − 200 y1 ) = 0 ,025 (−20 ( −3,75 ) − 200 (0 .9375 )) = − 2,8125 k 2 = hf (( y1 + k1 , z1 + l1 , t1 ) = h ( z1 + l1 ) = 0 ,025 (−3,75 − 2,8125 ) = 0,1640625 l 2 = hg (( y1 + k1 , z1 + l1 , t1 ) = h[(−20 ( z1 + l1 ) − 200 ( y1 + k1 )] = 0,025[ −20 (− 3,75 − 2,8125 ) − 200 (0 ,9375 − 0,093750 )] = −0 ,9375 1 y2 = y1 + ( −0,09375 − 0,1640625 ) = 0 .80859 2 1 z 2 = z1 + ( −2,8125 − 0,9375 ) = −5,625 2 (b) dan (c) pada bagian ini komputasi dikerjakan dengan menggunakan program, dan memberikan hasil setelah 50 langkah pada 0,75 detik, seperti terlihat pada tabel berikut: t (detik) 0 0,05
(b) (B=10) y (meter) 1,000 0,823
(c) (B=0) y (meter) 1,000 0,760 Äfisika-komputasi ⊇
98
0,1 0,15 0,2 0,25 0,3 0,35 0,4 0,45 0,5 0,55 0,6 0,65 0,7 0,75
0,508 0,238 0,066 -0,016 -0,042 -0,038 -0,025 -0,013 -0,004 0,000 0,001 0,001 0,001 0,000
0,155 -0,523 -0,951 -0,923 -0,45 0,235 0,810 0,996 0,705 0,075 -0,590 -0,973 -0,889 -0,378
Program dalam C # include <stdio.h> #include <stdlib.h> #include <math.h> /*
waktu: t y, z: y, y’ kount: jumlah langkah antara 2 jalur yang dicetak k, m, b: k, M(massa), B(koefisien damping)
*/
main() { int kount, n, kstep=0; float bm, k1, k2, km, l1, l2; static float time, k=100.0, m=0.5, b=10.0, z=0.0 static float y= 1.0, h=0,001; printf(“ Hasil Komputasi Metode Runge Kutta Orde Kedua \n”); printf(“
t
printf(“
%12.6f
y %12.5e
z\n”); %12.5e \n”, time, y, z);
km= k/m; bm=b/m; for ( n=1; n<=20; n++) { for ( kount =1; kount<=50; kount++) { kstep=kstep+1; time=h*kstep; Äfisika-komputasi ⊇
99
k1=h*z;
l1=-h*(bm*z+km*y);
k2=h*(z+11); l2=-h*(bm*(z+l1)+km*(y+k1)); y=y+ (k1 + k2)/2;
z=z+(11+12)/2;
} printf(“
%12.6f
%12.5e
%12.5e
\n”, time, y, z );
} exit (0) }
D. SOAL-SOAL (3.1) Selesaikan persamaan diferensial biasa berikut dalam 0
t 5 menggunakan
metode Euler, dengan h=0.5 y'+ty = 1, y'+3 y = e
−t
y' = (t 2 − y ), (3.2)
y (0 ) = 1 y (0 ) = 1 y(0 ) = 0 ,5
tanki kanonik berisi air dengan ketinggian 0,5m dari dasar. Tanki memiliki lobang dengan radius 0,02m pada dasarnya. Radius tanki dalam y diberikan oleh r=0,25y, dimana r adalah radius dan y adalah ketinggian yang diukur dari dasar tanki. Kecepatan air yang mengalir melalui lubang diberikan oleh v2=2gy dimana g=9,8 m/det2. Menggunakan metode Euler ( gunakan h=0,001 detik), tentukan berapa menit tanki menjadi kosong.
(3.3)
Hitunglah y(2) untuk persamaan berikut dengan metode Runge Kutta orde pertama, dengan h=1. y' = −
(3.4)
y , y(0) = 1 t + y2
Masalah nilai awal pada persamaan diferensial biasa diberikan oleh y ' ' ' = − y , y (0 ) = 1 y' (0 ) = y' ' ( 0) = 0 Gunakan metode Runge Kutta orde kedua dengan h=0,2, hitunglah y(0,4) dan y(1).
(3.5)
Buatlah program untuk menyelesaikan persamaan diferensial biasa pada soal (5.2) Äfisika-komputasi ⊇
100
E. DAFTAR PUSTAKA Chapra, S.C., and Canale, R.P., Numerical Methods for Engineers , McGraw-Hill, 1998 Gear, C. W., Numerical Initial Value Problems in Ordinary Differential Equations, Prentice -Hall, 1971 Hall, G and J.M Watt, Modern Numerical Methods for Ordinary Differential Equations , Clarendon Press, 1976 Koonin, S.E., Computational Physics, Addison-Wesley Inc, 1986 Mathews, J.H., Numerical Methods for Mathematics, Science and Engineering, Prentice -Hall Inc., 1992 Morris,J.L., Computational Methods in Elementary Numerical Analysis, Wiley, 1983 Nakamura, S., Applied Numerical Methods in C , Prentice-Hall Inc. 1993
Äfisika-komputasi ⊇
101
