BAB VI PERSAMAAN DIFERENSIAL PARSIAL (PDP)

Bab VI

Supardi, M.Si

BAB VI PERSAMAAN DIFERENSIAL PARSIAL (PDP) 1 Pendahuluan Persamaan diferensial parsial memegang peranan penting di dalam penggambaran keadaan fisis, dimana besaran-besaran yang terlibat didalamnya berubah terhadap ruang dan waktu. Sebagai contoh, jika kita meninjau topik-topik fisika lanjut (advanced physics), seperti halnya mekanika klasik lanjut yang membicarakan tentang gelombang elektromagnetik, hidrodinamik dan mekanika kuantum (gelombang Schroedinger), maka kita akan menemukan penggunaan persamaan diferensial parsial yang digunakan untuk menggambarkan fenomena fisis yang berkaitan dengan masalah-masalah tersebut. Masalah-masalah tersebut dalam kenyataannya sulit untuk dipecahkan dengan cara analitik biasa, sehingga metode numerik

perlu diterapkan

untuk menyelesaikannya.

Penggunaan persamaan

diferensial tidak terbatas pada masalah fisika saja, tetapi lebih luas lagi dalam bidang sains dan teknologi.

2 Pendekatan Beda Hingga Untuk memahami dengan benar masalah persamaan diferensial ini, sebelumnya pada bab 5 kita sudah membahas bahwa suatu derivatif dapat didekati dengan beda hingga, sehingga persamaan diferensial dapat didekati dengan persamaan beda hingga pula. Dalam bab ini metode beda hingga yang telah dikenalkan sebelumnya akan diperluas lagi untuk kasus di dalam ruang multidimensi yang lebih tepat dikaji dengan menggunakan persamaan diferensial parsial. Pada bab V yang lalu, kita sudah menggunakan pendekatan beda hingga untuk mendekati ungkapan turunan pertama dan kedua. Namun pada pembahasan yang lalu kita masih membatasi pada pendekatan untuk turunan pada ruang dimensi satu. Saat ini, kita masih akan menggunakan kaidah-kaidah pendekatan tersebut namun ditingkatkan untuk ruang dimensi dua. Pada pembahasan tentang persamaan diferensial biasa di bab 7 yang lalu, kita telah melakukan pendekatan beda hingga pada penyelesaiannya. Nah di bab ini, kita juga Persamaan Diferensial Parsial

Bab VI

Supardi, M.Si

akan melakukan hal sama pada bentuk derivatif parsialnya. Mengapa? Karena untuk masalah-masalah yang melibatkan dua atau lebih variabel bebas, prinsip-prinsip tersebut masih tetap dapat diterapkan. Di dalam pembahasan tentang persamaan diferensial biasa, variabel bebas yang terlibat dalam masalah hanya satu, sedangkan untuk persamaan diferensial parsial variabel bebas berjumlah lebih dari satu. Tentu saja, hal ini saja membuat permasalahan

akan

semakin

kompleks.

Untuk

memberikan

ilustrasi

dan

mempermudah pemahaman kita tentang masalah ini, sekarang marilah kita tinjau sebuah jaring kotak yang menggambarkan dua variabel bebas x dan y seperti terlihat pada gambar 8.1. Setiap kotak dalam jaring tersebut memiliki lebar ∆ x dan ∆ y . Oleh karena itu, panjang variabel bebas x setelah langkah ke i dinyatakan oleh xi = i ( ∆ x )

i = 0,1,......, N x

(8-5)

dan panjang variabel y setelah langkah ke j adalah yj = j( ∆ t)

j = 0,1,....., N t

(8-6)

Dengan menggunakan titik-titik jaring pada gambar 8.1, diferensial orde pertama dan

Δy

Δy

kedua dapat didekati oleh:

Persamaan Diferensial Parsial Δx

Δx

Gambar 8.1. Jaring titik-titik hitungan pada pendekatan beda hingga dengan variabel bebas x dan y.

Bab VI

Supardi, M.Si

∂ u ui + 1, j − ui , j = ∂x ∆x

(8-7)

∂ u ui , j − ui − 1, j = ∂x ∆x

(8-8) (8-9)

∂ u ui + 1, j − ui − 1, j = ∂x 2∆ x ∂ 2u ui − 1, j − 2ui , j + ui + 1, j = 2 ∂ x2 ( ∆ x)

(8-10)

∂ 2u ui + 1, j + 1 − ui − 1, j + 1 − ui + 1, j − 1 + ui − 1, j − 1 = ∂ x∂ y 4( ∆ x) ( ∆ y)

(8-11)

Dalam beberapa masalah fisika dan teknik persamaan diferensial ada yang dinyatakan dalam turunan pertama terhadap waktu dan turunan kedua terhadap ruang, misalnya pada persamaan difusi. Untuk persamaan diferensial parsial yang mengandung variabel ruang dan waktu ini, pendekatan beda hingga dapat dinyatakan

Δt

dalam jaring (jaring) bidang x dan t (lihat gambar 8.2).

ui , j + 1

ui − 1, j

ui , j

ui , j + 1

Δx Δx Gambar 8.2. Jaring titik-titik hitungan pada pendekatan beda hingga dengan variabel bebas t dan x. Persamaan Diferensial Parsial

Bab VI

Supardi, M.Si

Jaring kotak yang menyatakan variabel ruang dan waktu dibagi menjadi piaspias dengan interval ruang dan waktu ∆ x dan ∆ t . Panjang variabel ruang x setelah interval ke i dinyatakan sebagai xi = i ( ∆ x )

i = 0,1,......, N x

(8-12)

Sedangkan untuk variabel waktu t setelah interval waktu ke j adalah tj = j( ∆ t)

j = 0,1,....., N t

(8-13)

Bentuk turunan pertama terhadap waktunya dapat dituliskan sebagai ∂ u ui , j + 1 − ui , j ≈ ∂t ∆t

(8-14)

Ungkapan (8-14) dapat pula dituliskan sebagai ∂ u uij + 1 − uij ≈ ∂t ∆t

(8-15)

dengan indeks bawah menyatakan harga u pada langkah waktu, dan indeks atas menunjukkan harga u pada langkah ruang. Sedangkan untuk derivatif kedua terhadap variabel ruang seperti dinyatakan pada persamaan (8-10) dapat dituliskan kembali ∂ 2u uij− 1 − 2uij + uij+ 1 = 2 ∂ x2 ( ∆ x)

(8-16)

8.1 Klasifikasi Persamaan Diferensial Parsial Persamaan diferensial parsial dibagi menjadi tiga jenis, yaitu persamaan diferensial eliptik, parabolik dan hiperbolik. Untuk membedakan ketiga jenis persamaan diferensial parsial tersebut, marilah sekarang kita meninjau sebuah persamaan diferensial parsial orde dua dalam dua variabel ruang x dan waktu t, A

∂ 2u ∂ 2u ∂ 2u ∂u ∂u  + B + C + D  x, t , u , ,  = 0 2 2 ∂x ∂ x∂ t ∂t ∂x ∂t  

(6-1)

dimana A,B dan C merupakan fungsi dari x dan t, dan D adalah fungsi dari u dan derivatif

∂u ∂u dan , serta x dan t. Kita juga akan memperkenalkan variabel baru ∂x ∂t

Persamaan Diferensial Parsial

Bab VI

Supardi, M.Si

sedemikian hingga suku-suku yang mengandung derivatif campuran akan sama dengan nol. Selanjutnya, pembedaan atas tiga klas persamaan diferensial parsial tersebut didasarkan pada harga diskriminan B 2 − 4 AC dari persamaan (8-1) tersebut. Pertama, jika kita meninjau pada suatu titik

( x0 , t0 ) dan

di titik tersebut

memenuhi syarat bahwa harga diskriminan B 2 ( x0 , t0 ) − 4 A ( x0 , t0 ) C ( x0 , t0 ) > 0

(8-2)

maka persamaan diferensial parsial tersebut dikatakan hiperbolik pada titik ( x0 , t0 ) . Selajutnya, jika persamaan tersebut hiperbolik pada seluruh titik di dalam ranah (domain) yang ditinjau, maka persamaan diferensial parsial tersebut dikatakan sebagai persamaan hiperbolik. Sebagai contoh, jika kita meninjau persamaan gelombang yang mengambil bentuk Persamaan gelombang

∂ 2u 2 ∂ 2u −c = 0 ∂ t2 ∂ x2

Dalam persamaan gelombang tersebut harga A = − c 2 , B = 0 dan C = 1 , sehingga harga diskriminannya berharga positip. Ini berarti persamaan gelombang benar-benar masuk dalam klasifikasi persamaan diferensial hiperbolik. Persamaan (8-1) tersebut memiliki dimensi ruang satu dengan c adalah kecepatan gelombang cahaya di ruang hampa. Persamaan tersebut menjelaskan dengan sederhana bahwa derivatif kedua dari penyelesaiannya berbanding lurus dengan derivatif kedua lainnya dengan konstanta kesebandingan c 2 . Kedua, Jika pada suatu titik ( x0 , t0 ) memenuhi persyaratan B 2 ( x0 , t0 ) − 4 A ( x0 , t0 ) C ( x0 , t0 ) = 0

(8-3)

maka persamaan tersebut dikatakan parabolik pada titik ( x0 , t0 ) . Dan jika di seluruh titik dipenuhi harga diskriminan (8-3), maka persamaan tersebut disebut persamaan parabolik. Contoh dari persamaan diferensial parabolik adalah persamaan difusi yang mengambil bentuk Persamaan difusi


∂u ∂ 2u −κ = 0 ∂t ∂ x2

Bab VI

Supardi, M.Si

dengan A = − κ , B = 0 dan C = 0 . Oleh sebab itu, harga deskriminannya sama dengan nol. Persamaan ini dikenal dengan persamaan panas, yang menggambarkan aliran (difusi) panas melalui sebuah penghantar. Dalam kasus ini κ

adalah

konduktivitas termal yang merupakan kebalikan dari R yang merupakan hambatan termal. Di dalam ilmu fisika persamaan diferensial yang mirip dengan persamaan difusi adalah persamaan Schroedinger yaitu, Persamaan Schroedinger

 h2 2  ∂u  − 2 m ∇ + V ( x, y , z )  u − i h ∂ t = 0  

Persamaan Schroedinger ini memegang peran penting di dalam mekanika kuantum. Ketiga, jika pada suatu titik ( x0 , t0 ) berlaku syarat B 2 ( x0 , t0 ) − 4 A ( x0 , t0 ) C ( x0 , t0 ) < 0

(8-4)

maka persamaan tersebut dikatakan eliptik pada titik ( x0 , t0 ) , dan jika di seluruh titik dipenuhi syarat tersebut, maka persamaan tersebut masuk dalam klas persamaan eliptik. Contoh dari persamaan eliptik adalah persamaan Poisson dan Laplace yang di dalam ruang dimensi dua masing-masing mengambil bentuk Persamaan Poisson

∂ 2u ∂ 2u + = S ( x, y ) ∂ x2 ∂ y2

Persamaan Laplace

∂ 2u ∂ 2u + = 0 ∂ x2 ∂ y2

Persamaan Poisson memperkenalkan sumber panas ke dalam sistem yang ditinjau. Sedangkan persamaan Laplace merupakan kasus khusus dari persamaan Poisson tanpa sumber. Disamping itu, persamaan Laplace juga bisa diturunkan dari persamaan difusi. Jika sebuah objek diisolasi dari lingkungan, maka akan dicapai distribusi suhu dalam keadaan mantap, suatu kondisi setimbang yang digambarkan oleh derivatif waktu sama dengan nol pada persamaan difusi. Keadaan mantap suatu aliran panas ditunjukkan oleh kuantitas yang sama antara panas yang keluar dan masuk suatu tampang lintang. Dari kenyataan bahwa derivatif waktu pada persamaan difusi sama dengan nol, maka diperoleh persamaan Laplace. Oleh karena tidak ada Persamaan Diferensial Parsial

Bab VI

Supardi, M.Si

variabel waktu yang gayut, maka penyelesaian untuk persamaan Laplace maupun Poisson tersebut adalah tak gayut waktu. Persamaan menarik lain yang menggambarkan persamaan eliptik dan agak mirip dengan persamaan Poisson adalah persamaan Helmholtz yaitu, ∂ 2u ∂ 2u + + λu= 0 ∂ x2 ∂ y2

Persamaan Helmotz 8.1 Persamaan Beda Hingga 8.1.1 Persamaan Hiperbolik Persamaan Gelombang

Contoh klasik dari persamaan hiperbolik adalah persamaan gelombang yang dinyatakan oleh 2 ∂ 2u 2 ∂ u = c ∂ t2 ∂ x2

(8-17)

Persamaan ini muncul dalam berbagai masalah dari elastisitas dan akustik sampai hidraulika. Oleh sebab itu, dari tiga bentuk persamaan diferensial parsial yang kita ketahui, persamaan hiperbolik merupakan persamaan yang paling banyak dikaji oleh ilmuwan komputasi. Jika persamaan gelombang (8-17) didekati menggunakan pendekatan beda hingga, maka dapat dituliskan sebagai uij + 1 − 2uij + uij − 1

( ∆ t)

2

− c2

uij+ 1 − 2uij + uij− 1

( ∆ x)

2

= 0

(8-18)

dengan uij = u ( xi , t j )

(8-19)

j+ 1 Dengan memecahkannya untuk variabel ui maka kita memperoleh

( ∆ t ) c2 u j + u j + i− 1 ) 2 ( i+ 1 ( ∆ x) 2

ui

j+ 1

=

2  ∆ t ) c2  j ( 2 1−  ui − uij − 1 2  ( ∆ x )  

(8-20)

Persamaan ini menjelaskan kepada kita bahwa apabila kita mengetahui u pada seluruh xi pada saat-saat t j dan t j − 1 , maka kita dapat menentukan harga u pada


Bab VI

Supardi, M.Si

seluruh xi pada langkah waktu berikutnya. Hal ini disebut dengan metode eksplisit. Tetapi, ada sedikit masalah pada permulaan perhitungan, karena secara umum kita tidak mengetahui harga u pada dua waktu berturut-turut. Sedangkan, kita harus mengetahui harga u ( xi , 0 ) dan derivatif ∂ u ( xi , 0 ) ∂ t di seluruh harga xi . Oleh sebab itu, dengan mengetahui ungkapan ∂ u ( xi , t ) ∂t

t= 0

ui1 − ui− 1 = 2( ∆ t)

(8-21)

atau ui− 1 = ui1 − 2 ( ∆ t )

∂ u ( xi , t ) ∂t

(8-22) t= 0

1 maka, kita dapat menyatakan ui sebagai

( ∆ t ) c 2 u 0 + u 0 +  1 − ( ∆ t ) c 2  u 0 + ∆ t ∂ u ( xi , 0 ) ( ) i− 1 ) 2 ( i+ 1 2   i ∂t 2( ∆ x) ∆ x ( )   2

u = 1 i

2

(8-23)

Persamaan Adveksi Persamaan adveksi merupakan satu-satunya persamaan di dalam dinamika fluida yang munculnya lebih sering dibandingkan persamaan difusi. Persamaan ini memerikan cara suatu besaran kekal (conserved) seperti halnya suhu potensial ataupun momentum dibawa bersama aliran udara atau air. Untuk menjelaskan secara fisika tentang masalah adveksi ini, sekarang misalnya ada seorang pengamat berdiri di suatu lapangan dengan membawa sebuah termometer. Di tempat tersebut bertiup angin dari arah barat membawa udara lebih hangat menuju ke arah timur yang bersuhu udara lebih dingin. Dalam hal ini sebut saja bahwa arah barat ke timur adalah x . Selajutnya, apa yang dilihat oleh pengamat tersebut dengan termometer yang dibawanya? Ternyata angka yang ditunjukkan oleh termometer semakin besar, yang berarti bahwa keadaan suhu di tempat tersebut semakin hangat. Hal ini disebabkan oleh pergantian udara yang terjadi di tempat tersebut, yaitu dari keadaan udara yang dingin diganti dengan udara yang lebih hangat.


Bab VI

Supardi, M.Si

Jika yang terjadi adalah bahwa angin yang berhembus ke arah pengamat tersebut tidak mengalami perubahan suhu, maka pengamat tersebut tidak dapat memberi informasi bahwa terjadi kenaikan suhu. Nah, karena kenyataannya terjadi perubahan suhu maka ada yang disebut gradien suhu. Laju perubahan suhu yang terjadi di tempat itu bergantung kepada besarnya gradien maupun laju perpindahan udara, yaitu Laju perubahan suhu = -(Laju perpindahan udara) x (Gradien suhu) Tanda minus menyatakan bahwa suhu hanya akan naik apabila gradien suhu turun, atau dengan kata lain udara akan menjadi lebih hangat jika kita bergerak ke arah x atau dari arah timur ke barat, yakni bergerak ke arah berlawanan dengan arah angin. Dalam bahasa matematika, pernyataan di atas dapat diungkapkan dalam bentuk ∂u ∂u = −c ∂t ∂x

(8-24)

dengan u menyatakan suhu potensial yang merupakan besaran kekal yang dalam hal ini merupakan variabel yang diadveksi. Dalam kaitannya dengan masalah ini, maka kita hanya akan membahas untuk harga c konstan. Penyelesaian umum untuk persamaan (8-24) adalah u = F ( x − ct ) dengan F merupakan fungsi sembarang bernilai tunggal.

Gambar 8.3. Angin bertiup dari arah barat ke timur membawa udara hangat Persamaan Diferensial Parsial

(8-25)

Bab VI

Supardi, M.Si

Persamaan adveksi diatas merupakan contoh yang sangat bagus bahwa antara pendekatan numerik dengan analitis tidak selalu menemukan hasil yang sama. Di dalam pasal ini kita akan membahas beberapa pendekatan numerik yang dapat digunakan untuk mendekati persamaan (8-24) tersebut dan setiap metode akan kita kaji stabilitas dan akurasinya Metode FTCS (Forward-Time Centered-Space) Untuk menyelesaikan persamaan (8-24) kita akan mengimplementasikan sebuah metode dengan menggunakan pendekatan beda terpusat (metode Leap-Frog) untuk derivatif ruangnya dan metode Euler maju untuk derivatif waktunya. u nj + 1 − u nj ∆t

 ( u nj + 1 − u nj − 1 )  2 + O( ∆ t) + c  + O ( ∆ x)  = 0 2∆ x  

(8-26)

atau u nj + 1 ≈ u nj − dimana indeks bawah

j

c∆ t n u j + 1 − u nj − 1 ) ( 2∆ x menyatakan langkah ruang dan indeks atas

(8-27) n

menyatakan

langkah waktu. Dengan menggunakan analogi terhadap pembahasan tentang metode Euler dan metode Leap-Frog pada bab yang lalu, maka kita dapat menyimpulkan bahwa ketelitian untuk metode ini adalah orde pertama untuk t -nya dan orde kedua untuk x , Pendekatan beda hingga untuk persamaan adveksi (8-26) inilah yang disebut dengan forward in time, centered in space atau lebih dikenal dengan metode FTCS. Pertanyaan selanjutnya apakah metode ini stabil saat mendekati persamaan adveksi tersebut? Untuk mengetahui apakah metode yang kita gunakan untuk mendekati persamaan tersebut stabil atau tidak, maka kita perlu melakukan uji kestabilan dengan menggunakan analisa stabilitas Von Neuman. Ide dari bentuk analisis kestabilan ini, kita dapat membayangkan bahwa koefisien-koefisien dari persamaan beda berubah Persamaan Diferensial Parsial

Bab VI

Supardi, M.Si

sangat lambat ketika diperlakukan sebagai konstanta dalam ruang dan waktu. Dalam kasus demikian, penyelesaian bebasnya atau swamode dari persamaan beda mengambil bentuk u nj = ξ n exp ( ikj ∆ x )

(8-28)

dengan k menyatakan bilangan gelombang ruang real yang dapat berharga sembarang, sedangkan ξ = ξ ( k ) adalah bilangan komplek yang bergantung pada k. Jika kita mensubstitusikan persamaan (8-28) ke persamaan hampiran (8-27), maka dengan mudah diperoleh

ξ = 1− i

c∆ t sin ( k ∆ x ) 2∆ x

(8-29)

Dari persamaan (8-29) dapat diketahui modulus dari ξ yaitu

ξ

2

Persamaan

 c∆ t  = 1+  sin ( k ∆ x )   2∆ x  (8-30)

memberi

2

(8-30) arti

bahwa

penguatan

(amplification)

penyelesaiannya berhrga ≥ 1 , ini berarti bahwa metode FTCS tidak stabil mutlak untuk mendekati persamaan adveksi. Skema untuk metode FTCS dapat diilustrasikan seperti gambar 8.4

Gambar 8.4. gambaran tentang metode FTCS. Dalam gambar (8-4) tersebut bulatan kosong menggambarkan titik baru yang akan ditentukan nilainya, sedangkan bulatan hitam merupakan harga-harga fungsi yang sudah diketahui yang akan digunakan untuk memperoleh penyelesaian pada bulatan kosong. Garis sambung


menghubungkan antara titik-titik yang akan

Bab VI

Supardi, M.Si

digunakan untuk menghitung derivatif ruang, sedangkan garis putus-putus menghubungkan titik-titik yang akan digunakan untuk menghitung derivatif waktu Metode BTCS (Backward-Time Centered-Space) Dengan menggunakan pendekatan beda mundur untuk langkah waktunya dan beda terpusat untuk langkah ruangnya, maka persamaan adveksi dapat didekati dengan u nj − u nj + 1 ∆t

 ( u nj ++ 11 − u nj −+ 11 )  2 + O( ∆ t) + c  + O ( ∆ x)  = 0 2∆ x  

(8-31)

atau dapat disusun kembali menjadi u nj ≈ u nj + 1 −

c∆ t n + 1 n + 1 ( u j+ 1 − u j− 1 ) 2∆ x

(8-32)

Penggunaan analisa stabilitas Von Nouman pada pendekatan BTCS untuk persamaan adveksi ini menghasilkan 1= ξ +

c∆ t ξ eik ∆ x − ξ e − ik ∆ x 2∆ x

(

)

(8-33)

atau

ξ =

1 1+ i

c∆ t sin ( k ∆ x ) 2∆ x

(8-34)

Persamaan (8-34) menunjukkan bahwa faktor penguatannya adalah

ξ =

1 c∆ t 1+ sin ( k ∆ x ) 2∆ x

≤1

(8-35)

yang berarti, skema (8-31)) adalah stabil mutlak. Metode Centered-Time Centered-Space (CTCS) Untuk persamaan adveksi, penggunaan metode Euler maju untuk langkah waktu (forward-time) tidak stabil mutlak, apakah ini berarti dengan menggunakan pendekatan beda terpusat (centered-space) akan stabil? Untuk menjawab pertanyaan ini, marilah kita lakukan pendekatan persamaan adveksi tersebut dengan skema CTCS ini. Persamaan Diferensial Parsial

Bab VI

Supardi, M.Si

Dengan menggunakan skema CTCS, maka persamaan adveksi dapat didekati menjadi u nj + 1 − u nj − 1 2∆ t

+ O( ∆ t) − 2

u nj + 1 − u nj − 1 2∆ x

+ O ( ∆ x) = 0 2

(8-36)

Persamaan (8-36) dapat disusun kembali menjadi bentuk u nj + 1 ≈ u nj − 1 −

c∆ t n u j + 1 − u nj − 1 ) ( ∆x

(8-37)

Stabilitas Kita dapat mengetes stabilitas dari skema ini dengan analisa stabilitas Von Nouman. Dengan mensubstitusi mode Fourier adveksi yang didefinisikan (8-28) pada persamaan (8-37) maka diperoleh

ξ 2 = 1 − iξ

c∆ t sin ( k ∆ x ) ∆x

(8-38)

Persamaan (8-38) merupakan persamaan kuadrat dalam ξ , sehingga harga-harga untuk ξ dapat dinyatakan oleh

ξ 12 =

c∆ t −i sin ( k ∆ x ) ± ∆x

2

 c∆ t  − sin ( k ∆ x )  + 4  ∆x  2

(8-39)

Modulus dari masing-masing akar adalah 1, sedangkan syarat stabil adalah ξ

2

≤ 1,

ini berarti bahwa metode CTCS stabil untuk menyelesaikan persamaan adveksi. 8.7 Metode Lax Metode

Lax

merupakan

sebuah

metode

yang

dimaksudkan

untuk

memodifikasi metode FTCS dari sisi perbaikan terhadap stabilitasnya. Caranya adalah j dengan mengganti un dalam derivatif waktu dengan rerata ruangnya

u nj →

1 n ( u j + 1 + u nj − 1 ) 2

(8-40)

sehingga persamaan adveksi menjadi u nj + 1 =

1 n c∆ t n u j + 1 + u nj − 1 ) − u j + 1 − u nj − 1 ) ( ( 2 2∆ x


(8-42)

Bab VI

Supardi, M.Si

Gambar 8.5. Deskripsi untuk skema beda Lax Dengan mensubstitusi bentuk mode Fourier ke persamaan (8-28) ke persamaan beda (8-42) diperoleh

ξ = cos k ∆ x − i

c∆ t sin k ∆ x ∆x

(8-43)

Modulus dari ξ adalah

ξ

2

2

 c∆ t  2 = cos ( k ∆ x ) +   sin ( k ∆ x )  ∆x  2

(8-44)

Pernyataan (8-44) mengisyaratkan kepada kita bahwa metode Lax stabil untuk c∆ t c∆ t ≤ 1 . Untuk harga < 1 faktor penguatannya berkurang. Faktor penguatan ini ∆x ∆x dinyatakan oleh 2

ξ = Untuk harga

 c∆ t  2 cos 2 ( k ∆ x ) +   sin ( k ∆ x )  ∆x 

(8-45)

υ∆t = 1 , penyelesaiannya adalah eksak karena faktor penguatannya ∆x

berharga 1 atau tidak mengalami penguatan, sehingga u nj + 1 = u nj − 1 Kriteria stabilitas

(8-46)

c∆ t ≤ 1 dikenal dengan syarat Courant. Secara intuitif, syarat ∆x

stabilitas ini dapat dideskripsikan seperi pada gambar (8.6). Gambar tersebut n+ 1 menerangkan bahwa kuantitas u j dalam persamaan (8-42) dapat diketahui setelah

diperoleh informasi titik-titik j − 1 dan j + 1 pada saat n . Dengan kata lain, x j − 1 dan


Gambar 8.6 Daerah dibawah garis putus-putus secara fisis adalah menurut

Bab VI

Supardi, M.Si

x j + 1 merupakan batas yang memungkinkan untuk memberikan informasi pada n+ 1 besaran u j .

n Hasil yang mengagumkan pada pendekatan Lax adalah bahwa penggantian u j

dengan reratanya seperti terlihat pada ungkapan (8-41) dapat menstabilkan skema FTCS. Skema Lax pada (8-42) selajutnya dapat ditampilkan dalam bentuk u nj + 1 − u nj ∆t

 u nj + 1 − u nj − 1  1  u nj + 1 − 2u nj + u nj − 1  = − c +   2∆ x  2  ∆t    

(8-47)

yang merupakan representasi dari metode FTCS ∂u ∂ u ( ∆ x ) ∂ 2u = −c + ∂t ∂x 2∆ t ∂ x 2 2

(8-48)

Dalam persamaan (8-48) ini, kita memiliki suku difusi. Oleh sebab itu, skema Lax ini dikatakan memiliki disipasi numerik. 8.8 Skema Lax-Wendroff Skema Wendroff merupakan metode dengan akurasi orde kedua terhadap n+ 1 2 waktu. Jika kita mendefinisikan suatu harga intermediet u j + 1 2 pada langkah waktu

tn + 1 2 dan langkah ruang x j + 1 2 . Jika ini dihitung dengan menggunakan metode Lax, maka akan diperoleh u nj ++ 11 22 =

1 n ∆t u j + 1 + u nj − Fjn+ 1 − Fjn 2 2∆ x

(

)

(

)

(8-49)

n+ 1 Sedangkan, harga terbaru untuk u j dapat dihitung dengan pernyataan terpusat

sebagai Persamaan Diferensial Parsial

Bab VI

Supardi, M.Si ∆ t n+ 1 2 Fj + 1 2 − Fjn−+1122 ∆x

(

u nj + 1 = u nj −

)

(8-50)

Gambar 8.7. Titik-titik jaring pada skema LaxWendroff

Gambar 8.8. Deskripsi skema Lax-Wendorf

Selanjutnya, kita akan mengkaji stabilitas dari metode ini untuk persamaan adveksi dengan mensubstitusi F = cu . Dengan mensubstitusi pernyataan (8-49) ke ungkapan (8-50), maka diperoleh u nj + 1 = u nj − −

c∆ t  1 n 1 c∆ t n n u + u − u j + 1 − u nj j + 1 j  ∆x 2 2 ∆x

(

)

(

1 n 1 c∆ t n  u j + u nj − 1 + u j − u nj − 1  2 2 ∆x 

(

)

(

)

) (8-51)

Dengan menggunakan uji stabilitas Von Nouman, maka dengan mudah diperoleh


Bab VI

Supardi, M.Si 2

c∆ t  c∆ t  ξ = 1− i sin k ∆ x −   ( 1 − cos k ∆ x ) ∆x  ∆x 

(8-52)

Harga modulus dari ξ adalah 2

ξ

2

ξ

2

2   c∆ t  2   c∆ t  =  1−  1 − cos k ∆ x )  +  sin k ∆ x    ∆ x  (   ∆x   

(8-53)

atau  c∆ t  = 1−    ∆x 

2

  1 − 

 c∆ t     ∆x 

2

 2  ( 1 − cos k ∆ x ) 

Kriteria stabilitas yang harus dipenuhi adalah ξ

2

(8-54)

≤ 1 , hal ini mensyaratkan harga

2

 c∆ t    ≤ 1 atau lebih dikenal sebagai kriteria Courant.  ∆x  8.1.2 Persamaan Parabolik Persamaan difusi, konduksi panas dan persamaan Schroedinger gayut waktu merupakan contoh dari persamaan diferensial parabolik. Persamaan parabolik memilki kemiripan dengan persamaan hiperbolik yakni batasnya yang terbuka. Di dalam Geofisika, persamaan difusi merupakan salah satu persamaan yang sangat penting yang muncul dalam berbagai konteks yang berbeda-beda. Di bawah ini diberikan bebarapa contoh persamaan diferensial parabolik yang dinyatakan dalam ungkapan matematis a. Persamaan netron transien dalam ruang satu dimensi

ρc

∂ 2T ( x, t ) ∂T = k + Q ( x) ∂t ∂ x2

b. Persamaan konduksi panas transien dalam ruang satu dimensi 1 ∂ ∂ 2ψ ψ ( x, t ) = D 2 ∑ a ψ + υ ∑ f ψ + S υ ∂t ∂x dengan ψ menyatakan fluks netron. c. Persamaan difusi untuk transpot konvektif spesies kimia


Bab VI

Supardi, M.Si ∂ ∂ ∂2 ψ = − u ( x) ψ + D 2 ψ ∂t ∂t ∂t dengan ψ menyatakan rapat fluks spesies kimia, u ( x ) adalah kecepatan aliran dan D adalah konstanta difusi.

8.1.2.1 Metode Eksplisit (Euler Maju) Marilah kita ditinjau sebuah persamaan difusi yang mengambil bentuk ∂u ∂ 2u −κ = 0 ∂t ∂ x2

(8-55)

Dengan mengimpementasikan metode Euler maju untuk derivatif waktu seperti yang telah kita bahas pada bab persamaan diferensial biasa yang lalu, serta menggunakan pendekatan derivatif orde kedua terpusat pada turunan kedua terhadap variabel ruangnya, maka diskritisasi terhadap ungkapan (8-55) tersebut mengambil bentuk u nj + 1 − u nj ∆t

=κ

(u

n j+ 1

− 2u nj + u nj − 1 )

( ∆ x)

2

(8-56)

atau dapat dituliskan kembali sebagai u nj + 1 = u nj +

κ∆t

( ∆ x)

2

(u

n j+ 1

− 2u nj + u nj − 1 )

(8-57)

n Skema ini disebut sebagai metode eksplisit, karena jika ui diketahui untuk seluruh tn n+ 1 pada titik-titik jaring, maka kita dapat menghitung ui pada waktu tn + 1 tanpa

menyelesaikan melalui persamaan simultan. Deskripsi skema ini dapat dilihat pada gambar 8.9.

Gambar 8.9 Deskripsi metode eksplisit pada persamaan difusi


Bab VI

Supardi, M.Si

Apabila pendekatan penyelesaian persamaan difusi (8-57) dilakukan uji stabilitas menggunakan prosedur analisa stabilitas Von Nueman, maka dengan mudah dapat diperoleh bahwa

ξ = 1+ 2

k∆ t

( ∆ x)

2

( cos ( kx ) − 1)

(8-58)

1  sin 2  k ∆ x  2 

(8-59)

atau

ξ = 1− 4

κ∆t

( ∆ x)

2

Dari hasil analisa stabilitas dapat ketahui bahwa metode yang kita gunakan untuk mendekati persamaan difusi tersebut stabil karena syarat stabil ξ ≤ 1 dipenuhi. Metode Implisit (Euler Mundur) Untuk memberikan gambaran tentang pendekatan metode implisit pada persamaan difusi yang kita miliki, sekarang marilah kita mengingat kembali tentang kemungkinan pendekatan persamaan tersebut dengan beda mundur. Jika persamaan difusi tersebut kita dekati dengan beda mundur, maka diperoleh uin − uin − 1 u n − 2uin + uin− 1 = κ i+ 1 2 ∆t ( ∆ x)

(8-60)

yang dapat disusun kembali menjadi ungkapan uin −

κ∆t

( ∆ x)

2

(u

n i+ 1

)

− 2uin + uin− 1 = uin − 1

(8-61)

Ungkapan (8-61) sebenarnya mengikuti suatu perjanjian, bahwa kuantitas yang belum diketahui harganya ditempatkan di ruas kiri, sedangkan besaran yang sudah diketahui ditempatkan diruas kanan. Dalam kasus ini, harga-harga u pada langkah waktu n dianggap tidak dketahui, harga-harga yang diketahui adalah pada langkah waktu ke n − 1 . Deskripsi skema implisit ini dapat dilihat pada gambar 8.10.

Persamaan Diferensial Parsial Gambar 8.10 Deskripsi metode implisit pada persamaan difusi

Bab VI

Supardi, M.Si

Dengan mengambil

κ∆t

α ≡

( ∆ x)

(8-62)

2

maka untuk setiap titik ruang x j dengan j = 1, 2,3,..., N − 1 , kita memperoleh −α ψ

n i− 1

+ ( 1 + 2α ) ψ

n i

−αψ

n i+ 1

=ψ

n− 1 i

(8-63)

Jika syarat batas pada ujung-ujungnya diberikan yaitu u0 dan u N , maka kita persamaan (8-63) dapat ditampilkan dalam bentuk persamaan simultan linier sebagai berikut AgΨ

n

= Ψ

(8-64)

n− 1

dengan  1 −α   0 A=   .  .   .

0 1 + 2α . . . .

0 −α . . . .

. 0 . . −α 0

. . 0 . 1 + 2α 0

0 0 . . −α 1

        

(8-65)

Kita juga akan menggunakan analisa stabilitas Von Nouman untuk meyakinkan apakah skema implisit ini stabil atau tidak stabil. Jika kita mensubstitusikan mode Fourier ke persamaan (8-61), maka dengan mudah diperoleh 1−

κ∆t

( ∆ x)

2

( cos k ∆ x − 2 ) = ξ − 1

(8-66)

atau dapat disusun kembali menjadi

ξ =

1 1+

κ∆t

( ∆ x)

2

1 sin k ∆ x 2


(8-67)

Bab VI

Supardi, M.Si

Faktor penguatan yang memiliki bentuk semacam ini, tentunya harus berharga ≤ 1. Ini menunjukkan bahwa skema implisit yang kita gunakan untuk mendekati persamaan difusi adalah stabil mutlak. 8.1.2.2 Metode Dufort-Frankle Metode ini merupakan salah satu dari beberapa metode yang digunakan untuk mengatasi masalah stabilitas yang ditemukan pada metode Euler maju atau FTCS. Metode Dufort-Frankle merupakan satu teknik yang memanfaatkan stabilitas tak bersyarat dari metode intrinsic untuk persamaan diferensial sederhana. Selanjutnya kita dapat memodifikasi persamaan (8-61) menggunakan metode Dufort-Frankle sebagai berikut u nj + 1 = u nj − 1 −

2κ ∆ t

( ∆ x)

2

 u nj + 1 − ( u nj + 1 + u nj − 1 ) + u nj − 1   

(8-68)

Gambar 8.11. Deskripsi metode Dufort-Frankle Jika diambil β =

2κ ∆ t

( ∆ x)

2

, maka persamaan (8-68) dapat disusun kembali menjadi

bentuk u nj + 1 =

1 − α n− 1 α uj − ( u nj + 1 + u nj − 1 ) 1+ α 1+ α

(8-69)

Pengujian stabilitas terhadap pendekatan Dufort-Frankle menggunakan analisa Von Nouman memunculkan persamaan kuadrat dalam ξ , hal ini dikarenakan


Bab VI

Supardi, M.Si

munculnya tiga pangkat konskutif pada ξ ketika prosedur Von Nueman disubstitusi ke dalam persamaan tersebut. Persamaan kuadrat tersebut adalah

ξ 2+ξ

2α 1− α cos k ∆ x − = 0 1+ α 1+ α

(8-70)

Selanjutnya persamaan (8-70) memiliki dua penyelesaian yaitu

(

1 α cos k ∆ x ± 1 − α 2 sin 2 k ∆ x 1+ α

ξ =

)

(8-71)

Untuk mengetahui kestabilan skema ini, maka kita dapat mengecek bagaimana modulus dari ξ tersebut. Dengan menganggap α 2 sin 2 k ∆ x ≥ 1 dan α 2 sin 2 k ∆ x < 1 , maka kita akan memperoleh bahwa ξ

2

≤ 1 . Ini menunjukkan bahwa skema Dufort-

Frankle tersebut stabil mutlak. Metode Cranck-Nicolson Pendekatan metode Cranck-Nicolson untuk menyelesaikan persamaan diferensial parabolik didasarkan pada metode Euler termodifikasi seperti yang telah dibahas pada bab yang lalu. Dengan menggunakan metode ini, maka pendekatan pada persamaan difusi selanjutnya dapat ditulis kembali menjadi

ψ

n+ 1 i

−ψ ∆t

n i

=

κ 2( ∆ x)

2

(

n+ 1 i− 1

− 2ψ

n+ 1 i

+ψ

n+ 1 i+ 1

(

n+ 1 i− 1

− 2ψ

n+ 1 i

+ψ

n+ 1 i+ 1

 ψ 

) + (ψ

n i− 1

− 2ψ

n i

+ψ

n i+ 1

) 

(8-72)

) + (ψ

n i− 1

− 2ψ

n i

+ψ

n i+ 1

) 

(8-73)

atau

ψ

n+ 1 i

=ψ

n i

+

κ∆t 2( ∆ x)

2

 ψ 

Gambar 8.12. Deskripsi skema Cranck-Nicolson


Bab VI

Supardi, M.Si

Dengan mendefinisikan γ =

κ∆t 2( ∆ x)

, maka ungkapan (8-73) juga dapat dinyatakan

2

dalam bentuk persamaan simultan sebagai berikut −γψ

n+ 1 j− 1

+ ( 1 + 2γ ) ψ

n+ 1 j

− γψ

n+ 1 j+ 1

= γψ

n j− 1

+ ( 1 − 2γ ) ψ

n j

+ γψ

n j+ 1

(8-74)

atau Aψ g

n+ 1

=B ψ g

n

(8-75)

dengan matriks A dan B didefinisikan sebagai  1 −γ   0 A=   .  .   .

0 1 + 2γ . . . .

0 −γ . . . .

. 0 . . −γ 0

. . 0 . 1 + 2γ 0

0 0 . . −γ 1

        

(8-76)

0 0  .  . γ  1 

(8-77)

dan 0 1  γ 1 − 2γ  0 . B=  . . . .  .  .

0 . . γ 0 . . . 0 . . . . γ 1 − 2γ . 0 0

Dengan menggunakan analisa stabilitas Von Nouman seperti yang telah kita terapkan pada metode-metode sebelumnya, maka diperoleh faktor penguatannya sebesar

ξ =

1 − 2γ sin 2 ( k ∆ x 2 )

1 + 2γ sin 2 ( k ∆ x 2 )

(8-78)

Faktor penguatan tersebut menunjukkan bahwa harganya selalu ≤ 1 . Ini menunjukkan bahwa skema ini stabil mutlak. Lebih lanjut lagi, karena pendekatan beda yang digunakan dalam metode ini adalah metode Euler termodifikasi, maka ketelitian metode ini lebih tinggi dibanding metode Euler maju ataupun mundur.


Bab VI

Supardi, M.Si

Contoh penggunaan skema Cranck-Nicolsan adalah pada penyelesaian persamaan Schroedinger. 8.13 Persamaan Schroedinger Jika kita mengkaji secara serius ilmu fisika, maka kadang-kadang kita menemukan suatu masalah yang mengandung kendala (constraint), sebagai contoh persamaan Scrhoedinger gayut waktu di dalam Mekanika Persamaan ini termasuk ke dalam persamaan diferensial parabolik untuk evolusi besaran kompleks ψ . Untuk persamaan diferensial parsial yang memerikan hamburan paket gelombang yang disebabkan oleh potensial V ( x ) dalam ruang 1D, maka persamaannya memiliki bentuk i

∂ψ h2 ∂ 2ψ = − + V ( x) ψ ∂t 2m ∂ x 2

(8-79)

Jika kita menggunakan satuan universal, sedemikian hingga konstanta Planck

 = 1 dan massa partikel m = 1 / 2 , maka persamaan Schroedinger (8-79) akan mengmbil bentuk i

∂ψ ∂ 2ψ = − + V ( x )ψ ∂t ∂ x2

(8-80)

Pengenaan syarat batas untuk masalah di atas adalah harga ψ pada saat awal atau ψ ( x, t = 0 ) bersama dengan x → ± ∞ yaitu ψ → 0 . Selanjutnya langkah diskritisasi

untuk persamaan gelombang (8-71) dapat dinyatakan dalam bentuk ψ i 

ψ − ψ nj   = − ∆t  

n+ 1 j

n+ 1 j+ 1

− 2ψ

n+ 1 j 2

( ∆ x)

+ 2ψ

n+ 1 j− 1

  + V jψ 

n+ 1 j

(8-81)

Skema yang ditunjukkan pada persamaan beda (8-81) menggunakan skema implisit atau metode BTCS. Oleh sebab itu, factor penguatannya adalah

ξ =

1  4∆ t  2  k∆ x  1+ i  sin + V ∆ t    j 2  2   ( ∆ x ) 

atau


(8-82)

Bab VI

Supardi, M.Si

ξ

Dengan harga ξ

2

1

=

2

 4∆ t   k∆ x  1+  sin 2   + Vj∆ t  2  2   ( ∆ x ) 

(8-83)

di atas menunjukkan bahwa skema ini stabil mutlak. Sayangnya,

skema ini tidak uniter. Mengapa harus uniter? Hal ini disebabkan oleh suatu syarat bahwa probabilitas total suatu partikel ditemukan dalam suatu range daerah yang terbentang dari − ∞ sampai ∞ adalah satu. ∞

∫

ψ

2

dx = 1

(8-84)

−∞

Persamaan (8-84) mensyaratkan fungsi gelombang awal ψ ( x, 0 ) ternormalisir. Jika ungkapan persamaan Schroedinger (8-80) dinyatakan dalam bentuk i

∂ψ = Hψ ∂t

(8-85)

dengan H adalah operator hamiltonian yang mengambil bentuk ∂2 H = − 2 + V ( x) ∂x

(8-86)

maka penyelesaian persamaan (8-85) tersebut secara analitik adalah

ψ ( x, t ) = e− iHtψ ( x, 0 )

(8-87)

Implementasi algoritma FTCS untuk mendekati persamaan (8-87) berbentuk

ψ

n+ 1 j

= ( 1 − iH ∆ t ) ψ

n j

(8-88)

dimana H dinyatakan oleh pendekatan beda hingga terpusat dalam x . Sedangkan, penggunaan skema implisit BTCS akan mengambil bentuk berbeda yaitu

ψ

n+ 1 j

= ( 1 + iH ∆ t ) ψ −1

n j

(8-89)

Dua metode yang digunakan di atas memiliki akurasi orde pertama dalam waktu, seperti telah dibahas di depan. Dengan kenyataan bahwa metode eksplisit maupun implisit bukan metode yang baik untuk menyelesaikan persamaan Schroedinger gayut waktu ini, maka kita


Bab VI

Supardi, M.Si

akan menggunakan bentuk Cayleys untuk menyatakan wakilan beda hingga e − iHt yang memiliki akurasi orde dua dan uniter yaitu e − iHt

1 iH ∆ t 2 ; 1 1 + iH ∆ t 2 1−

(8-90)

dengan kata lain, 1    1 + iH ∆ t  ψ 2  

n+ 1 j

1   =  1 − iH ∆ t  ψ 2  

n j

(8-91)

Selanjutnya dari persamaan (8-91), maka kita memiliki sistem tridiagonal. Skema tersebut adalah stabil, uniter dan memiliki akurasi orde kedua. Nah cara ini disebut sebagai metode Crank-Nicolson. Contoh source code untuk menyelesaikan persamaan difusi Program Difusi Integer*4 maxn, maxnplot parameter( maxn = 300, maxnplot = 500 ) integer*4 n, i, j, iplot, nlangkah, plot_langkah, nplot, ilangkah real*8 tau, l, h, kappa, koef, tt(maxn), tt_baru(maxn) real*8 xplot(maxn), tplot(maxnplot), ttplot(maxn,maxnplot) C initialisasi parameter (langkah waktu, pias, dll). write(*,*) ‘masukkan langkah waktu: ' read(*,*) tau write(*,*) ‘masukkan jumlah jaring: ' read(*,*) n l = 1. h = l/(n-1) kappa = 1. koef = kappa*tau/h**2 if( koef .lt. 0.5 ) then write(*,*) 'penyelesaian diharapkan stabil' else write(*,*) 'warning: apakah penyelesaian diharapkan tidak stabil’ endif C set syarat awal dan syarat batas. do i=1,n tt(i) = 0.0 Persamaan Diferensial Parsial

Bab VI

Supardi, M.Si tt_baru(i) = 0.0 enddo tt(n/2) = 1/h iplot = 1 nlangkah = 300 plot_langkah = 6 nplot = nlangkah/plot_langkah + 1 do i=1,n xplot(i) = (i-1)*h - l/2 enddo do ilangkah=1,nlangkah do i=2,(n-1) tt_baru(i) = tt(i) + koef*(tt(i+1) + tt(i-1) - 2*tt(i)) enddo do i=2,(n-1) tt(i) = tt_baru(i) enddo if( mod(ilangkah,plot_langkah) .lt. 1 ) then do i=1,n ttplot(i,iplot) = tt(i) enddo tplot(iplot) = ilangkah*tau iplot = iplot+1 endif enddo nplot = iplot-1

open(11,file='tplot.txt',status='unknown') open(12,file='xplot.txt',status='unknown') open(13,file='ttplot.txt',status='unknown') do i=1,nplot write(11,*) tplot(i) enddo do i=1,n write(12,*) xplot(i) do j=1,(nplot-1) write(13,1001) ttplot(i,j) enddo write(13,*) ttplot(i,nplot) enddo 1001 format(e12.6,', ',$) stop end 8.2 Persamaan Eliptik Persamaan Diferensial Parsial

Bab VI

Supardi, M.Si

Contoh umum dari persamaan diferensial eliptik adalah persamaan Poisson yang berbentuk ∂ 2u ∂ 2u + = − ρ ( x, y ) ∂ x2 ∂ y2

(8-92)

Jika ρ ( x, y ) = 0 , maka disebut persamaan Laplace yang berbentuk ∂ 2u ∂ 2u + = 0 ∂ x2 ∂ y2

(8-93)

Untuk menyelesaikan persamaan eliptik dibutuhkan syarat batas di ujungujungnya. Oleh sebab itu penyelesaian persamaan eliptik masuk dalam kategori masalah nilai batas. Metode

penyelesaian

numerik

untuk

persamaan

diferensial

eliptik

diklasifikasikan dalam dua kategori, yaitu metode beda hingga dan elemen hingga. Tetapi dalam pasal ini kita hanya akan menggunakan metode beda hingga untuk menangani persamaan ini. Metode beda hingga diturunkan dari jaring kotak. Penggunaan metode ini untuk menyelesaikan masalah diferensial eliptik memiliki banyak keuntungan. Adapun keuntungan metode elemen hingga diantaranya adalah bahwa persamaan diskritnya tidak terganggu oleh bentuk geometri yang rumit, sehingga metode ini fleksibel untuk diterapkan dalam bentuk geometri apapun. Namun akhir-akhir ini, metode beda hingga juga telah dikembangkan untuk mengatasi masalah geometri ini yaitu dengan cara transformasi koordinat. Persamaan Beda dalam Geometri Rectangular Dalam pasal ini kita tidak akan membahas metode beda hingga dalam geometri yang rumit, tetapi kita hanya akan membahas metode tersebut di dalam gometri kotak saja. Untuk memudahkan pemahaman kita tentang metode ini, sekarang marilah kita tinjau sebuah persamaan Laplace dalam koordinat kartesan seperti terlihat pada persamaan (8-93). Untuk mempermudah pemahaman kita tentang masalah yang kita bahas ini, sekarang ditinjau untuk domain 0 ≤ x ≤ Lx dan 0 ≤ y ≤ Ly seperti terlihat pada gambar 8.10. Syarat batas yang dikenakan pada sisi-sisinya adalah


Bab VI

Supardi, M.Si

Batas kiri

∂u = 0 ∂x

(syarat batas Neumann)

Batas kanan

u= 0

(syarat Dirichlet)

Batas atas

u= 0

Batas bawah

∂u = 0 ∂y

Untuk menurunkan persamaan beda hingga pada persamaan Laplace, maka kita perlu membuat jaring pada kotak tersebut. Jika kita mengasumsikan bahwa lebar pias ∆ = ∆ x = ∆ y , maka persamaan Poisson tersebut dapat didekati dengan pendekatan beda terpusat yang mengambil bentuk 1 δ i2ui , j + δ j2ui , j = ρ i , j ∆2

(

)

(8-94)

atau secara eksplisit dapat ditunjukkan dalam bentuk deskrit 1 1 u − 2ui , j + ui − 1, j  + 2  ui , j + 1 − 2ui , j + ui , j − 1  = ρ i , j 2  i + 1, j ∆ ∆

(8-95)

4

… jmax

dengan ui , j ≡ u ( xi , y j )

∂u = 0 ∂x

y= 0

x= 0

3

u= 0

2

∆

j=1

Ly

u= 0 ∂u ∂y

L

∆

0 dalamx menangani Masalah yang =timbul persamaan beda (8-95) adalah i i=1

max

2

3

4

…

bagaimana cara memberikan perlakuan pada titik-titik di syarat batas pada sisi kiri Gambar 8.14. Persamaan Laplace dalam geometri kotak bersama dengan syarat batasnya Persamaan Diferensial Parsial

Bab VI

Supardi, M.Si

dan bawah. Sedangkan persamaan beda di sisi atas dan kanan tidak penting, karena kita sudah mengetahui harga dari u. Sekarang kita akan meninjau untuk persamaan beda di perbatasan sisi kiri. Pada gambar tersebut terlihat bahwa syarat batas di sebelah kiri adalah

Derivatif kedua dari u atau

∂ 2u ∂ x2

∂u = 0. ∂x

pada titik-titik batas sebelah kiri selanjutnya dapat kita

dekati dengan  ∂ 2u   2  ∂ x  1, j

 ∂u  ∂u   1 −   ∂ x  i+ , j  ∂ y  i, j 2 ≈ ∆ 2

(8-96)

Sekarang, bagaimana cara kita mendekati derivatif pertama u terhadap x pada suku pertama pembilang pada persamaan beda (8-96) tersebut? Lihatlah gambar 8.15. Pada gambar tersebut tampak tiga titik terdekat yang mengelilingi titik (1,j). Dengan kenyataan itu, kita memang tidak bisa menerapkan persamaan beda untuk derivatif kedua seperti terlihat pada persamaan (8-95), karena kita harus memiliki empat titik yang mengitari titik pusat. Oleh sebab itu, kita berharap dengan pendekatan beda (896) tersebut permasalahan tersebut teratasi.

i,j+1 i-1,j

i,j

i+1,j

i,j-1 Gambar 8.15. Titik (i,j) dikelilingi empat titik terdekat. Pendekatan beda untuk derivatif kedua u terhadap x dinyatakan pada persamaan (88). 1,j+1 1,j Persamaan Diferensial Parsial

2,j

1,j-1

Gambar 8.16. Titik-titik di perbatasan kiri

Bab VI

Supardi, M.Si

Selajutnya, suku pertama pembilang pada persamaan (8-96) dapat didekati dengan u2, j − u1, j  ∂u   1 ≈ ∆  ∂ x  1+ , j

(8-97)

2

Oleh karena suku kedua pembilang pada persamaan beda (8-96) berharga sama dengan nol, maka pendekatan beda

∂ 2u ∂ x2

adalah

2u2, j − 2u1, j  ∂ 2u   2 ≈ ∆  ∂ x  1, j

(8-98)

Selanjutnya, persamaan beda derivatif kedua u terhadap x di perbatasan sisi kiri mengambil bentuk 1 1  2u2, j − 2u1, j  + 2  u1, j + 1 − 2u1, j + u1, j − 1  = ρ 1, j 2  ∆ ∆

(8-99)

Dengan menggunakan cara yang sama, maka kita dapat mendekati derivatif pertama u terhadap x di perbatasan sisi bawah kotak adalah 2ui ,2 − 2ui ,1  ∂ 2u   2 ≈ ∆  ∂ y  i ,1 sehingga derivatif kedua u dapat didekati dengan persamaan beda 1 1  u − 2ui ,1 + ui − 1,1  + 2  2ui ,2 − 2ui ,1  = ρ i ,1 2  i + 1,1 ∆ ∆

I,j+1 Persamaan Diferensial Parsial i-1,j

I,j

i+1,j

Gambar 8.16. Titik-titik di perbatasan kiri

(8-100)

(8-101)

Bab VI

Supardi, M.Si

Dengan menggunakan persamaan beda () dan () untuk derivatif pertama u terhadap x, maka kita dapat menentukan pendekatan derivatif kedua dari u yang mengambil bentuk 1 1  u − u1,1  + 2  u1,2 − u1,1  = ρ 1,1 2  2,1 ∆ ∆

(8-102)

Contoh Ditinjau sebuah persamaan Laplace dalam ruang dimensi dua dengan domain ∂ 2ϕ ∂ 2ϕ + = 0 ∂ x2 ∂ y 2

0 ≤ x ≤ 8 dan 0 ≤ y ≤ 6 mengambil bentuk

dengan syarat batas yang diberikan adalah Batas kiri

∂u = 0 ∂x

u= 0

Ly

Batas kanan u = 1 Batas atas

u= 0

u= 0

∂u = 0 Batas bawah ∂y

6

∂u = 0 ∂x

y= 0

∂u = 0 ∂y

x= 0

8

4

Penyelesaian

Lx

3

Untuk menyelesaikan persamaan di atas, maka kita membuat jaring dengan 2

lebar pias sama yaitu 2. Harga titik-titik di perbatasan kotak atas dan kanan berharga

ϕ = 0 , sedangkan titik-titik di perbatasan kiri dan bawah memenuhi ∂ ϕ ∂ x = 0 dan 2 j=1

∂ ϕ ∂ y = 0 . Dengan syarat batas yang diberikan tadi, kita akan menghitung titik-titik


i=1

2

2

3

4

5

Bab VI

Supardi, M.Si

yang lain kecuali pada perbatasan atas dan kanan, karena di perbatasan ini harga ϕ sudah diketahui.

Dengan menggunakan persamaan beda hingga (8-95), (8-99) dan (8-101), maka kita dapat menuliskan persamaan simultan dalam ϕ yaitu 1. Titik (1,1)

− 4ϕ 1,1 + 2ϕ 2,1 + 2ϕ 1,2 = 0

2. Titik (2,1) ϕ 1,1 − 4ϕ 2,1 + ϕ 3,1 + 2ϕ 2,2 = 0 3. Titik (3,1) ϕ 2,1 − 4ϕ 3,1 + ϕ 4,1 + 2ϕ 3,2 = 0 4. Titik (4,1) ϕ 3,1 − 4ϕ 4,1 + 2ϕ 4,2 = − 1 5. Titik (1,2) ϕ 1,1 − 4ϕ 1,2 + 2ϕ 2,2 + ϕ 1,3 = 0 6. Titik (2,2) ϕ 2,1 + ϕ 1,2 + ϕ 3,2 − 4ϕ 2,2 + ϕ 2,3 = 0 7. Titik (3,2) ϕ 3,1 + ϕ 2,2 − 4ϕ 3,2 + ϕ 4,2 + ϕ 3,3 = 0 8. Titik (4,2) ϕ 4,1 + ϕ 3,2 − 4ϕ 4,2 + ϕ 4,3 = − 1 9. Titik (1,3) u1,2 − 4u1,3 + 2u2,3 = 0 10. Titik (2,3) u2,2 + u1,3 − 4u2,3 + u3,3 = 0 11. Titik (3,3) u3,2 + u2,3 − 4u3,3 + u4,3 = 0 12. Titik (4,3) ϕ 4,2 + ϕ 3,3 − 4ϕ 4,3 = − 1

Jika persamaan simultan di atas dinyatakan dalam bentuk matriks, maka bentuknya


Bab VI

Supardi, M.Si

−4 2 0 0 2 0 0 0 0 0 0 0   1 −4 1 0 0 2 0 0 0 0 0 0     0 1 −4 1 0 0 2 0 0 0 0 0     0 0 1 −4 0 0 0 2 0 0 0 0   1 0 0 0 −4 2 0 0 1 0 0 0     0 1 0 0 1 −4 1 0 0 1 0 0   0 0 1 0 0 1 −4 1 0 0 1 0     0 0 0 1 0 0 1 −4 0 0 0 1   0 0 0 0 1 0 0 0 −4 2 0 0     0 0 0 0 0 1 0 0 1 −4 1 0     0 0 0 0 0 0 1 0 0 1 −4 1   0 0 0 0 0 0 0 1 0 0 1 − 4 

 ϕ 11   0  ϕ   0   21     ϕ 31   0       ϕ 41   − 1  ϕ 12   0       ϕ 22   0  ϕ  =  0   32     ϕ 42   − 1 ϕ   0   13     ϕ 23   0       ϕ 33   0   ϕ 43   − 1

Jika persamaan simultan linier yang dinyatakan dalam bentuk matriks tersebut dinyakan oleh AgX = b maka untuk menemukan harga setiap elemen matriks X dapat dilakukan dengan cara X = A − 1 gb Dengan menggunakan kaidah ini, maka elemen-elemen matriks X maka kita dapat menentukan harga ϕ pada setiap titik yaitu

ϕ 11 = 0.3377, ϕ 21 = 0.3799, ϕ 31 = 0.5118 ϕ 41 = 0.7379, ϕ 12 = 0.2955, ϕ 22 = 0.3351 ϕ 32 = 0.4647, ϕ 42 = 0.7199, ϕ 13 = 0.1740 ϕ 23 = 0.2003, ϕ 33 = 0.2920, ϕ 43 = 0.5030 Metode Iteratif Jacobi Sesuai dengan namanya, ide dari metode iteratif Jacobi adalah menemukan harga setiap titik-titik dalam kotak melalui jalan iterasi hingga ditemukan harga yang optimum. Iterasi awal dimulai dengan memberikan nilai tebakan pada variabelvariabelnya. Iterasi dilakukan terus menerus hingga selisih harga elemen kini dan sebelumnya melebihi toleransi yang diberikan. Untuk lebih jelasnya, sekarang kita akan meninjau kembali persamaan Laplace seperti pada contoh 8.1 tetapi dengan syarat batas sebagai berikut Batas kiri

u= 0


Bab VI

Supardi, M.Si u= 1

Batas kanan Batas atas

u= 0

Batas bawah

u=1

Dengan syarat batas yang diberikan tersebut, maka kita dapat membentuk persamaan simultan baru sebagai berikut 4ϕ 2,2 − ϕ 1,2 − ϕ 3,2 − ϕ 2,1 − ϕ 2,3 = 0 1. Titik (2,2)

ϕ 2,2 =

atau

ϕ 3,2 + ϕ 2,3 4

4ϕ 3,2 − ϕ 2,2 − ϕ 4,2 − ϕ 3,1 − ϕ 3,3 = 0 2. Titik (3,2)

ϕ 3,2 =

atau

ϕ 2,2 + ϕ 4,2 + ϕ 3,3 4

4ϕ 4,2 − ϕ 3,2 − ϕ 5,2 − ϕ 4,1 − ϕ 4,3 = 0 3. Titik (4,2)

ϕ 4,2 =

atau

1 + ϕ 3,2 + ϕ 4,3 4

4ϕ 2,3 − ϕ 1,3 − ϕ 3,3 − ϕ 2,2 − ϕ 2,4 = 0 4. Titik (2,3)

ϕ 2,3 =

atau

1 + ϕ 2,2 + ϕ 3,3 4

4ϕ 3,3 − ϕ 2,3 − ϕ 4,3 − ϕ 3,2 − ϕ 3,4 = 0 5. Titik (3,3)

ϕ 3,3 =

atau

1 + ϕ 2,3 + ϕ 4,3 + ϕ 3,2 4

4ϕ 4,3 − ϕ 3,3 − ϕ 5,3 − ϕ 4,2 − ϕ 4,4 = 0 6. Titik (4,3)

0

ϕ 4,3 =

atau

1

1

2 + ϕ 3,3 + ϕ 4,2 4 1 1

0

1

0

1

0

0


0

0

1

Bab VI

Supardi, M.Si

Sebagai langkah awal, kita berikan tebakan awal seluruh titik sama dengan nol, kecuali pada batas-batas yang telah kita tentukan. Dari langkah ini, kita memiliki 0 0 0 0 harga-harga pada setiap titik antara lain ϕ 2,2 = 0 , ϕ 3,2 = 0 , ϕ 4,2 = 0, 25 , ϕ 2,3 = 0, 25 , 0 0 ϕ 3,3 = 0, 25 , ϕ 4,3 = 0,50 . Dengan menggunakan bahasa pemrograman, maka arga

titik-titik pada iterasi berikutnya dapat kita temukan sampai toleransi yang diberikan. Program Iterasi_jacobi dimension pa(5,4),pb(5,4) real phip character*10 fname write(*,5) read 9,fname 9 format(15a) 5 format(23x,'nama file output:',\) open(8,file=fname) c tebakan awal untuk seluruh titik diberikan sama dengan nol c kecuali pada batas-batas yang telah ditentukan c syarat batas pada titik-titik jaring adalah pa(2,4)=1. pa(3,4)=1. pa(4,4)=1. pa(5,4)=1. pa(5,2)=1. pa(5,3)=1. c pb(2,4)=1. pb(3,4)=1. pb(4,4)=1. pb(5,2)=1. pb(5,3)=1. pb(5,4)=1. c do 25 iter=1,100 write(8,90)iter do 30 i=2,4 do 40 j=2,3 pb(i,j)=(pa(i-1,j)+pa(i+1,j)+pa(i,j-1)+pa(i,j+1))/4. Persamaan Diferensial Parsial

Bab VI

Supardi, M.Si 40 continue 30 continue do 35 i=2,4 do 45 j=2,3 pa(i,j)=pb(i,j) 45 continue 35 continue

c write (8,22) pb(1,4),pb(2,4),pb(3,4),pb(4,4),pb(5,4) write (8,22) pb(1,3),pb(2,3),pb(3,3),pb(4,3),pb(5,3) write (8,22) pb(1,2),pb(2,2),pb(3,2),pb(4,2),pb(5,2) write (8,22) pb(1,1),pb(2,1),pb(3,1),pb(4,1),pb(5,1) write (*,22) pb(1,4),pb(2,4),pb(3,4),pb(4,4),pb(5,4) write (*,22) pb(1,3),pb(2,3),pb(3,3),pb(4,3),pb(5,3) write (*,22) pb(1,2),pb(2,2),pb(3,2),pb(4,2),pb(5,2) write (*,22) pb(1,1),pb(2,1),pb(3,1),pb(4,1),pb(5,1) 25 continue 90 format(i4) 22 format(5f10.6) close(8) stop end Tabel 8.1 Contoh eksekusi program iterasi Jacobi Iterasi ke-1 .000000 1.000000 .000000 .250000 .000000 .000000 .000000 .000000 Iterasi ke-4 .000000 1.000000 .000000 .406250 .000000 .148438 .000000 .000000 Iterasi ke-7 .000000 1.000000 .000000 .451660 .000000 .190125 .000000 .000000 Iterasi ke-10 .000000 1.000000 .000000 .460958 .000000 .200234 .000000 .000000 Iterasi ke-13 .000000 1.000000


1.000000 .250000 .000000 .000000

1.000000 .500000 .250000 .000000

1.000000 1.000000 1.000000 .000000

1.000000 .574219 .265625 .000000

1.000000 .738281 .480469 .000000

1.000000 1.000000 1.000000 .000000

1.000000 .634094 .330688 .000000

1.000000 .784973 .523438 .000000

1.000000 1.000000 1.000000 .000000

1.000000 .648406 .343851 .000000

1.000000 .794291 .533567 .000000

1.000000 1.000000 1.000000 .000000

1.000000

1.000000

1.000000

Bab VI

Supardi, M.Si .000000 .463182 .000000 .202281 .000000 .000000 Iterasi ke-16 .000000 1.000000 .000000 .463632 .000000 .202770 .000000 .000000 Iterasi ke-19 .000000 1.000000 .000000 .463740 .000000 .202869 .000000 .000000 Iterasi ke-22 .000000 1.000000 .000000 .463762 .000000 .202892 .000000 .000000 Iterasi ke-25 .000000 1.000000 .000000 .463767 .000000 .202897 .000000 .000000 Iterasi ke-27 .000000 1.000000 .000000 .463768 .000000 .202898 .000000 .000000 Iterasi ke-28 .000000 1.000000 .000000 .463768 .000000 .202898 .000000 .000000

.651300 .346998 .000000

.796516 .535614 .000000

1.000000 1.000000 .000000

1.000000 .651992 .347634 .000000

1.000000 .796966 .536103 .000000

1.000000 1.000000 1.000000 .000000

1.000000 .652132 .347786 .000000

1.000000 .797073 .536202 .000000

1.000000 1.000000 1.000000 .000000

1.000000 .652165 .347817 .000000

1.000000 .797095 .536226 .000000

1.000000 1.000000 1.000000 .000000

1.000000 .652172 .347824 .000000

1.000000 .797100 .536230 .000000

1.000000 1.000000 1.000000 .000000

1.000000 .652173 .347825 .000000

1.000000 .797101 .536231 .000000

1.000000 1.000000 1.000000 .000000

1.000000 .652174 .347826 .000000

1.000000 .797101 .536232 .000000

1.000000 1.000000 1.000000 .000000

8.10 Metode Relaksasi Konsep dari metode relaksasi didasarkan pada suatu ide bahwa konvergensi ke suatu penyelesaian dari pemberian terkaan awal tertentu dapat dicapai dengan cara mengulang-ulang iterasi setiap titiknya. Konsep dari iterasi berasal dari suatu ide bahwa perubahan perlahan-lahan (evolusi) terhadap waktu dapat dilihat ketika persamaan diferensial parsial eliptik dinyatakan dalam bentuk persamaan diferensial parabolik. 8.10.1 Metode RelaksasiGauss-Seidel


Bab VI

Supardi, M.Si

Metode relaksasi Gauss-Seidel telah terbukti memperoleh sukses besar dalam keberhasilannya menyelesaikan persamaan diferensial parsial eliptik. Untuk lebih jelasnya, sekarang kita akan menyatakan persamaan eliptik sebagai persamaan difusi ∂ 2u ∂ 2u + = ρ ( x, y ) ∂ x2 ∂ y2

(8-72)

menjadi ∂ u ∂ 2u ∂ 2u = + − ρ ( x, y ) ∂ t ∂ x2 ∂ y 2

(8-73)

Apabila pada t = 0 terdapat distribusi awal, maka kita dapat mengatakan bahwa bahwa ketika t → ∞ penyelesaian sudah merelaks ke arah keadaan setimbang. Saat t → ∞ tersebut, maka dipenuhi ∂ u / ∂ t → 0 . Jika persamaan (8-73) kita lakukan diskritisasi menggunakan metode FTCS, maka ungkapan tersebut akan menjadi bentuk u nj ,+ℓ1 = u nj ,ℓ +

∆t n ( u j − 1,ℓ + u nj + 1,ℓ + u nj ,ℓ− 1 + u nj ,ℓ+ 1 − 4u nj ,ℓ ) − ρ ∆ x2

j ,ℓ

∆t

(8-74)

dengan indeks atas n mewakili variabel waktu, sedangkan indeks bawah menyatakan variabel ruang. Dengan mengingat kembali bahwa di dalam ruang 1D metode FTCS stabil 2 hanya jika dipenuhi ∆ t / ∆ ≤

1 1 2 , dan stabil dalam ruang 2D hanya jika ∆ t / ∆ ≤ , 2 4

maka ungkapan (8-74) dapat dinyatakan kembali dalam bentuk u nj ,+ℓ1 =

1 n ∆2 n n n u + u + u + u − ( j − 1,ℓ j + 1,ℓ j ,ℓ− 1 j ,ℓ+ 1 ) 4 ρ 4

j ,ℓ

(8-75)

Dari ungkapan (8-75), kita dapat menemukan harga terbaru dari u

pada

langkah ( n + 1) dengan menggunakan empat harga lama yang mengelilinginya pada langkah n dan suku sumbernya.

Prosedur menemukan harga terbaru tersebut

dilakukan dengan cara menyapu titik-titik yang diawali dari baris demi baris titik dan menghitung harga baru u dengan mengunakan ungkapan (8-75). Prosedur ini diulangulang hingga ketelitian yang diharapkan dicapai. Metode ini disebut dengan iterasi


Bab VI

Supardi, M.Si

Jacobi seperti yang telah dibahas di atas. Sayangnya, metode ini masih cukup lambat mencapai konvergen. Satu metode yang barangkali lebih baik dibandingkan dengan metode iterasi Jacobi membuat algoritma tersebut menjadi bentuk semi implisit u

n+ 1 j ,ℓ

1 n ∆2 n+ 1 n n+ 1 = ( u j + 1,ℓ + u j − 1,ℓ + u j ,ℓ+ 1 + u j ,ℓ− 1 ) − ρ 4 4

(8-76)

j ,ℓ

Dalam skema ini, harga-harga baru dari u digunakan segera setelah hargaharga tersebut ada, artinya bahwa titik-titik yang sudah ter-update akan digunakan segera dalam perhitungan untuk memperoleh harga terbaru u pada titik berikutnya. Skema yang diperlihatkan pada (8-76) tersebut dikenal dengan metode relaksasi Gauss-Seidel. Sayangnya, metode ini juga masih lambat konvergensinya. 8.10.2 Metode Over-Relaksasi Simultan Untuk memperoleh metode relaksasi lebih baik dalam hal kecepatan konvergensi, maka kita perlu mengkoreksi secara over metode Gauss-Seidel. Kita akan melakukan generalisasi terhadap skema (8-76) sehingga setiap langkah relaksasi

ϕ

j ,l

akan digantikan dengan kombinasi linier antara harga lamanya dan harga

terupdatenya. Jadi u nj ,+ℓ1 = ( 1 − ω ) u nj ,l +

ω  ( u nj + 1,ℓ + u nj −+ 1,1 ℓ + u nj ,ℓ+ 1 + u nj ,+ℓ1− 1 ) − ∆ 2 ρ 4

j ,ℓ

 

(8-76)

dimana ω merupakan parameter over relaksasi. Metode ini konvergen hanya dalam ranah 0 < ω < 2 . Untuk harga 0 < ω < 1 , maka skema (8-76) disebut dengan under relaxation , sedangkan untuk ranah 1 < ω < 2 skema tersebut dikenal dengan over relaxation. Untuk harga ω dalam ranah 1 < ω < 2 memberikan konvergensi lebih cepat dibandingkan dengan metode Gauss-Seidel. Contoh source code untuk menyelesaikan Persamaan Laplace menggunakan Iterasi Gauss-Seidel dan over relaksasi Program Laplace Integer max real omega Persamaan Diferensial Parsial

Bab VI

Supardi, M.Si

Parameter(max1=4,max2=5,omega=1.25) Real*8 x, p(max2,max1),phip Integer i, j, iter, y c membuka file output Open(8, File='laplace.dat', Status='Unknown') c sisi-sisi jaring dengan potensial konstan Do 10 i=1, max2 p(i,1)=0.0 p(i,4)=1.0 10 Continue Do 11 j=1, max1 p(1,j)=0.0 p(5,j)=1.0 11 Continue c

c c c c c c c c

algoritma iterasi Do 20 iter=1, 100 write(8,21)iter Do 30 i=2,(max2-1) Do 40 j=2,(max1-1) menentukan harga titik-titik pada jaring dengan metode Gauss-Seidel p(i,j)=0.25*(p(i+1,j) * +p(i-1,j)+p(i,j+1)+p(i,j-1)) menentukan harga titik-titik pada jaring dengan metode over relaksasi dengan parameter relaksasi omega=1.25 phip=0.25*(p(i+1,j) * +p(i-1,j)+p(i,j+1)+p(i,j-1)) p(i,j)=(1.-omega)*p(i,j)+omega*phip 40 Continue 30 Continue Write (8,22) p(1,4),p(2,4),p(3,4),p(4,4),p(5,4) Write (8,22) p(1,3),p(2,3),p(3,3),p(4,3),p(5,3) Write (8,22) p(1,2),p(2,2),p(3,2),p(4,2),p(5,2) Write (8,22) p(1,1),p(2,1),p(3,1),p(4,1),p(5,1) Write (*,22) p(1,4),p(2,4),p(3,4),p(4,4),p(5,4) Write (*,22) p(1,3),p(2,3),p(3,3),p(4,3),p(5,3) Write (*,22) p(1,2),p(2,2),p(3,2),p(4,2),p(5,2)


Bab VI

Supardi, M.Si Write (*,22) p(1,1),p(2,1),p(3,1),p(4,1),p(5,1)

20 continue 21 Format(i4) 22 Format(5f10.6) Close(8) Stop 'data tersimpan dalam laplace.dat’ End Contoh eksekusi untuk penyelesaian persamaan Laplace menggunakan metode iterasi Gauss-Seidel Iterasi ke-1 .000000 1.000000 .000000 .250000 .000000 .000000 .000000 .000000 Iterasi ke-3 .000000 1.000000 .000000 .415039 .000000 .125000 .000000 .000000 Iterasi ke-5 .000000 1.000000 .000000 .457157 .000000 .191956 .000000 .000000 Iterasi ke-7 .000000 1.000000 .000000 .462890 .000000 .201444 .000000 .000000 Iterasi ke-9 .000000 1.000000 .000000 .463652 .000000 .202706 .000000 .000000 Iterasi ke-11 .000000 1.000000 .000000 .463753 .000000 .202873 .000000 .000000 Iterasi ke-13 .000000 1.000000 .000000 .463766 .000000 .202895 .000000 .000000 Iterasi ke-15 .000000 1.000000 .000000 .463768 .000000 .202898 .000000 .000000


1.000000 .312500 .000000 .000000

1.000000 .640625 .250000 .000000

1.000000 1.000000 1.000000 1.000000

1.000000 .609619 .277344 .000000

1.000000 .778870 .505859 .000000

1.000000 1.000000 1.000000 1.000000

1.000000 .646527 .338470 .000000

1.000000 .794691 .532238 .000000

1.000000 1.000000 1.000000 1.000000

1.000000 .651425 .346585 .000000

1.000000 .796782 .535702 .000000

1.000000 1.000000 1.000000 1.000000

1.000000 .652074 .347661 .000000

1.000000 .797059 .536162 .000000

1.000000 1.000000 1.000000 1.000000

1.000000 .652161 .347804 .000000

1.000000 .797096 .536223 .000000

1.000000 1.000000 1.000000 1.000000

1.000000 .652172 .347823 .000000

1.000000 .797101 .536231 .000000

1.000000 1.000000 1.000000 1.000000

1.000000 .652174 .347826 .000000

1.000000 .797101 .536232 .000000

1.000000 1.000000 1.000000 1.000000

Bab VI

Supardi, M.Si Iterasi ke-16 .000000 1.000000 .000000 .463768 .000000 .202898 .000000 .000000

1.000000 .652174 .347826 .000000

1.000000 .797101 .536232 .000000

1.000000 1.000000 1.000000 1.000000

Contoh eksekusi untuk penyelesaian persamaan Laplace menggunakan metode relaksasi dengan parameter relaksasi omega =1.25 Iterasi .000000 .000000 .000000 .000000 Iterasi .000000 .000000 .000000 .000000 Iterasi .000000 .000000 .000000 .000000 Iterasi .000000 .000000 .000000 .000000 Iterasi .000000 .000000 .000000 .000000 Iterasi .000000 .000000 .000000 .000000 Iterasi .000000 .000000 .000000 .000000

ke-1 1.000000 .300000 .000000 .000000 ke-3 1.000000 .468270 .167400 .000000 ke-5 1.000000 .463300 .202917 .000000 ke-7 1.000000 .463740 .202854 .000000 ke-9 1.000000 .463770 .202902 .000000 ke-10 1.000000 .463768 .202899 .000000 ke-11 1.000000 .463768 .202898 .000000

1.000000 .390000 .000000 .000000

1.000000 .807000 .300000 .000000

1.000000 1.000000 1.000000 1.000000

1.000000 .659511 .363960 .000000

1.000000 .799566 .538470 .000000

1.000000 1.000000 1.000000 1.000000

1.000000 .651879 .348104 .000000

1.000000 .796899 .536191 .000000

1.000000 1.000000 1.000000 1.000000

1.000000 .652166 .347799 .000000

1.000000 .797104 .536242 .000000

1.000000 1.000000 1.000000 1.000000

1.000000 .652174 .347827 .000000

1.000000 .797102 .536232 .000000

1.000000 1.000000 1.000000 1.000000

1.000000 .652174 .347826 .000000

1.000000 .797101 .536232 .000000

1.000000 1.000000 1.000000 1.000000

1.000000 .652174 .347826 .000000

1.000000 .797101 .536232 .000000

1.000000 1.000000 1.000000 1.000000

SOAL LATIHAN


Bab VI

Supardi, M.Si

1. Jelaskan dengan singkat perbedaan antara persamaan diferensial hiperbolik, parabolik dan eliptik serta berikan contoh masing-masing. Apakah perbedaan fisis terpenting antara persamaan hiperbolik dan parabolik di satu sisi dan persamaan eliptik di sisi lain. 2.

Apakah

persamaan-persamaan

diferensial

berikut

merupakan

persamaan

hiperbolik, parabolik atau elipti? ∂2f ∂2f ∂2f + 2 − = 0 ∂ t2 ∂ x2 ∂ t∂ x ∂f ∂ 2f ∂ 2f ∂2f 7 + 2 − −2 2 = 0 ∂t ∂t ∂ t∂ x ∂x 2 2 2 ∂ f ∂ f ∂ f ∂f 3 2 + 2 +6 + = ex ∂x ∂t ∂ t∂ x ∂ t 2 2 ∂ f ∂ f ∂ 2f ∂f ∂f + − 2 + + 2 = 0 2 2 ∂t ∂x ∂ t∂ x ∂ x ∂t ∂  ∂f ∂f  ∂  ∂f ∂f  3 +3 − f+  2 −  = sin ( x + t ) ∂t ∂t ∂x  ∂x ∂t ∂x

a.

3

b. c. d. e. 3. Persamaan

∂f ∂3f = β ∂t ∂ x3 dapat dinyatakan dalam persamaan beda hingga sebagai berikut f n(

m + 1)

= f n(

m)

(

+ µ f n(+ 2) − 3 f n(+ 1) + 3 f n( m

m

m)

)

− f n(− 1) µ = m

βδ t δ x3

tentukan kesalahan pembulatan dari persamaan beda tersebut. 4. Jelaskan dengan singkat, apa yang saudara ketahui tentang masalah nilai awal dan masalah nilai batas. Apa yang mmbedakan keduanya, dan berikan contohnya masing-masing. 5. Skema Lax ditulis sebagai u in + 1 =

(

)

(

1 n 7 n u i + 1 + u in− 1 − u i + 1 + u in− 1 2 2

)

dengan γ ≡ α ∆ t / ∆ x . Tunjukkan bahwa skema tersebut stabil jika 0 < γ < 1 . 6. Persamaan diferensial parsial diberikan oleh


Bab VI

Supardi, M.Si ∂u ∂u + a ( x, t ) = s ( x, t ) ∂t ∂x

dengan a( x, t ) = 3 x + 0.1

s( x, t ) = 1 − x 2 + 0.1t

Dengan mengasumsikan syarat awal diberikan oleh u ( x, t ) = 1 untuk t = 0 , tentukan penyelesaiaan untuk masalah tersebut. 7. Persamaan difusi dalam suatu ruang, dimana konstanta difusi D berubah terhadap ruang D = D ( x ) dinyatakan oleh ∂f ∂  ∂f  = D  ∂t ∂x ∂x

atau

∂f ∂ 2f ∂D ∂f = D 2 + ∂t ∂x ∂x ∂x

Tunjukkan bahwa persamaan beda dinyatakan sebagai f nm + 1 − f nm f nm+ 1 − 2 f nm + f nm− 1  Dn + 1 − Dn − 1   f nm+ 1 − f nm− 1  = Dn +   δt δ x2 2δ x    2δ x  tidak konservatif, yaitu bahwa

∫

fdx tidak kekal.

8. Buatlah suatu skema alternatif yang menunjukkan bahwa persamaan difusi yang dinyatakan pada soal nomor 1 tersebut konservatif. 9. Tunjukkan bahwa skema Lax untuk penyelesaiam persamaan adveksi ekivalen dengan ∂f ∂ f  δ x2 1 2  ∂ 2 f = −u +  − u δ t  2 + suku orde lebih tinggi ∂t ∂ x  2δ t 2  ∂x 10. Ujilah perilaku penyelesaian like- gelombang f = exp ( i ( kx − ω t ) ) dalam skema Lax dan jelaskan perilaku dalam suku-suku difusi.


BAB VI PERSAMAAN DIFERENSIAL PARSIAL (PDP)

Recommend Documents