PENYELESAIAN PERSAMAAN GELOMBANG AIR DANGKAL DENGAN BEBERAPA METODE NUMERIS

PLAGIAT MERUPAKAN TINDAKAN TIDAK TERPUJI

PENYELESAIAN PERSAMAAN GELOMBANG AIR DANGKAL DENGAN BEBERAPA METODE NUMERIS

SKRIPSI

Diajukan untuk Memenuhi Salah Satu Syarat Memperoleh Gelar Sarjana Matematika Program Studi Matematika

Disusun oleh: Ilga Purnama Sari NIM: 123114023

PROGRAM STUDI MATEMATIKA JURUSAN MATEMATIKA FAKULTAS SAINS DAN TEKNOLOGI UNIVERSITAS SANATA DHARMA YOGYAKARTA 2016


PENYELESAIAN PERSAMAAN GELOMBANG AIR DANGKAL DENGAN BEBERAPA METODE NUMERIS

SKRIPSI

Diajukan untuk Memenuhi Salah Satu Syarat Memperoleh Gelar Sarjana Matematika Program Studi Matematika

Disusun oleh: Ilga Purnama Sari NIM: 123114023

PROGRAM STUDI MATEMATIKA JURUSAN MATEMATIKA FAKULTAS SAINS DAN TEKNOLOGI UNIVERSITAS SANATA DHARMA YOGYAKARTA 2016

i


SOLUTION TO THE SHALLOW WATER WAVE EQUATIONS WITH SOME NUMERICAL METHODS

A THESIS

Presented as Partial Fulfillment of the Requirements to Obtain the Degree of Sarjana Matematika Mathematics Study Program

Written by: Ilga Purnama Sari Student ID: 123114023

MATHEMATICS STUDY PROGRAM DEPARTMENT OF MATHEMATICS FACULTY OF SCIENCE AND TECHNOLOGY SANATA DHARMA UNIVERSITY YOGYAKARTA 2016

ii


iii


iv


v


vi


ABSTRAK

Persamaan gelombang air dangkal adalah persamaan yang memodelkan aliran air di tempat terbuka. Persamaan gelombang air dangkal seringkali disebut sistem Saint-Venant yang diturunkan dari hukum kekekalan massa dan momentum. Dalam skripsi ini, persamaan gelombang air dangkal diselesaikan menggunakan beberapa metode, yaitu metode volume hingga, metode beda hingga grid kolokasi dan metode beda hingga grid selang-seling. Metode volume hingga bekerja dengan cara membagi domain ruang menjadi sebanyak berhingga sel, kemudian dihitung rata-rata kuantitas untuk masing-masing sel. Metode beda hingga grid kolokasi dikerjakan secara implisit yaitu membagi domain ruang menjadi sebanyak berhingga titik, kemudian persamaan gelombang air dangkal didiskretkan dan membentuk sebuah sistem persamaan linear. Terakhir, metode beda hingga grid selang-seling bekerja dengan cara membagi domain perhitungan ruang secara selang-seling. Kedalaman air dihitung pada grid dengan indeks bilangan bulat dan kecepatan air dihitung pada grid dengan indeks pecahan. Penggunaan metode yang tepat akan menghasilkan solusi yang akurat untuk persamaan gelombang air dangkal. Penelitian ini menguji beberapa metode numeris yang bisa digunakan untuk menyelesaikan persamaan gelombang air dangkal satu dimensi. Pengujian dilakukan menggunakan simulasi numeris. Analisis hasil simulasi dilakukan dengan observasi hasil simulasi di setiap metode dan membandingkannya dengan hasil solusi eksak.

vii


ABSTRACT The shallow water wave equations model water flows in an open channel. The shallow water wave equations are often called the Saint-Venant system derived from the conservations of mass and momentum. In this thesis, the shallow water wave equations are solved using several methods, the Lax-Friedrichs finite volume method, collocation grid finite difference method and staggered grid finite difference method. The finite volume method works by dividing the spatial domain into a finite number of cells, then calculating the average quantity for each cell. The collocation grid finite difference method divides the spatial domain into a finite number of computational points for the shallow water equation discretization and forms a linear system of equations. Finally, the staggered grid finite difference method works by discretising the computational domain into staggered spatial partitions. The staggered finite diference means that we approximate the quantities of interest of the shallow water equations on different cells. In the staggered formulation, water depth is calculated at full grid points and water velocity is calculated at half grid points. The appropriate method will produce an accurate solution for the shallow water wave equations. This study examines several numerical methods for solving the one dimensional shallow water wave equations. We investigate the performance of these numerical methods using numerical simulations. Analysis of simulation results are done by observing the results of each method and comparing the numerical results with the exact solution.

viii


DAFTAR ISI HALAMAN JUDUL.............................................................................................. i HALAMAN JUDUL DALAM BAHASA INGGRIS .......................................... ii HALAMAN PERSETUJUAN PEMBIMBING .................................................. iii HALAMAN PENGESAHAN .............................................................................. iv HALAMAN PERNYATAAN KEASLIAN KARYA .......................................... v LEMBAR PERNYATAAN PERSETUJUAN PUBLIKASI............................... vi ABSTRAK .......................................................................................................... vii ABSTRAK DALAM BAHASA INGGRIS ....................................................... viii DAFTAR ISI ........................................................................................................ ix BAB I PENDAHULUAN .................................................................................... 1 A. Latar Belakang ......................................................................................... 1 B. Rumusan Masalah .................................................................................... 4 C. Pembatasan Masalah ................................................................................ 5 D. Tujuan Penulisan ...................................................................................... 5 E. Metode Penulisan ..................................................................................... 5 F. Manfaat Penulisan .................................................................................... 6 G. Sistematika Penulisan .............................................................................. 6 BAB II PERSAMAAN DIFERENSIAL ........................................................... 8

ix


A. Integral ...................................................................................................... 8 B. Klasifikasi Persamaan Diferensial .......................................................... 11 C. Nilai Eigen dan Vektor Eigen ................................................................. 15 D. Persamaan Diferensial Hiperbolik .......................................................... 16 E. Penurunan Numeris ................................................................................. 18 F. Karakteristik Persamaan Gelombang Air Dangkal ................................. 25 BAB III METODE NUMERIS UNTUK PERSAMAAN GELOMBANG AIR DANGKAL ........................................................................................................ 28 A. Solusi Eksak Persamaan Gelombang Air Dangkal ................................. 28 B. Penurunan Persamaan Gelombang Air Dangkal Satu Dimensi .............. 29 C. Masalah Bendungan Bobol ..................................................................... 35 D. Metode Volume Hingga Lax-Friedrichs ................................................. 37 E. Metode Beda Hingga Grid Kolokasi ....................................................... 44 F. Metode Beda Hingga Grid Selang-Seling ............................................... 53 BAB IV PERBANDINGAN BEBERAPA METODE NUMERIS DALAM PENYELESAIAN MODEL GELOMBANG AIR DANGKAL .................... 57 A. Metode Volume Hingga Lax-Friedrichs ................................................. 58 B. Metode Beda Hingga Grid Kolokasi ....................................................... 62 C. Metode Beda Hingga Grid Selang-Seling ............................................... 66

BAB V PENUTUP ............................................................................................. 71

x


A. Kesimpulan ............................................................................................ 71 B. Saran ........................................................................................................ 72 DAFTAR PUSTAKA ........................................................................................ 73 LAMPIRAN ....................................................................................................... 74

xi


BAB I PENDAHULUAN

Dalam bab ini akan dijelaskan latar belakang, rumusan dan pembatasan masalah, tujuan dan manfaat penulisan, juga akan disertakan sistematika penulisan. A. Latar Belakang Persamaan diferensial merupakan persamaan yang menghubungkan suatu fungsi beserta turunan-turunannya. Banyak masalah fisis yang dapat diselesaikan dengan persamaan diferensial. Masalah fisis merupakan masalah yang berhubungan dengan hukum alam yang dibahas dalam ilmu fisika. Masalah fisis seperti fluida dapat diselesaikan dengan menggunakan persamaan diferensial parsial. Fluida merupakan zat yang dapat mengalir. Zat itu dapat berupa gas atau cairan. Aliran fluida merupakan salah satu masalah fisis yang sering dijumpai dalam kehidupan sehari-hari. Masalah yang sudah pernah ada adalah terjadinya bencana alam seperti bobolnya bendungan air atau tsunami. Bencana alam tersebut disebabkan oleh aliran air dalam skala besar. Aliran tersebut dapat dimodelkan secara matematis (Crowhurst, 2013; Crowhurst dan Li, 2013). Model gelombang air yang sudah ada salah satunya adalah model gelombang air dangkal atau Shallow Water Wave Equations. Dangkal dalam arti matematis adalah amplitudo gelombang jauh lebih kecil dibandingkan panjang gelombangnya.

1

PLAGIAT MERUPAKAN TINDAKAN TIDAK TERPUJI 2

Dengan demikian, nilai perbandingan antara amplitudo gelombang (𝑎) dengan panjang gelombang (𝜆) ditulis

𝑎 𝜆

≪ 1.

Penyelesaian persamaan gelombang air dangkal memiliki dua komponen penting yang tidak diketahui yaitu kedalaman dan kecepatan air, dengan ℎ(𝑥, 𝑡) adalah kedalaman air dan 𝑢(𝑥, 𝑡) adalah kecepatan air. Di sini, 𝑡 adalah variabel yang menyatakan waktu dan 𝑥 adalah variabel yang menyatakan ruang satu dimensi. Persamaan gelombang air dangkal dalam bentuk sistem persamaan diferensial parsial dinyatakan oleh dua persamaan simultan, yaitu ℎ𝑡 + (𝑢ℎ)𝑥 = 0

(1.1)

1 (𝑢ℎ)𝑡 + (𝑢2 ℎ + 𝑔ℎ2 )𝑥 = −𝑔ℎ𝑧𝑥 2

(1.2)

dan

dengan 𝑧(𝑥) adalah ketinggian tanah, dan 𝑔 adalah konstanta percepatan gravitasi. Ilustrasi gelombang air dangkal ditunjukkan dalam Gambar 1.1. 𝜆

y

𝑤(𝑥, 𝑡)

x

Gambar 1.1: Gelombang air dangkal dan variabel-variabel terkait dalam persamaan gelombang air dangkal.


Pada skripsi ini akan diselesaikan persamaan gelombang air dangkal terkait dengan masalah bendungan bobol. Masalah bendungan bobol memiliki beberapa asumsi, syarat awal dan syarat batas yang akan dibahas lebih lanjut pada Bab III. Ilustrasi bendungan air dapat dilihat pada Gambar 1.2. Permukaan air

ℎ1

Bendungan air

Permukaan air ℎ0 𝑥

Gambar 1.2: Bendungan air Secara umum, solusi persamaan gelombang air dangkal tersebut cukup sulit untuk dicari secara analitis, sehingga diperlukan cara lain untuk memecahkannya. Metode numeris adalah salah satu cara untuk memperoleh solusi persamaan gelombang air dangkal tersebut. Banyak metode numeris yang telah dikembangkan sebelumnya untuk memecahkan solusi persamaan tersebut, mulai dari metode karakteristik, metode beda hingga, metode elemen hingga, metode volume hingga dan sebagainya. Metode beda hingga dikembangkan berdasarkan diskritisasi langsung dari persamaan diferensial yang dipandang (LeVeque, 1992). Metode beda hingga unggul dalam kemudahan komputasi. Dalam tugas akhir ini akan dibahas penyelesaian persamaan gelombang air dangkal dengan metode beda hingga dan metode volume hingga, karena perumusan kedua metode tersebut sederhana.


Metode beda hingga terbagi atas dua model yaitu model grid kolokasi dan model grid selang-seling. Pada grid kolokasi ditentukan nilai pendekatan untuk semua variabel ℎ dan 𝑢 yang tidak diketahui secara bersamaan. Pada grid selangseling ditentukan pendekatan variabel ℎ dan 𝑢 secara selang-seling. Salah satu referensi tentang grid kolokasi adalah LeVeque (1992). Salah satu referensi tentang grid selang-seling adalah Stelling dan Duinmeijer (2003). Penelitian ini akan membandingkan hasil perhitungan terbaik antara metode volume hingga, metode beda hingga grid kolokasi dan metode beda hingga grid selang-seling. Dari ketiga metode tersebut diperoleh hasil perhitungan terbaik yang dapat memperbaiki metode beda hingga dengan tidak ada getaran semu (artificial oscillation) pada hasil simulasi aliran air. Fokus penelitian ini adalah mengembangkan metode numeris dengan metode beda hingga dan volume hingga untuk menyimulasikan aliran air.

B. Rumusan Masalah Permasalahan yang dirumuskan dalam skripsi ini ada empat, yaitu: 1. Bagaimana memodelkan gelombang air dangkal? 2. Bagaimana menyelesaikan dan menyimulasikan persamaan gelombang air dangkal dengan menggunakan metode volume hingga? 3. Bagaimana menyelesaikan dan menyimulasikan persamaan gelombang air dangkal dengan menggunakan metode beda hingga grid kolokasi? 4. Bagaimana menyelesaikan dan menyimulasikan persamaan gelombang air dangkal dengan menggunakan metode beda hingga grid selang-seling?


C. Pembatasan Masalah Pembahasan masalah dalam skripsi ini akan dibatasi pada memodelkan gelombang air dangkal satu dimensi dan mencari penyelesaian persamaan gelombang air dangkal dengan metode beda hingga dan volume hingga.

D. Tujuan Penulisan Skripsi ini mempunyai dua tujuan, yaitu: 1. Merumuskan dan menyelesaikan persamaan gelombang air dangkal dengan beberapa metode numeris, yaitu metode volume hingga Lax-Friedrichs, metode beda hingga grid kolokasi dan metode beda hingga grid selang-seling. 2. Membandingkan beberapa hasil simulasi numeris menggunakan metode volume hingga Lax-Friedrichs, metode beda hingga grid kolokasi dan metode beda hingga grid selang-seling. Dari hasil simulasi tersebut kemudian dipilih hasil yang terbaik yang memuat galat yang paling kecil. 3. Menyimulasikan metode volume hingga Lax-Friedrichs, metode beda hingga grid kolokasi dan metode beda hingga grid selang-seling untuk masalah bendungan bobol.

E. Metode Penulisan Metode penulisan yang digunakan adalah studi pustaka dari buku-buku dan jurnal serta praktek simulasi numeris.


F. Manfaat Penulisan Dengan memodelkan persamaan gelombang air, kita dapat 1. Menyimulasikan terjadinya banjir. 2. Memperkirakan daerah mana saja yang akan tenggelam. 3. Memperkirakan kecepatan dan kedalaman air pada daerah tersebut.

G. Sistematika Penulisan Sistematika penulisan tugas akhir ini adalah sebagai berikut:

BAB I

PENDAHULUAN

A. Latar Belakang Masalah B. Rumusan Masalah C. Batasan Masalah D. Metode Penulisan E. Tujuan Penulisan F. Manfaat Penulisan BAB II

PERSAMAAN DIFERENSIAL

A. Integral B. Klasifikasi Persamaan Diferensial C. Nilai eigen dan vektor eigen D. Persamaan Diferensial Hiperbolik E. Penurunan Numeris F. Karakteristik Persamaan Gelombang Air Dangkal BAB III METODE NUMERIS UNTUK PERSAMAAN GELOMBANG AIR DANGKAL A. Solusi Eksak Persamaan Gelombang Air Dangkal B. Penurunan Persamaan Gelombang Air Dangkal Satu Dimensi C. Masalah Bendungan Bobol D. Metode Volume Hingga Lax-Friedrichs


E. Metode Beda Hingga Grid Kolokasi F. Metode Beda Hingga Grid Selang-Seling BAB IV PERBANDINGAN BEBERAPA METODE NUMERIS DALAM PENYELESAIAN MODEL GELOMBANG AIR DANGKAL A. Metode Volume Hingga Lax-Friedrichs B. Metode Beda Hingga Grid Kolokasi C. Metode Beda Hingga Grid Selang-Seling BAB V

PENUTUP

A. Kesimpulan B. Saran DAFTAR PUSTAKA LAMPIRAN


BAB II PERSAMAAN DIFERENSIAL

Landasan teori skripsi ditulis dalam bab ini. Landasan teori tersebut meliputi: integral, klasifikasi persamaan diferensial, nilai eigen dan vektor eigen, persamaan diferensial hiperbolik, penurunan numeris, dan karakteristik persamaan gelombang air dangkal.

A. Integral Pada bagian ini dibahas mengenai integral yang meliputi definisi dan contoh dari integral tentu dan aturan Leibniz. Definisi 2.1 Jika diberikan suatu fungsi 𝑓(𝑥) pada suatu interval 𝐼 dan berlaku 𝐹 ′ (𝑥) = 𝑓(𝑥), untuk suatu 𝐹(𝑥), maka 𝐹(𝑥) adalah anti turunan dari 𝑓(𝑥). Dengan kata lain 𝐹 ′ (𝑥) = 𝑓(𝑥). Contoh 2.1 Carilah suatu anti turunan dari 𝑓(𝑥) = 2𝑥 2 pada (−∞, ∞). Penyelesaian: Fungsi 𝐹(𝑥) = 2𝑥 3 bukan anti turunannya karena turunan 2𝑥 3 adalah 6𝑥 2 . Tetapi 2

2

hal ini menyarankan 𝐹(𝑥) = 3 𝑥 3 , yang memenuhi 𝐹 ′ (𝑥) = 3 3𝑥 2 = 2𝑥 2 . Dengan 2

demikian, suatu anti turunan dari 𝑓 adalah 3 𝑥 3 .

8


Anti turunan dinotasikan dengan ∫ … 𝑑𝑥. Notasi tersebut menunjukkan anti turunan terhadap 𝑥. Anti turunan biasanya disebut integral tak tentu. Integral Tentu Perhatikan Gambar 2.1 berikut ini.

𝑦

𝑦 = 𝑓(𝑥)

𝑥

𝑎

𝑏

Gambar 2.1: Ilustrasi fungsi satu variabel Untuk menghitung luasan dibawah kurva 𝑦 = 𝑓(𝑥) pada interval [𝑎, 𝑏], dapat dihitung dengan cara aproksimasi yaitu dengan membagi interval [𝑎, 𝑏] menjadi 𝑛 subinterval. Subinterval tersebut memiliki panjang yang sama yaitu

𝑏−𝑎 𝑛

untuk 𝑛 >

0. Setelah membagi interval menjadi 𝑛 subinterval kemudian menghitung total jumlah luasan dari masing-masing persegi panjang yang dibentuk oleh masingmasing subinterval tersebut. Hal ini diperoleh dengan memilih 𝑥0 , 𝑥1 , … , 𝑥𝑛 dengan 𝑎 = 𝑥0 , 𝑏 = 𝑥𝑛 , dan 𝑥𝑖 − 𝑥𝑖−1 =

𝑏−𝑎 𝑛


untuk 𝑖 = 1,2, … , 𝑛. Andaikan panjang masing-masing subinterval yaitu

𝑏−𝑎 𝑛

dinotasikan dengan ∆𝑥 = 𝑥𝑖 − 𝑥𝑖−1 . Luas daerah dibawah kurva diaproksimasikan dengan total luas daerah yang dibentuk oleh masing-masing subinterval, aproksimasi luas di bawah kurva adalah 𝐴1 + 𝐴2 + ⋯ + 𝐴𝑛 . Artinya total luas tersebut yang dapat ditulis 𝑛

𝑓(𝑢1 )∆𝑥 + 𝑓(𝑢2 )∆𝑥 + ⋯ + 𝑓(𝑢𝑛 )∆𝑥 = ∑ 𝑓(𝑢𝑖 )∆𝑥 𝑖=1

yang disebut jumlahan Riemann fungsi 𝑓 pada interval [𝑎, 𝑏], sebagai pendekatan luas daerah di bawah kurva 𝑦 = 𝑓(𝑥) dan diatas sumbu 𝑥. Di sini, 𝑢𝑖 ∈ [𝑥𝑖−1 , 𝑥𝑖 ]. Semakin banyak subinterval yang digunakan, artinya ∆𝑥 → 0 maka semakin baik pula aproksimasi luasan tersebut dan semakin dekat dengan luasan yang sebenarnya. Dengan demikian, Luas daerah = lim ∑ 𝑓(𝑢𝑖 )∆𝑥. ∆𝑥→0

𝑖

Definisi 2.2 Andaikan 𝑓 fungsi yang terdefinisi pada [𝑎, 𝑏]. Integral tentu 𝑓 dari 𝑎 sampai 𝑏

𝑏 dinotasikan ∫𝑎 𝑓 (𝑥)𝑑𝑥, adalah 𝑏

∫ 𝑓(𝑥)𝑑𝑥 = lim ∑ 𝑓(𝑢𝑖 )∆𝑥. 𝑎

∆𝑥→0

𝑖


Aturan Leibniz Teorema 2.1 Aturan Leibniz untuk satu variabel: Jika 𝑓 adalah fungsi kontinu pada interval [𝑎, 𝑏] dan jika 𝑢(𝑥) dan 𝑣(𝑥) adalah fungsi yang dapat diturunkan terhadap 𝑥 yang nilainya terletak di interval [𝑎, 𝑏], maka 𝑑 𝑣(𝑥) 𝑑𝑣 𝑑𝑢 ∫ 𝑓(𝑡)𝑑𝑡 = 𝑓(𝑣(𝑥)) − 𝑓(𝑢(𝑥)) 𝑑𝑥 𝑢(𝑥) 𝑑𝑥 𝑑𝑥 Teorema 2.2 Aturan Leibniz untuk dua variabel: Jika 𝑓(𝑥, 𝑡) adalah fungsi sedemikian sehingga turunan parsial dari 𝑓 terhadap 𝑡 ada dan kontinu, maka 𝑣(𝑥) 𝑑 𝑣(𝑥) 𝜕𝑓 𝑑 𝑑 ∫ 𝑓(𝑥, 𝑡)𝑑𝑡 = ∫ 𝑑𝑡 + 𝑓(𝑣(𝑥), 𝑥) 𝑣(𝑥) − 𝑓(𝑢(𝑥), 𝑥) 𝑢(𝑥). 𝑑𝑥 𝑢(𝑥) 𝑑𝑥 𝑑𝑥 𝑢(𝑥) 𝜕𝑥

Bukti dapat dilihat pada buku karangan David. B dan George. C yang berjudul Basic Partial Differential Equations.

B. Klasifikasi Persamaan Diferensial Berikut ini dibahas tentang klasifikasi persamaan diferensial. Klasifikasi persamaan diferensial yang dibahas meliputi definisi dan contoh persamaan


diferensial, persamaan diferensial biasa, persamaan diferensial parsial, orde persamaan diferensial dan kelinearan suatu persamaan diferensial.

Definisi 2.3 Persamaan diferensial adalah persamaan yang melibatkan variabel-variabel tak bebas dan turunan-turunannya terhadap variabel-variabel bebas. Contoh 2.2 Persamaan di bawah ini merupakan contoh persamaan diferensial: 𝑑2𝑦 𝑑𝑦 2 + 𝑥𝑦 ( ) =0 𝑑𝑥 2 𝑑𝑥

(2.1)

𝑑4 𝑥 𝑑2𝑥 + 5 + 3𝑥 = sin 𝑡 𝑑𝑡 4 𝑑𝑡 2

(2.2)

𝜕𝑣 𝜕𝑣 + =𝑣 𝜕𝑠 𝜕𝑡

(2.3)

𝜕 2𝑢 𝜕 2𝑢 𝜕 2𝑢 + + = 0. 𝜕𝑥 2 𝜕𝑦 2 𝜕𝑧 2

(2.4)

Definisi 2.4 Persamaan diferensial biasa adalah persamaan diferensial yang melibatkan turunan biasa beserta satu atau lebih variabel tak bebas terhadap satu variabel bebas. Contoh 2.3 Persamaan (2.1) dan (2.2) adalah persamaan diferensial biasa. Pada persamaan (2.1) variabel 𝑥 adalah suatu variabel bebas, dan variabel 𝑦 adalah variabel tak bebas. Pada persamaan (2.2), variabel 𝑡 adalah variabel bebas, dengan 𝑥 adalah variabel tak bebasnya.


Definisi 2.5 Persamaan diferensial parsial adalah persamaan diferensial yang melibatkan turunan parsial dari satu atau lebih variabel tak bebas terhadap lebih dari satu variabel bebas. Contoh 2.4 Persamaan (2.3) dan (2.4) adalah persamaan diferensial parsial. Pada persamaan (2.3), variabel 𝑠 dan 𝑡 adalah variabel bebas dan 𝑣 adalah variabel tak bebasnya. Pada persaman (2.4) terdapat tiga variabel bebas yaitu 𝑥, 𝑦, dan 𝑧. Pada persamaan (2.4) variabel tak bebasnya adalah 𝑢.

Definisi 2.6 Orde dari persamaan diferensial adalah tingkat tertinggi dari turunan yang terkandung dalam persamaan diferensial. Contoh 2.5 Persamaan diferensial biasa (2.1) adalah persamaan diferensial orde kedua, karena tingkat tertinggi dari turunan pada persamaan tersebut adalah dua. Persamaan (2.2) adalah persamaan diferensial biasa orde keempat. Persamaan (2.3) termasuk persamaan diferensial parsial orde pertama. Persamaan (2.4) merupakan persamaan diferensial parsial orde kedua.

Definisi 2.7 Suatu persamaan diferensial biasa orde ke-𝑛 𝐹(𝑥, 𝑦, 𝑦 ′ , 𝑦 ′′ , … , 𝑦 (𝑛) ) = 0


dikatakan linear jika 𝐹 merupakan suatu fungsi linear dari variabel 𝑦, 𝑦 ′ , 𝑦 ′′ , … , 𝑦 (𝑛) ; definisi yang sama juga berlaku untuk persamaan diferensial parsial. Secara umum persamaan diferensial biasa linear orde 𝑛 dituliskan sebagai 𝑎0 (𝑥)𝑦 (𝑛) + 𝑎1 (𝑥)𝑦 (𝑛−1) + ⋯ + 𝑎𝑛 (𝑥)𝑦 = 𝑏(𝑥)

(2.5)

dengan 𝑎0 tidak sama dengan nol. Contoh 2.6 Persamaan diferensial biasa berikut keduanya linear. Pada kedua persamaan berikut, variabel 𝑦 adalah variabel tak bebas. Perhatikan bahwa 𝑦 dan turunanturunannya terjadi dengan pangkat satu saja dan tidak ada perkalian dari 𝑦 dan/ atau turunan dari 𝑦. 𝑑2𝑦 𝑑𝑦 + 5 + 6𝑦 = 0 𝑑𝑥 2 𝑑𝑥

(2.6)

𝑑4𝑦 𝑑3𝑦 𝑑𝑦 2 + 𝑥 + 𝑥3 = 𝑥𝑒 𝑥 . 4 3 𝑑𝑥 𝑑𝑥 𝑑𝑥

(2.7)

Definisi 2.8 Suatu persamaan diferensial biasa yang tidak memiliki bentuk (2.5) dinamakan persamaan diferensial biasa tak linear. Contoh 2.7 Persamaan diferensial biasa berikut semuanya tak linear: 𝑑2𝑦 𝑑𝑦 +5 + 6𝑦 2 = 0 2 𝑑𝑥 𝑑𝑥

(2.8)


𝑑2 𝑦 𝑑𝑦 3 + 5 ( ) + 6𝑦 = 0 𝑑𝑥 2 𝑑𝑥

(2.9)

𝑑2𝑦 𝑑𝑦 + 5𝑦 + 6𝑦 = 0 2 𝑑𝑥 𝑑𝑥

(2.10)

Persamaan (2.8) tak linear karena variabel tak bebas 𝑦 terdapat pada pangkat kedua dalam bentuk 6𝑦 2. Persamaan (2.9) juga tak linear karena terdapat bentuk 𝑑𝑦 3

5 (𝑑𝑥 ) yang melibatkan pangkat tiga pada turunan pertama. Persamaan (2.10) tak 𝑑𝑦

linear karena pada bentuk 5𝑦 𝑑𝑥 melibatkan perkalian terhadap variabel tak bebas dan turunan pertamanya.

C. Nilai Eigen dan Vektor Eigen Berikut dibahas mengenai definisi dan contoh dari nilai eigen dan vektor eigen. Definisi 2.9 Jika 𝐴 adalah matriks 𝑛 × 𝑛, maka vektor tak nol 𝐱 di ℝ𝑛 disebut vektor eigen dari 𝐴 jika 𝐴𝒙 merupakan perkalian skalar dengan 𝐱 atau dapat ditulis 𝐴𝐱 = 𝜆𝐱 untuk suatu skalar 𝜆. Skalar 𝜆 disebut nilai eigen dari 𝐴 dan 𝐱 disebut vektor eigen yang bersesuaian dengan 𝝀. Contoh 2.8 1 Vektor 𝐱 = [ ] adalah vektor eigen dari 2


3 0 𝐴=[ ] 8 −1 yang bersesuaian dengan nilai eigen 𝜆 = 3, karena 𝐴𝐱 = [

3 0 1 3 ] [ ] = [ ] = 3𝐱. 8 −1 2 6

Secara geometri, perkalian matriks 𝐴 dengan vektor 𝐱 memiliki kelipatan 3 terhadap vektor 𝐱. Ilustrasi secara geometri ditunjukkan dalam Gambar 2.2.

𝑦 3𝑥

6

2

𝑥

1

3

𝑥

Gambar 2.2: Ilustrasi geometri vektor eigen. D. Persamaan Diferensial Hiperbolik Persamaan diferensial hiperbolik dapat digunakan untuk memodelkan banyak fenomena yang melibatkan pergerakan gelombang. Perhatikan bentuk persamaan diferensial berikut 𝑞𝑡 (𝑥, 𝑡) + 𝐴𝑞𝑥 (𝑥, 𝑡) = 0, Di sini 𝑞𝑡 =

𝜕𝑞 𝜕𝑡

𝜕𝑞

dan 𝑞𝑥 = 𝜕𝑥 .

(2.11)


Dalam kasus yang paling sederhana yaitu koefisien konstan dan linear. Dalam hal ini 𝑞 ∶ ℝ × ℝ → ℝ𝑛 adalah vektor dengan 𝑛 komponen yang menyatakan fungsi yang tidak diketahui (tekanan, kecepatan, dan sebagainya) yang ingin ditentukan, dan 𝐴 suatu konstan merupakan matriks real berukuran 𝑛 × n. Andaikan 𝐴 = 𝑢̅ suatu konstan yang menyatakan kecepatan perambatan (aliran pada pipa satu dimensi misalnya), maka persamaan (2.11) menjadi 𝑞𝑡 (𝑥, 𝑡) + 𝑢̅𝑞𝑥 (𝑥, 𝑡) = 0,

(2.12)

persamaan ini disebut persamaan adveksi. Persamaan 𝑞𝑡 (𝑥, 𝑡) + 𝐴𝑞𝑥 (𝑥, 𝑡) = 0 adalah persamaan hiperbolik, jika matriks 𝐴 memiliki nilai eigen real dan berkorespondensi dengan 𝑛 vektor eigen yang bebas linear. Artinya, semua vektor dalam ℝ𝑛 dapat secara tunggal diuraikan sebagai kombinasi linear dari nilai-nilai eigen tersebut. Secara formal definisi persamaan diferensial hiperbolik sebagai berikut.

Definisi 2.10 Suatu sistem linear dengan bentuk 𝑞𝑡 + 𝐴𝑞𝑥 = 0 dikatakan hiperbolik jika matriks 𝐴 yang berukuran 𝑛 × n dapat didiagonalkan dengan nilai eigen real. Secara khusus untuk persamaan adveksi, diketahui bahwa 𝐴 = 𝑢̅, yang merupakan suatu konstanta real. Jadi 𝐴 dapat didiagonalkan oleh nilai 𝐴 itu sendiri dan nilai eigen dari 𝐴 adalah 𝐴 itu sendiri. Dengan demikian, persamaan adveksi merupakan persamaan diferensial hiperbolik. Keterangan lengkap tentang


persamaan diferensial hiperbolik dapat ditemukan dalam buku karangan LeVeque (2004).

E. Penurunan Numeris Pada subbab ini dibahas mengenai penurunan numeris beserta contohnya dan penjelasan mengenai tiga hampiran dalam menghitung turunan numerik yaitu hampiran beda maju, hampiran beda mundur dan hampiran beda pusat. Definisi 2.11 Suatu turunan fungsi didefinisikan dengan 𝑓(𝑥 + ∆𝑥) − 𝑓(𝑥) . ∆𝑥→0 ∆𝑥

𝑓 ′ (𝑥) = lim

Seringkali fungsi 𝑓(𝑥) tidak diketahui secara eksplisit, tetapi hanya diketahui beberapa titik data saja. Seringkali 𝑓(𝑥) diketahui secara eksplisit tetapi karena bentuknya yang sangat rumit sehingga untuk menentukan fungsi turunannya juga sulit, misalnya pada fungsi-fungsi berikut ini:

(a). 𝑓(𝑥) =

√cos(2𝑥 2 ) + 𝑥 tan(3𝑥) , 2𝑥 𝑥 sin(𝑥) + 𝑒 − cos(𝑥)

(b). 𝑓(𝑥) = 𝑥𝑒 (2𝑥+2) ln(4𝑥 2 ). Perhitungan nilai turunan pada fungsi (a) dan (b) dapat dikerjakan secara numerik. Nilai turunan yang diperoleh merupakan nilai hampiran dan diharapkan nilai galatnya sekecil mungkin.


Tiga Hampiran dalam Menghitung Turunan Numerik Turunan adalah limit dari hasil bagi pengurangan dua buah nilai yang besar 𝑓(𝑥 + ∆𝑥) − 𝑓(𝑥) dan membaginya dengan bilangan yang kecil (∆𝑥). Misal diberikan nilai-nilai 𝑥 di 𝑥0 − ∆𝑥, 𝑥0 , dan 𝑥0 + ∆𝑥, serta nilai fungsi untuk nilainilai 𝑥 tersebut. Titik-titik yang diperoleh adalah (𝑥−1 , 𝑓−1 ), (𝑥0 , 𝑓0 ), dan (𝑥1 , 𝑓1 ), yang dalam hal ini 𝑥−1 = 𝑥0 − ∆𝑥 dan 𝑥1 = 𝑥0 + ∆𝑥. Terdapat tiga hampiran dalam menghitung nilai 𝑓 ′ (𝑥0 ): 1. Hampiran Beda Maju Diketahui fungsi 𝑦 = 𝑓(𝑥0 ). Akan ditunjukkan 𝑓 ′ (𝑥0 ) dengan hampiran beda maju. 𝑓(𝑥0 + ∆𝑥) − 𝑓(𝑥0 ) ∆𝑥→0 ∆𝑥

𝑓 ′ (𝑥0 ) = lim ≈

𝑓(𝑥0 + ∆𝑥) − 𝑓(𝑥0 ) ∆𝑥 =

𝑓1 − 𝑓0 . ∆𝑥

2. Hampiran Beda Mundur Diketahui fungsi 𝑦 = 𝑓(𝑥0 ). Akan ditunjukkan 𝑓 ′ (𝑥0 ) dengan hampiran beda mundur. 𝑓(𝑥0 ) − 𝑓(𝑥0 − ∆𝑥) ∆𝑥→0 ∆𝑥

𝑓 ′ (𝑥0 ) = lim ≈

𝑓(𝑥0 ) − 𝑓(𝑥0 − ∆𝑥) ∆𝑥 =

𝑓0 − 𝑓1 . ∆𝑥


3. Hampiran Beda Pusat Diketahui fungsi 𝑦 = 𝑓(𝑥0 ). Akan ditunjukkan 𝑓 ′ (𝑥0 ) dengan hampiran beda pusat.

𝑓(𝑥0 + ∆𝑥) − 𝑓(𝑥0 − ∆𝑥) ∆𝑥→0 2∆𝑥

𝑓 ′ (𝑥0 ) = lim ≈

𝑓(𝑥0 + ∆𝑥) − 𝑓(𝑥0 − ∆𝑥) 2∆𝑥 =

𝑓1 − 𝑓−1 . 2∆𝑥

Tafsiran geometri dari ketiga pendekatan di atas diperlihatkan pada Gambar 2.3. 𝑦

𝑦

𝑦1 𝑦0

𝑦0

𝑦 = 𝑓(𝑥)

𝑦 = 𝑓(𝑥) 𝑦−1

∆𝑥 𝑥−1

𝑥0 (𝐚)

𝑥1

∆𝑥

𝑥

𝑥−1

𝑥0

(𝐛)

𝑥1

𝑥

𝑦

𝑦1

𝑦 = 𝑓(𝑥)

𝑦−1

(𝐜)2∆𝑥 𝑥−1

𝑥0

𝑥1

𝑥

Gambar 2.3: Tiga pendekatan dalam perhitungan numeris; (a) Hampiran beda maju, (b) Hampiran beda mundur, dan (c) Hampiran beda pusat.


Penurunan Rumus Turunan dengan Deret Taylor Misalkan diberikan titik-titik (𝑥𝑖 , 𝑓𝑖 ), 𝑖 = 0,1,2,3, … , 𝑛, yang dalam hal ini 𝑥𝑖 = 𝑥0 + 𝑖∆𝑥 dan 𝑓𝑖 = 𝑓(𝑥𝑖 ). Selanjutnya akan dihitung 𝑓 ′ (𝑥) yang dalam hal ini 𝑥 = 𝑥0 + 𝑠∆𝑥, 𝑠 ∈ 𝑅 dengan ketiga pendekatan yang disebutkan di atas (maju, mundur, pusat). 1.

Hampiran Beda Maju

Uraikan 𝑓(𝑥𝑖+1 ) di sekitar 𝑥𝑖 : (𝑥𝑖+1 − 𝑥𝑖 ) ′ (𝑥𝑖+1 − 𝑥𝑖 )2 ′′ 𝑓(𝑥𝑖+1 ) = 𝑓(𝑥𝑖 ) + 𝑓 (𝑥𝑖 ) + 𝑓 (𝑥𝑖 ) + ⋯ 1! 2! dengan mensubtitusikan (𝑥𝑖+1 − 𝑥𝑖 ) = ∆𝑥 dan penulisan 𝑓(𝑥𝑖+1 ) dapat ditulis 𝑓𝑖+1 diperoleh 𝑓𝑖+1 = 𝑓𝑖 + ∆𝑥𝑓𝑖 ′ +

∆𝑥 2 ′′ 𝑓 +⋯ 2 𝑖

(2.13)

atau dapat ditulis ∆𝑥𝑓𝑖 ′ = 𝑓𝑖+1 − 𝑓𝑖 −

∆𝑥 2 ′′ 𝑓 −⋯ 2 𝑖

kedua ruas dibagi dengan ∆𝑥 sehingga diperoleh 𝑓𝑖 ′ =

𝑓𝑖+1 − 𝑓𝑖 ∆𝑥 ′′ − 𝑓 −⋯ ∆𝑥 2 𝑖

karena

∆𝑥 2

𝑓𝑖 ′′ − ⋯ merupakan bilangan yang sangat kecil dan tidak begitu

mempengaruhi nilai 𝑓𝑖 ′ sehingga dapat ditulis


𝑓𝑖 ′ =

𝑓𝑖+1 − 𝑓𝑖 + 𝑂(∆𝑥) ∆𝑥

yang dalam hal ini, 𝑂(∆𝑥) =

∆𝑥 2

𝑓 ′′ (𝑡), 𝑥𝑖 < 𝑡 < 𝑥𝑖+1 .

Untuk nilai-nilai 𝑓 di 𝑥0 dan 𝑥1 persamaan rumusnya menjadi: 𝑓0 ′ =

𝑓1 − 𝑓0 + 𝑂(∆𝑥). ∆𝑥

Dalam hal ini 𝑂(∆𝑥) =

∆𝑥 2

𝑓 ′′ (𝑡), 𝑥𝑖 < 𝑡 < 𝑥𝑖+1 menyatakan penurunan numeris

secara beda maju memiliki tingkat keakuratan tingkat satu atau ditulis 𝑂(∆𝑥).

2. Hampiran Beda Mundur Uraikan 𝑓(𝑥𝑖−1 ) di sekitar 𝑥𝑖 : 𝑓(𝑥𝑖−1 ) = 𝑓(𝑥𝑖 ) +

(𝑥𝑖+1 − 𝑥𝑖 ) ′ (𝑥𝑖+1 − 𝑥𝑖 )2 ′′ 𝑓 (𝑥𝑖 ) + 𝑓 (𝑥𝑖 ) + ⋯ 1! 2!

dengan mensubtitusikan (𝑥𝑖+1 − 𝑥𝑖 ) = ∆𝑥 dan penulisan 𝑓(𝑥𝑖−1 ) dapat ditulis 𝑓𝑖−1 diperoleh 𝑓𝑖−1

∆𝑥 2 ′′ = 𝑓𝑖 − ∆𝑥𝑓𝑖 + 𝑓 −⋯ 2 𝑖 ′

(2.14)

atau dapat ditulis ∆𝑥 2 ′′ ∆𝑥𝑓𝑖 = 𝑓𝑖 − 𝑓𝑖−1 + 𝑓 −⋯ 2 𝑖 ′

kedua ruas dibagi dengan ∆𝑥 sehingga diperoleh 𝑓𝑖 ′ =

𝑓𝑖 − 𝑓𝑖−1 ∆𝑥 ′′ − 𝑓 +⋯ ∆𝑥 2 𝑖

karena

∆𝑥 2

𝑓𝑖 ′′ − ⋯ merupakan bilangan yang sangat kecil dan tidak begitu

mempengaruhi nilai 𝑓𝑖 ′ sehingga dapat ditulis


𝑓𝑖 ′ =

𝑓𝑖 − 𝑓𝑖−1 + 𝑂(∆𝑥) ∆𝑥

yang dalam hal ini, 𝑂(∆𝑥) = −

∆𝑥 2

𝑓 ′′ (𝑡), 𝑥𝑖−1 < 𝑡 < 𝑥𝑖 .

Untuk nilai-nilai 𝑓 di 𝑥0 dan 𝑥−1 persamaan rumusnya menjadi: 𝑓0 ′ =

𝑓0 − 𝑓−1 + 𝑂(∆𝑥) ∆𝑥

Dalam hal ini 𝑂(∆𝑥) = −

∆𝑥 2

𝑓 ′′ (𝑡), 𝑥𝑖+1 < 𝑡 < 𝑥𝑖 menyatakan penurunan numeris

secara beda mundur memiliki tingkat keakuratan tingkat satu atau ditulis 𝑂(∆𝑥).

3.

Hampiran Beda Pusat

Kurangkan persamaan (2.13) dengan persamaan (2.14) diperoleh: 𝑓𝑖+1 − 𝑓𝑖−1

∆𝑥 2 ′′ ∆𝑥 2 ′′ ′ = 𝑓𝑖 + ∆𝑥𝑓𝑖 + 𝑓 + ⋯ − (𝑓𝑖 − ∆𝑥𝑓𝑖 + 𝑓 − ⋯) 2 𝑖 2 𝑖 ′

dengan menggunakan operasi penjumlahan dan pengurangan diperoleh 𝑓𝑖+1 − 𝑓𝑖−1 = 2∆𝑥𝑓𝑖 ′ +

∆𝑥 3 ′′′ 𝑓 +⋯ 3 𝑖

atau dapat ditulis 2∆𝑥𝑓𝑖 ′ = 𝑓𝑖+1 − 𝑓𝑖−1 −

∆𝑥 3 ′′′ 𝑓 −⋯ 3 𝑖

Kedua ruas dibagi dengan 2∆𝑥 sehingga diperoleh 𝑓𝑖+1 − 𝑓𝑖−1 ∆𝑥 2 ′′′ 𝑓𝑖 = − 𝑓 −⋯ 2∆𝑥 6 𝑖 ′

karena

∆𝑥 2 6

𝑓𝑖 ′′′ − ⋯ merupakan bilangan yang sangat kecil dan tidak begitu

mempengaruhi nilai 𝑓𝑖 ′ sehingga dapat ditulis 𝑓𝑖 ′ =

𝑓𝑖+1 − 𝑓𝑖−1 + 𝑂(∆𝑥 2 ), 2∆𝑥


yang dalam hal ini, 𝑂(∆𝑥 2 ) = −

∆𝑥 2 6

𝑓 ′′′ (𝑡), 𝑥𝑖−1 < 𝑡 < 𝑥𝑖+1 .

Untuk nilai-nilai 𝑓 di 𝑥−1 dan 𝑥1 persamaan rumusnya menjadi: 𝑓0 ′ =

𝑓1 − 𝑓−1 + 𝑂(∆𝑥 2 ) 2∆𝑥

Dalam hal ini 𝑂(∆𝑥 2 ) = −

∆𝑥 2 6

𝑓 ′′′ (𝑡), 𝑥𝑖−1 < 𝑡 < 𝑥𝑖+1 menyatakan penurunan

numeris secara beda pusat yang memiliki tingkat keakuratan tingkat dua atau ditulis 𝑂(∆𝑥 2 ). Perhatikan bahwa hampiran beda pusat lebih baik daripada dua hampiran sebelumnya, sebab orde galatnya adalah 𝑂(∆𝑥 2 ).

Menentukan Orde Galat Pada penurunan rumus turunan numeris dengan deret Taylor, rumus galat dalam penurunan rumus turunan numeris tersebut dapat langsung diperoleh. Tetapi dengan polinom interpolasi harus dicari rumus galat tersebut dengan bantuan deret Taylor. Contoh 2.9 Tentukan rumus galat dan orde dari rumus turunan numeris hampiran beda pusat: 𝑓 ′ (𝑥0 ) =

𝑓1 − 𝑓−1 +𝐸 2∆𝑥

Nyatakan 𝐸 (galat) sebagai ruas kiri persamaan, lalu ekspansi rusa kanan dengan deret Taylor di sekitar 𝑥0 : 𝐸 = 𝑓 ′ (𝑥0 ) −

𝑓1 − 𝑓−1 2∆𝑥


1

=𝑓0 ′ − 2∆𝑥 [(𝑓0 + ∆𝑥𝑓0 ′ + ∆𝑥 2 2

𝑓0 ′′ −

∆𝑥 3 6

1

=− =−

∆𝑥 2 6 ∆𝑥 2 6

2

𝑓0 ′′ +

∆𝑥 3 6

𝑓0 ′′′ + ⋯ ) − (𝑓0 − ∆𝑥𝑓0 ′ +

𝑓0 ′′′ + ⋯ )]

=𝑓0 ′ − 2∆𝑥 (2∆𝑥𝑓0 ′ + =𝑓0 ′ − 𝑓0 ′ −

∆𝑥 2

∆𝑥 2 6

∆𝑥 3 3

𝑓0 ′′′ + ⋯ )

𝑓0 ′′′ + ⋯

𝑓0 ′′′ + ⋯ 𝑓0 ′′′ , 𝑥−1 < 𝑡 < 𝑥1

= 𝑂(∆𝑥 2 ). Jadi, hampiran beda pusat memiliki galat 𝐸 = −

∆𝑥 2 6

𝑓0 ′′′ , 𝑥−1 < 𝑡 < 𝑥1 , dengan

orde 𝑂(∆𝑥 2 ).

F. Karakteristik Persamaan Gelombang Air Dangkal Dipandang persamaan gelombang air dangkal ℎ𝑡 + (𝑢ℎ)𝑥 = 0

(2.15)

1 (ℎ𝑢)𝑡 + (𝑢2 ℎ + 𝑔ℎ2 )𝑥 = 0. 2

(2.16)

Gabungan persamaan (2.15) dan (2.16) dalam suatu sistem persamaan gelombang air dangkal yaitu

[

ℎ𝑢 ℎ ] + [ 2 1 2] = 0 ℎ𝑢 + 𝑔ℎ ℎ𝑢 𝑡 2 𝑥

(2.17)

Jika diasumsikan ℎ dan 𝑢 halus (smooth), maka persamaan (2.16) dapat disederhanakan dengan memperluas turunan-turunannya dan menggunakan (2.15)


untuk menggantikan bentuk ℎ𝑡 . Kemudian dengan menghilangkan beberapa bentuk, persamaan (2.16) menjadi 1 𝑢𝑡 + [ 𝑢2 + 𝑔ℎ]𝑥 = 0. 2

(2.18)

Pada persamaan (2.15) dan (2.18) memiliki bentuk yang bergantung dengan konstanta 𝑔. Bentuk tersebut dapat disubtitusi dengan variabel 𝜑 = 𝑔ℎ. Sehingga sistem persamaan air dangkal menjadi 𝑢2 𝑢 [𝜑 ] + [ 2 + 𝜑 ] = 0 𝑡 𝑢𝜑 𝑥

(2.19)

Sistem persamaan tersebut ekuivalen dengan sistem persamaan (2.17) untuk solusi yang halus. Namun ada catatan penting bahwa manipulasi yang dilakukan di atas bergantung pada kehalusan pada masalah. Kedua sistem dari hukum konservasi tidak ekuivalen dalam menghitung shock waves. Sistem yang tepat untuk digunakan adalah sistem persamaan (2.17) yang berasal dari persamaan integral asli. Untuk mempelajari shock waves digunakan persamaan (2.17) dan diambil 𝑞2 ℎ𝑢 𝑞1 2 ℎ 1 𝑞(𝑥, 𝑡) = [ ] = [𝑞 ] , 𝑓(𝑞) = [ 2 1 2 ] = [(𝑞2 ) ]. + 𝑔(𝑞1 )2 ℎ𝑢 + 𝑔ℎ 2 ℎ𝑢 𝑞1 2 2 Untuk solusi halus (smooth), persamaan tersebut dapat ditulis secara ekuivalen dalam bentuk quasilinear 𝑞𝑡 + 𝑓 ′ (𝑞)𝑞𝑥 = 0 Dengan matriks Jacobian 𝑓 ′ (𝑞) adalah


𝑓

′ (𝑞)

0 2 (𝑞 ) 1 =[ 2 + 𝑔(𝑞1 )2 𝑞1 2

1 𝑞2 ] = [ 2 0 −𝑢 + 𝑔ℎ 2 𝑞1

1 ]. 2𝑢

(2.20)

Nilai eigen dari 𝑓 ′ (𝑞) adalah 𝜆1 = 𝑢 − √𝑔ℎ, 𝜆2 = 𝑢 + √𝑔ℎ

(2.21)

Dengan vektor eigen 𝑟1 = [

1 1 ] , 𝑟2 = [ ]. 𝑢 − √𝑔ℎ 𝑢 + √𝑔ℎ

(2.22)

Nilai eigen dan vektor eigen adalah fungsi 𝑞 untuk sistem nonlinear. Jika diinginkan gelombang dengan amplitudo yang sangat kecil, maka persamaan (2.17) dapat dilinearkan terlebih dahulu untuk mendapatkan sistem yang linear. Keterangan lengkap tentang karakteristik persamaan air dangkal dapat ditemukan dalam buku karangan LeVeque (2004).


BAB III METODE NUMERIS UNTUK PERSAMAAN GELOMBANG AIR DANGKAL

Dalam bab ini akan dijelaskan metode volume hingga, metode beda hingga grid kolokasi dan metode beda hingga grid selang-seling. Metode tersebut digunakan untuk menyelesaikan masalah bendungan air bobol, terkait dengan persamaan gelombang air dangkal.

A.

Solusi Eksak Persamaan Gelombang Air Dangkal Galat perhitungan pada penurunan numeris diperoleh dengan mengurangkan

solusi numeris dengan solusi eksak. Berikut adalah solusi eksak yang akan digunakan dalam perhitungan simulasi numeris pada MATLAB: ℎ1 , ℎ3 = ℎ(𝑥) = ℎ2 = {

jika 𝑥 ≤ −𝑡√𝑔ℎ1

4 𝑥 2 (√𝑔ℎ1 − ) , 9𝑔 2𝑡

jika − 𝑡√𝑔ℎ1 < 𝑥 ≤ 𝑡(𝑢2 − √𝑔ℎ2 )

ℎ0 8𝜉̇ 2 (√1 + − 1) , 2 𝑔ℎ0

jika 𝑡(𝑢2 − √𝑔ℎ2 ) < 𝑥 < 𝑡𝜉̇ jika 𝑥 ≥ 𝑡𝜉̇

ℎ0 ,

dan

28


0, 2 𝑥 𝑢3 = (√𝑔ℎ1 + ) , 3 𝑡 2

𝑢(𝑥) =

𝑢2 = 𝜉̇ −

𝑔ℎ0 8𝜉̇ 2 (√1 + ), 𝑔ℎ0 4𝜉̇

{

0,

jika 𝑥 ≤ −𝑡√𝑔ℎ1 jika − 𝑡√𝑔ℎ1 < 𝑥 ≤ 𝑡(𝑢2 − √𝑔ℎ2 ) jika 𝑡(𝑢2 − √𝑔ℎ2 ) < 𝑥 < 𝑡𝜉̇ jika 𝑥 ≥ 𝑡𝜉̇

dengan ℎ(𝑥) adalah kedalaman air pada titik 𝑥 dan 𝑢(𝑥) adalah kecepatan air pada titik 𝑥. Notasi 𝜉̇ adalah konstanta kecepatan shock untuk 𝑡 > 0, yaitu 1 2

𝜉̇ = 2√𝑔ℎ1 +

𝑔ℎ0 8𝜉̇ 2 8𝜉̇ 2 (1 + √1 + ) − [2𝑔ℎ0 √1 + − 2𝑔ℎ0 ] . 𝑔ℎ0 𝑔ℎ0 4𝜉̇

Solusi eksak ini diambil dari Thesis karangan Mungkasi dengan judul Finite Volume Methods for the One Dimensional Shallow Water Equations (2008).

B.

Penurunan Persamaan Gelombang Air Dangkal Satu Dimensi Persamaan gelombang air dangkal dideskripsikan dari gerak fluida. Ada dua

jenis gerak fluida yang dideskripsikan, yaitu Langrangian dan Eulerian. Deskripsi Langrangian berpusat pada partikel individu, dan pergerakannya diamati sebagai fungsi dari waktu. Posisi, kecepatan dan percepatan setiap partikel dinotasikan dengan 𝑠(𝑥0 , 𝑡), 𝑢(𝑥0 , 𝑡), dan 𝑎(𝑥0 , 𝑡), kemudian kuantitasnya misalnya massa, momentum dan energi dapat dihitung. Pada kasus ini 𝑥0 merupakan titik awal atau penamaan partikel. Deskripsi Eulerian merupakan sebuah alternatif yang dapat diikuti setiap partikel fluida secara terpisah. Kemudian dilakukan pengamatan untuk kecepatan partikel yang melewati setiap titik identifikasi pada domain ruang yang diamati.


Laju perubahan kecepatan ketika partikel melewati setiap titik dapat diamati dengan 𝜕𝑢(𝑥,𝑡) 𝜕𝑥

, dan perubahan kecepatan terhadap waktu pada setiap titik tertentu dapat

diamati oleh

𝜕𝑢(𝑥,𝑡) 𝜕𝑡

. Dalam deskripsi Eulerian, sifat aliran (seperti kecepatan)

merupakan fungsi dari ruang dan waktu. Persamaan air dangkal ini terdiri dari dua persamaan. Persamaan pertama diturunkan dari hukum konservasi massa dan persamaan kedua diturunkan dari hukum konservasi momentum. Berikut ini akan diuraikan penurunan persamaan gelombang air dangkal atau biasa disebut Shallow Water Wave Equations. a. Hukum Kekekalan Massa Hukum kekekalan massa berarti massa tersebut tidak dapat diciptakan atau dimusnahkan. Hal ini berarti massa total pada keseluruhan sistem sama setiap saat. Terdapat beberapa asumsi yang terlibat dalam penurunan persamaan hukum kekekalan massa. Pertama, aliran air diasumsikan tenang artinya tidak ada gangguan dari luar dan kecepatannya diabaikan. Kedua, densitas 𝜌 air pada setiap titik adalah konstan sehingga air mampat. Selain itu diasumsikan bahwa tempat air kedap atau tertutup rapat karena massa adalah kekal. Oleh karena itu, massa pada setiap volume kontrol (yaitu volume tertentu atau kolam air yang diamati) hanya dapat berubah ketika aliran melintasi batas-batas volume kontrol. Secara umum, aliran air dapat diilustrasikan pada Gambar 1.1. Notasi yang digunakan yaitu 𝑥 menyatakan variabel jarak sepanjang aliran air, 𝑡 menyatakan variabel waktu, 𝑧(𝑥) adalah topografi tanah, ℎ(𝑥, 𝑡) adalah kedalaman air di titik 𝑥 dan pada waktu 𝑡, 𝑤(𝑥, 𝑡) = 𝑧(𝑥) + ℎ(𝑥, 𝑡) adalah ketinggian air mutlak disebut


stage, dan 𝑢(𝑥, 𝑡) adalah kecepatan aliran air di titik 𝑥 dan pada waktu 𝑡. Massa total 𝑚 pada air di setiap volume kontrol [𝑥1 , 𝑥2 ] ditentukan oleh 𝑥2

(3.1)

𝑚 = ∫ 𝜌ℎ(𝑥, 𝑡)𝑑𝑥 𝑥1

Pernyataan ini dapat diperoleh sebagai berikut. Kepadatan massa terhadap kedalaman 𝜌̅ di sebarang titik (𝑥, 𝑡) adalah 𝜌ℎ(𝑥, 𝑡) yang dapat dihitung dengan mengintegralkan 𝜌 dari 𝑧(𝑥) ke 𝑤(𝑥, 𝑡), yaitu 𝜌̅ (𝑥, 𝑡) = ∫

𝑤(𝑥,𝑡)

𝜌 𝑑𝑦

𝑧(𝑥)

= 𝜌ℎ(𝑥, 𝑡). Akibatnya, pengintegralan 𝜌ℎ(𝑥, 𝑡) dari 𝑥1 ke 𝑥2 mengarah ke massa total di volume kontrol seperti yang dinyatakan dalam (3.1). Tingkatan aliran air yang melewati setiap titik (𝑥, 𝑡) terhadap kedalaman air disebut flux massa 𝑓1 , yaitu 𝑓1 = 𝜌̅ (𝑥, 𝑡)𝑢(𝑥, 𝑡)

(3.2)

= 𝜌ℎ(𝑥, 𝑡)𝑢(𝑥, 𝑡) Dengan menggunakan (3.2) dan asumsi bahwa massa dapat berubah hanya karena aliran yang melewati batas volume kontrol, dapat ditentukan bahwa 𝑥2

𝑥2

∫ 𝜌ℎ(𝑥, 𝑡 + ∆𝑡)𝑑𝑥 = ∫ 𝜌ℎ(𝑥, 𝑡)𝑑𝑥 𝑥1

𝑥1

𝑡+∆𝑡

+∫ 𝑡

(3.3)

𝑡+∆𝑡

𝜌ℎ(𝑥1 , 𝑠)𝑢(𝑥1 , 𝑠)𝑑𝑠 − ∫ 𝑡

𝜌ℎ(𝑥2 , 𝑠)𝑢(𝑥2 , 𝑠)𝑑𝑠

berlaku untuk setiap volume kontrol. Hal ini berarti bahwa massa pada setiap langkah 𝑡 + ∆𝑡 adalah sama dengan massa pada waktu 𝑡 ditambah pergerakan flux


yang masuk dan dikurangi dengan flux yang keluar dari volume kontrol selama periode ∆𝑡. Ilustrasi dari kontinuitas massa ditunjukkan pada Gambar 3.1.

𝑢 𝑢̅

𝑥1

𝑥2

Gambar 3.1: Aliran air yang masuk dan keluar dari volume kontrol. Misalkan ∆𝑥 dan ∆𝑡 adalah kuantitas yang sangat kecil, yaitu ∆𝑥 = 𝑥2 − 𝑥1 . Dengan menggunakan perluasan Taylor, persamaan (3.3) dapat ditulis 𝜌ℎ(𝑥, 𝑡 + ∆𝑡)∆𝑥 = 𝜌ℎ(𝑥, 𝑡)∆𝑥 + 𝜌ℎ (𝑥 − − 𝜌ℎ (𝑥 +

∆𝑥 ∆𝑥 , 𝑡) 𝑢 (𝑥 − , 𝑡) ∆𝑡 2 2

∆𝑥 ∆𝑥 , 𝑡) 𝑢 (𝑥 + , 𝑡) ∆𝑡 + 𝑂((∆𝑡)3 ) + 𝑂((∆𝑥)3 ). 2 2

Dengan mengabaikan bentuk 𝑂((∆𝑡)3 ) dan 𝑂((∆𝑥)3 ) persamaan terakhir di atas ekuivalen dengan 𝜌ℎ(𝑥, 𝑡 + ∆𝑡) − 𝜌ℎ(𝑥, 𝑡) =− ∆𝑡

(𝜌ℎ𝑢)|

(𝑥+

∆𝑥 ,𝑡) 2

− (𝜌ℎ𝑢)|

(𝑥−

∆𝑥 ,𝑡) 2

(3.4)

∆𝑥

Persamaan (3.4) dibagi dengan 𝜌 kemudian ∆𝑥 dan ∆𝑡 diaproksimasikan menuju nol sehingga persamaan (3.4) menjadi ℎ𝑡 + (𝑢ℎ)𝑥 = 0. Persamaan (3.5) disebut persamaan hukum kekekalan massa.

(3.5)


b. Hukum Kekekalan Momentum Hukum kedua Newton menyatakan bahwa perubahan momentum dari suatu sistem sama dengan total gaya yang bekerja. Berdasarkan hukum Newton tersebut, maka dapat ditulis 𝐹=

𝑑𝑝 . 𝑑𝑡

Gaya 𝐹 didefinisikan sebagai laju perubahan momentum 𝑝 terhadap waktu 𝑡. Momentum total dari perpindahan air pada volume kontrol dari 𝑥1 ke 𝑥2 pada waktu 𝑡 dinotasikan dengan 𝑝(𝑡), yaitu

𝑝(𝑡) = ∫

𝑥2 (𝑡)

𝜌ℎ(𝑥, 𝑡)𝑢(𝑥, 𝑡)𝑑𝑥

(3.6)

𝑥1 (𝑡)

Dengan mengasumsikan tekanan hidrostatik, gaya pada titik 𝑥1 dan 𝑥2 di atas kedalaman air pada waktu 𝑡 adalah 𝐹1 (𝑡) =

1 𝜌𝑔ℎ2 (𝑥1 (𝑡), 𝑡) 2

1 𝐹2 (𝑡) = − 𝜌𝑔ℎ2 (𝑥2 (𝑡), 𝑡) 2 Dengan 𝑔 > 0 adalah konstanta yang menyatakan percepatan gravitasi. Lebih lanjut, gaya pada ∆𝑧 seperti yang ditunjukkan pada Gambar 3.2, yaitu ∆𝐹3 = −𝜌𝑔ℎ(𝑥, 𝑡)∆𝑧 Atau dapat ditulis dengan ∆𝐹3 = −𝜌𝑔ℎ(𝑥, 𝑡)

∆𝑧 ∆𝑥 ∆𝑥


Dan karena gaya melalui dasar volume kontrol maka 𝑥2

𝐹3 = ∫ −𝜌𝑔ℎ(𝑥, 𝑡)𝑧𝑥 𝑑𝑥. 𝑥1

Oleh karena itu, gaya total di atas volume kontrol dinyatakan dengan 𝐹 yang merupakan jumlahan dari 𝐹1 , 𝐹2 , dan 𝐹3 , yaitu 𝐹=

𝑥2 1 1 𝑑𝑧 𝜌𝑔ℎ2 (𝑥1 (𝑡), 𝑡) − 𝜌𝑔ℎ2 (𝑥2 (𝑡), 𝑡) − ∫ 𝜌𝑔ℎ(𝑥, 𝑡) 𝑑𝑥 2 2 𝑑𝑥 𝑥1

(3.7)

Turunan pertama dari 𝑝 terhadap 𝑡 adalah 𝑑𝑝 𝑑 𝑥2 = ∫ 𝜌ℎ(𝑥, 𝑡)𝑢(𝑥, 𝑡)𝑑𝑥 𝑑𝑡 𝑑𝑡 𝑥1 Dengan menggunakan aturan Leibniz, turunan hasil integral pada persamaan terakhir di atas dapat ditulis 𝑥2 𝑑𝑝 𝜕 =∫ 𝜌ℎ(𝑥, 𝑡)𝑢(𝑥, 𝑡)𝑑𝑥 + 𝜌ℎ(𝑥2 (𝑡), 𝑡)𝑢2 (𝑥2 (𝑡), 𝑡) 𝑑𝑡 𝑥1 𝜕𝑡

(3.8)

− 𝜌ℎ(𝑥1 (𝑡), 𝑡)𝑢2 (𝑥1 (𝑡), 𝑡) Menurut hukum kedua Newton tentang gerak, hasil dari persamaan (3.8) sama dengan persamaan (3.7). Oleh karena itu, untuk periode-∆𝑡 dapat ditulis


𝑡+∆𝑡

∫ 𝑡

∫

𝑥2 (𝑡)

𝑥1 (𝑡)

(𝜌ℎ𝑢)𝑡 𝑑𝑥 𝑑𝑡 + ∫ 𝑡

𝑡+∆𝑡

−∫ 𝑡

𝑡

𝑡

(3.9)

𝑡+∆𝑡 1 1 𝜌𝑔ℎ2 (𝑥1 (𝑡), 𝑡) 𝑑𝑡 ∫ 𝜌𝑔ℎ2 (𝑥2 (𝑡), 𝑡) 𝑑𝑡 2 2 𝑡

𝑡+∆𝑡

−∫

𝜌ℎ(𝑥2 (𝑡), 𝑡)𝑢2 (𝑥2 (𝑡), 𝑡)𝑑𝑡

𝜌ℎ(𝑥1 (𝑡), 𝑡)𝑢2 (𝑥1 (𝑡), 𝑡)𝑑𝑡

𝑡+∆𝑡

=∫

𝑡+∆𝑡

∫

𝑥2 (𝑡)

𝑥1 (𝑡)

𝜌𝑔ℎ(𝑥, 𝑡)𝑧𝑥 𝑑𝑥 𝑑𝑡

Dengan cara yang sama seperti (3.3), persamaan (3.9) dapat ditulis 1 (ℎ𝑢)𝑡 + (ℎ𝑢2 + 𝑔ℎ2 ) = −𝑔ℎ𝑧𝑥 2 𝑥 Yang biasa disebut dengan persamaan kekekalan momentum. Dengan demikian, persamaan gelombang air dangkal seringkali disebut sistem Saint-Venant. Persamaan tersebut dapat ditulis sebagai dua persamaan simultan ℎ𝑡 + (𝑢ℎ)𝑥 = 0 1 { (ℎ𝑢)𝑡 + (ℎ𝑢2 + 𝑔ℎ2 ) = −𝑔ℎ𝑧𝑥 2 𝑥

(3.10)

di sini variabel x menyatakan arah aliran air.

C. Masalah Bendungan Bobol Diketahui persamaan gelombang air dangkal dengan topografi horizontal 𝑞𝑡 + 𝑓(𝑞)𝑥 = 0 Dengan kuantitas dan flux berturut-turut adalah

(3.11)


𝑢ℎ ℎ 1 𝑞 = [ ] dan 𝑓(𝑞) = [ 2 ] 𝑢 ℎ + 𝑔ℎ2 𝑢ℎ 2 dengan 𝑥 adalah variabel ruang, 𝑡 adalah variabel waktu, ℎ = ℎ(𝑥, 𝑡) adalah kedalaman air, 𝑢 = 𝑢(𝑥, 𝑡) adalah kecepatan air dan 𝑔 = 9,81 adalah percepatan gravitasi. Semua kuantitas diasumsikan dalam satuan SI. Akan disimulasikan solusi masalah bendungan bobol dengan metode volume hingga Lax-Friedrics, metode beda hingga grid kolokasi dan metode beda hingga grid selang-seling dengan menggunakan MATLAB (kondisi awal adalah "air yang tenang" seperti yang ditunjukkan pada Gambar 3.2). Pada kasus ini dianggap dinding bendungan ditarik ke atas secara instan.

Permukaan air

Bendungan air

ℎ1 = 10

Permukaan air ℎ0 = 4 𝑥

Gambar 3.2: Bendungan air Dinding bendungan air berada di titik 𝑥 = 0 dan kedalaman awal air adalah ℎ(𝑥, 0) = {

dan kecepatan awal aliran air adalah

ℎ1 , jika 𝑥 < 0 ℎ0 , jika 𝑥 > 0


𝑢(𝑥, 0) = 0, untuk semua 𝑥. Diambil domain ruang [−5,5]. Simulasi pada program dihentikan pada 𝑡 = 0.2. Pada kasus ini diasumsikan massa jenis konstan, tidak ada turbulen dan fluida ideal.

D. Metode Volume Hingga Lax-Friedrichs Pada bagian ini dibahas mengenai skema metode volume hingga, perhitungan flux secara numeris dalam metode volume hingga dan solusi numeris metode volume hingga Lax-Friedrichs. 1.1. Skema Metode Volume Hingga Persamaan diferensial parsial hukum kekekalan yang bersifat hiperbolik adalah 𝑞𝑡 + 𝑓(𝑞)𝑥 = 0 atau ditulis 𝜕 𝜕𝑡

𝜕

𝑞(𝑥, 𝑡) + 𝜕𝑥 𝑓(𝑞(𝑥, 𝑡)) = 0.

Misalkan domain pada ruang didiskretkan menjadi sebanyak berhingga kontrol volume (interval) atau sel sebagai berikut: 𝑥𝑖−3

𝑥𝑖−1

2

𝑥𝑖+3

𝑥𝑖+1

2

2

2

𝑥𝑖−1

𝑥𝑖

𝑥𝑖+1

dengan ∆𝑥 = 𝑥𝑖−1 − 𝑥𝑖 atau ∆𝑥 = 𝑥𝑖+1 − 𝑥𝑖−1 . 2

2

Domain waktu didiskretkan menjadi 𝑡 𝑛 = 𝑛 ∙ ∆𝑡


dengan 𝑛 = 0,1,2,3, … Selanjutnya misalkan 𝑄𝑖𝑛 adalah pendekatan dari rata-rata volume kuantitas 𝑞(𝑥, 𝑡) dalam interval ruang ke- 𝑖 di waktu 𝑡 𝑛 , yaitu 𝑄𝑖𝑛 ≈

1 𝑥𝑖+12 ∫ 𝑞(𝑥, 𝑡 𝑛 ) 𝑑𝑥. ∆𝑥 𝑥 1 𝑖− 2

𝑛 Misalkan pula 𝐹𝑖+ 1 adalah pendekatan dari rata-rata debit material (flux) 𝑓(𝑞(𝑥, 𝑡)) 2

di titik 𝑥𝑖+1 dalam interval waktu [𝑡 𝑛 , 𝑡 𝑛+1 ], yaitu 2

𝑛+1

𝐹𝑛 1 𝑖+ 2

1 𝑡 ≈ ∫ ∆𝑡 𝑡 𝑛

𝑓 (𝑞 (𝑥𝑖+1 , 𝑡)) 𝑑𝑡. 2

Dari hukum kekekalan, laju perubahan kuantitas (massa) dinyatakan oleh 𝑑 𝑥𝑖+12 ∫ 𝑞(𝑥, 𝑡 𝑛 ) 𝑑𝑥 = − [𝑓 (𝑞 (𝑥𝑖+1 , 𝑡)) − 𝑓 (𝑞 (𝑥𝑖−1 , 𝑡))]. 𝑑𝑡 𝑥 1 2 2 𝑖−

2

Dengan nilai-nilai pendekatan diperoleh

𝑄𝑖𝑛+1 − 𝑄𝑖𝑛 =− ∆𝑡

𝐹𝑛 1 − 𝐹𝑛 1 𝑖+

𝑖−

2

2

∆𝑥

atau ditulis

𝑄𝑖𝑛+1 = 𝑄𝑖𝑛 −

∆𝑡 𝑛 (𝐹 1 − 𝐹 𝑛 1 ) 𝑖− ∆𝑥 𝑖+2 2

Persamaan tersebut merupakan skema volume hingga untuk 𝑞𝑡 + 𝑓(𝑞)𝑥 = 0. Skema metode volume hingga tersebut konsisten dengan skema metode beda hingga sebab dari skema metode volume hingga diperoleh


𝑄𝑖𝑛+1 − 𝑄𝑖𝑛 =− ∆𝑡

𝐹𝑛 1 − 𝐹𝑛 1 𝑖+

𝑖−

2

2

∆𝑥

atau

𝑄𝑖𝑛+1 − 𝑄𝑖𝑛 + ∆𝑡

𝐹𝑛 1 − 𝐹𝑛 1 𝑖+

𝑖−

2

∆𝑥

2

=0

yang merupakan suatu bentuk diskrit dari 𝑞𝑡 + 𝑓(𝑞)𝑥 = 0. 1.2. Perhitungan Flux Secara Numeris dalam Metode Volume Hingga Dipandang persamaan diferensial parsial berbentuk hukum kekekalan 𝑞𝑡 + 𝑓(𝑞)𝑥 = 0. 𝑛 𝑛 Misalkan 𝑄𝑖𝑛 ≈ 𝑞(𝑥𝑖 , 𝑡 𝑛 ) dan 𝐹𝑖+ 1 = 𝑓 (𝑞 (𝑥 1 , 𝑡 )). 𝑖+ 2

2

Skema metode volume hingga untuk persamaan diferensial parsial di atas adalah


∆𝑡 𝑛 (𝐹 1 − 𝐹 𝑛 1 ). 𝑖− ∆𝑥 𝑖+2 2

Diketahui nilai 𝑄𝑖𝑛 yaitu kuantitas numeris di semua titik 𝑥𝑖 dan pada waktu 𝑡 𝑛 . Oleh karena itu, flux di titik 𝑥𝑖 pada waktu 𝑡 𝑛 juga diketahui, yaitu 𝐹𝑖𝑛 ≈ 𝑓(𝑞(𝑥𝑖 , 𝑡 𝑛 )) ≈ 𝑓(𝑄𝑖𝑛 ).


Metode Stabil dan Tidak Stabil Metode numeris dikatakan stabil artinya bahwa galat yang muncul pada setiap iterasi tidak membesar terlalu cepat pada iterasi-iterasi berikutnya. Jika galat yang muncul pada suatu iterasi membesar menuju tak terhingga pada iterasi-iterasi berikutnya maka metodenya disebut tidak stabil. Teori kestabilan tidak dibahas lebih lanjut dalam skripsi ini. Teori kestabilan dapat dilihat dalam buku-buku referensi misalnya LeVeque (1992, 2004). Di setiap titik 𝑥𝑖+1 , flux dapat dihitung menggunakan beberapa pendekatan. 2

1. Flux tak stabil 𝐹𝑛 1 ≈ 𝑖+

2

1 𝑛 [𝑓(𝑄𝑖𝑛 ) + 𝑓(𝑄𝑖+1 )] 2

Sehingga skema metode volume hingga menjadi 𝑄𝑖𝑛+1 = 𝑄𝑖𝑛 −

∆𝑡 1 1 1 1 𝑛 𝑛 ) − 𝑓(𝑄𝑖−1 ) − 𝑓(𝑄𝑖𝑛 )) ( 𝑓(𝑄𝑖𝑛 ) + 𝑓(𝑄𝑖+1 ∆𝑥 2 2 2 2


∆𝑡 𝑛 𝑛 ) − 𝑓(𝑄𝑖−1 )) (𝑓(𝑄𝑖+1 2∆𝑥

Namun, skema metode volume hingga ini tidak stabil. 2. Flux Lax-Friedrichs Skema Lax-Friedrichs memodifikasi skema metode volume hingga tak stabil di atas, yaitu 𝑄𝑖𝑛+1 ≈ Sehingga diperoleh

1 𝑛 𝑛 (𝑄 + 𝑄𝑖−1 ) 2 𝑖+1


1 𝑛 ∆𝑡 𝑛 𝑛 𝑛 )− [𝑓(𝑄𝑖+1 ) − 𝑓(𝑄𝑖−1 )] 𝑄𝑖𝑛+1 = (𝑄𝑖+1 + 𝑄𝑖−1 2 2∆𝑥 Skema Lax-Friedrichs ini stabil untuk ∆𝑡 yang cukup kecil. 3. Flux Upwind Metode upwind cocok untuk metode yang sudah diketahui arah rambatan gelombangnya. Contoh jika akan diselesaikan persamaan diferensial parsial 𝑞𝑡 + 𝑐𝑞𝑥 = 0 dengan 𝑐 konstan positif (arah rambat gelombang ke kanan) 𝐹 𝑛 1 ≈ 𝑓(𝑞(𝑥𝑖 , 𝑡 𝑛 )) 𝑖+

2

= 𝑐𝑞(𝑥𝑖 , 𝑡 𝑛 ) ≈ 𝑐𝑄𝑖𝑛 𝑛 𝑛 dan 𝐹𝑖− 1 ≈ 𝑐𝑄𝑖−1 . 2

Dengan demikian skema upwind adalah 𝑄𝑖𝑛+1 = 𝑄𝑖𝑛 − = 𝑄𝑖𝑛 −

∆𝑡 𝑛 (𝐹 1 − 𝐹 𝑛 1 ) 𝑖− ∆𝑥 𝑖+2 2

∆𝑡 𝑛 (𝑐𝑄𝑖𝑛 − 𝑐𝑄𝑖−1 ) ∆𝑥

= 𝑄𝑖𝑛 − 𝑐

∆𝑡 𝑛 𝑛 (𝑄 − 𝑄𝑖−1 ). ∆𝑥 𝑖

Solusi Numeris Metode Volume Hingga Lax-Friedrichs Masalah persamaan gelombang air dangkal dapat diselesaikan dengan menggunakan metode volume hingga. Dipandang sistem persamaan air dangkal (3.11) yaitu


ℎ𝑡 + [𝑢ℎ]𝑥 = 0

(3.12)

1 [𝑢ℎ]𝑡 + [𝑢2 ℎ + 𝑔ℎ2 ]𝑥 = 0. 2

(3.13)

dan

Persamaan (3.12) mempunyai skema metode volume hingga sebagai berikut: 𝑄𝑖𝑛+1 = 𝑄𝑖𝑛 −

∆𝑡 𝑛 (𝐹 1 − 𝐹 𝑛 1 ) 𝑖− ∆𝑥 𝑖+2 2

Dengan menggunakan flux Lax-Friedrichs diperoleh 𝐹 𝑛 1 = 𝑓(𝑄𝑖𝑛 ) 𝑖+

2

dengan asumsi solusi dan kecepatannya positif. Jadi, jika diketahui persamaan (3.12) maka dalam metode volume hingga diperoleh 𝑞 = ℎ dan 𝑓(𝑞) = 𝑢ℎ. 𝑛 𝑛 Sekarang dari persamaan (3.12) dicari flux 𝐹𝑖+ 1 dan 𝐹 1 . 𝑖− 2

𝐹𝑛 1 𝑖+

2

𝐹𝑛 1 𝑖−

2

=

1 ∆𝑥 𝑛 𝑛 [𝑓(𝑄𝑖𝑛 ) + 𝑓(𝑄𝑖+1 )] − (𝑄𝑖+1 − 𝑄𝑖𝑛 ) 2 2∆𝑡

=

1 ∆𝑥 𝑛 [(𝑢ℎ)𝑛𝑖 + (𝑢ℎ)𝑛𝑖+1 ] − (ℎ − ℎ𝑖𝑛 ) 2 2∆𝑡 𝑖+1

=

1 𝑛 𝑛 ∆𝑥 𝑛 𝑛 𝑛 (𝑢𝑖 ℎ𝑖 + 𝑢𝑖+1 )− (ℎ − ℎ𝑖𝑛 ). ℎ𝑖+1 2 2∆𝑡 𝑖+1

=

1 ∆𝑥 𝑛 𝑛 [𝑓(𝑄𝑖−1 ) + 𝑓(𝑄𝑖𝑛 )] − (𝑄 𝑛 − 𝑄𝑖−1 ) 2 2∆𝑡 𝑖

=

1 ∆𝑥 𝑛 𝑛 [(𝑢ℎ)𝑛𝑖−1 + (𝑢ℎ)𝑛𝑖 ] − (ℎ𝑖 − ℎ𝑖−1 ) 2 2∆𝑡

2


=

1 𝑛 𝑛 ∆𝑥 𝑛 𝑛 (𝑢𝑖−1 ℎ𝑖−1 + 𝑢𝑖𝑛 ℎ𝑖𝑛 ) − (ℎ𝑖 − ℎ𝑖−1 ). 2 2∆𝑡

Persamaan (3.13) memiliki skema metode volume hingga sebagai berikut 𝑄𝑖𝑛+1 = 𝑄𝑖𝑛 −

∆𝑡 𝑛 (𝐹 1 − 𝐹 𝑛 1 ) 𝑖− ∆𝑥 𝑖+2 2

Dengan menggunakan flux Lax-Friedrichs diperoleh 𝐹 𝑛 1 = 𝑓(𝑄𝑖𝑛 ) 𝑖+

2

dengan asumsi solusi dan kecepatannya positif. Jadi, jika diketahui persamaan (3.13)

maka

dalam

metode

volume

hingga

diperoleh

𝑞 = 𝑢ℎ

dan

1

𝑓(𝑞) = 𝑢2 ℎ + 2 𝑔ℎ2 . 𝑛 𝑛 Sekarang dari persamaan (3.13) dicari flux 𝐹𝑖+ 1 dan 𝐹 1 . 𝑖− 2

𝐹𝑛 1 𝑖+

2

=

2

1 ∆𝑥 𝑛 [𝑓(𝑄𝑖𝑛 ) + 𝑓(𝑄𝑖+1 )] − (𝑄 𝑛 − 𝑄𝑖𝑛 ) 2 2∆𝑡 𝑖+1

1 1 2 𝑛 1 2 𝑛 ∆𝑥 2 2 ((𝑢ℎ)𝑛𝑖+1 − (𝑢ℎ)𝑛𝑖 ) = [(𝑢 ℎ + 𝑔ℎ ) + (𝑢 ℎ + 𝑔ℎ ) ] − 2 2 2 2∆𝑡 𝑖 𝑖+1 =

1 1 1 𝑛 [((𝑢2 )𝑛𝑖 ℎ𝑖𝑛 + 𝑔(ℎ2 )𝑛𝑖 ) + ((𝑢2 )𝑛𝑖+1 ℎ𝑖+1 + 𝑔(ℎ2 )𝑛𝑖+1 )] 2 2 2 −

𝐹𝑛 1 𝑖−

2

∆𝑥 𝑛 𝑛 (𝑢 ℎ − 𝑢𝑖𝑛 ℎ𝑖𝑛 ). 2∆𝑡 𝑖+1 𝑖+1

=

1 ∆𝑥 𝑛 𝑛 [𝑓(𝑄𝑖−1 ) + 𝑓(𝑄𝑖𝑛 )] − (𝑄𝑖𝑛 − 𝑄𝑖−1 ) 2 2∆𝑡

=

𝑛 𝑛 1 1 1 ∆𝑥 ((𝑢ℎ)𝑛𝑖 − (𝑢ℎ)𝑛𝑖−1 ) [(𝑢2 ℎ + 𝑔ℎ2 ) + (𝑢2 ℎ + 𝑔ℎ2 ) ] − 2 2 2 2∆𝑡 𝑖−1 𝑖


=

1 1 1 𝑛 [((𝑢2 )𝑛𝑖−1 ℎ𝑖−1 + 𝑔(ℎ2 )𝑛𝑖−1 ) + ((𝑢2 )𝑛𝑖 ℎ𝑖𝑛 + 𝑔(ℎ2 )𝑛𝑖 )] 2 2 2 −

∆𝑥 𝑛 𝑛 𝑛 𝑛 (𝑢 ℎ − 𝑢𝑖−1 ). ℎ𝑖−1 2∆𝑡 𝑖 𝑖

Hasil simulasi penyelesaian model gelombang air dangkal dengan metode volume hingga dengan menggunakan program MATLAB ditunjukkan dalam Gambar 3.3. Pada hasil simulasi persamaan air dangkal, program simulasi berhenti pada saat 𝑡 = 0.2. Kedalaman awal air pada program ini adalah ℎ1 = 10, ℎ0 = 4.

Gambar 3.3: Hasil simulasi penyelesaian model gelombang air dangkal dengan metode volume hingga Lax-Friedrichs.


E. Metode Beda Hingga Grid Kolokasi Metode beda hingga digunakan dalam masalah komputasi numerik. Metode beda hingga dibagi menjadi tiga bagian, yaitu beda maju, beda mundur, dan beda pusat. Solusi Numeris Metode Beda Hingga Grid Kolokasi Diketahui nilai awal kedalaman air ℎ0 = 4 dan ℎ1 = 10, kecepatan air 𝑢0 = 𝑢1 = 0, dan percepatan gravitasi 𝑔 = 9.81. Banyaknya langkah untuk ruang adalah 𝑁 dengan 𝑁 ∈ ℤ+ dan index 𝑖 untuk menunjukkan ruang yaitu 1 ≤ 𝑖 ≤ 𝑁 + 1. Langkah ukuran dilambangkan dengan ∆=

𝑥𝑏 −𝑥𝑎 𝑁

. Jadi, 𝑥𝑖 = 𝑥𝑎 + (𝑖 − 1)∆, 𝑥1 = 𝑥𝑎 dan 𝑥𝑁+1 = 𝑥𝑏 , dan banyaknya

interval kecil 𝑁 adalah 𝑁 − 1 simpul (node). Begitu juga untuk waktu, jumlah langkah untuk waktu adalah 𝑀 dengan 𝑀 ∈ ℤ+ . Untuk 𝑡 𝑛 waktu, 𝑛 ∈ ℤ+ : 𝑛 − 1 < 𝑛 ≤ 𝑀. Berikut akan dihitung persamaan beda yaitu beda mundur untuk waktu dan beda pusat untuk ruang.

Beda mundur untuk variabel waktu 𝒕 Dipandang persamaan (3.12) dan (3.13). Kedua persamaan tersebut ditulis secara sistematis dan diasumsikan nilai ℎ(𝑥, 𝑡) dan 𝑢(𝑥, 𝑡) pada waktu 𝑡 𝑛−1 kemudian dicari nilai pada waktu 𝑡 𝑛 . 𝜕 ℎ(𝑥𝑖 , 𝑡 𝑛 ) − ℎ(𝑥𝑖 , 𝑡 𝑛−1 ) ℎ(𝑥, 𝑡) ≈ 𝜕𝑡 ∆𝑡 atau dapat ditulis


≈

ℎ𝑖𝑛 − ℎ𝑖𝑛−1 . ∆𝑡

dan 𝜕 𝜕 𝜕 [ℎ(𝑥, 𝑡)𝑢(𝑥, 𝑡)] = ℎ(𝑥, 𝑡) 𝑢(𝑥, 𝑡) + 𝑢(𝑥, 𝑡) ℎ(𝑥, 𝑡) 𝜕𝑡 𝜕𝑡 𝜕𝑡 atau dapat ditulis ≈ ℎ(𝑥𝑖 , 𝑡 𝑛−1 )

𝜕 𝜕 𝑢(𝑥𝑖 , 𝑡 𝑛−1 ) + 𝑢(𝑥𝑖 , 𝑡 𝑛−1 ) ℎ(𝑥𝑖 , 𝑡 𝑛−1 ) 𝜕𝑡 𝜕𝑡

dengan melakukan diskritisasi lengkap diperoleh

≈

𝑛 𝑛−1 𝑢𝑖 ℎ𝑖

𝑛 𝑛−1 − 𝑢𝑖𝑛−1 𝑛−1 ℎ𝑖 − ℎ𝑖 + 𝑢𝑖 . ∆𝑡 ∆𝑡

Beda pusat untuk variabel ruang 𝒙 Dipandang persamaan (3.12) dan (3.13). Kedua persamaan tersebut ditulis secara sistematis dan diasumsikan nilai ℎ(𝑥, 𝑡) dan

𝑢(𝑥, 𝑡) pada ruang 𝑥𝑖−1

kemudian dicari nilai pada waktu 𝑥𝑖 . 𝜕 𝜕 𝜕 [𝑢(𝑥, 𝑡)ℎ(𝑥, 𝑡)] = 𝑢(𝑥, 𝑡) ℎ(𝑥, 𝑡) + ℎ(𝑥, 𝑡) 𝑢(𝑥, 𝑡) 𝜕𝑥 𝜕𝑥 𝜕𝑥 ≈ 𝑢(𝑥𝑖 , 𝑡 𝑛−1 )

≈

dan

𝜕 𝜕 ℎ(𝑥𝑖 , 𝑡 𝑛 ) + ℎ(𝑥𝑖 , 𝑡 𝑛−1 ) 𝑢(𝑥𝑖 , 𝑡 𝑛 ) 𝜕𝑥 𝜕𝑥

𝑛 𝑛−1 ℎ𝑖+1 𝑢𝑖

𝑛 𝑛 𝑛 − ℎ𝑖−1 𝑛−1 𝑢𝑖+1 − 𝑢𝑖−1 + ℎ𝑖 2∆𝑥 2∆𝑥


𝜕 1 𝜕 1 𝜕 [ℎ(𝑥, 𝑡)𝑢2 (𝑥, 𝑡)] + 𝑔 ℎ2 (𝑥, 𝑡) [ℎ(𝑥, 𝑡)𝑢2 (𝑥, 𝑡) + 𝑔ℎ2 (𝑥, 𝑡)] = 𝜕𝑥 2 𝜕𝑥 2 𝜕𝑥 dengan 𝜕 𝜕 𝜕 [ℎ(𝑥, 𝑡)𝑢2 (𝑥, 𝑡)] = ℎ(𝑥, 𝑡) 𝑢2 (𝑥, 𝑡) + 𝑢2 (𝑥, 𝑡) ℎ(𝑥, 𝑡) 𝜕𝑥 𝜕𝑥 𝜕𝑥 = 2ℎ(𝑥, 𝑡)𝑢(𝑥, 𝑡)

≈ 2ℎ(𝑥𝑖 , 𝑡 𝑛−1 )𝑢(𝑥𝑖 , 𝑡 𝑛−1 )

≈

𝜕 𝜕 𝑢(𝑥, 𝑡) + 𝑢2 (𝑥, 𝑡) ℎ(𝑥, 𝑡) 𝜕𝑥 𝜕𝑥 𝜕 𝜕 𝑢(𝑥𝑖 , 𝑡 𝑛 ) + 𝑢2 (𝑥𝑖 , 𝑡 𝑛−1 ) ℎ(𝑥𝑖 , 𝑡 𝑛 ) 𝜕𝑥 𝜕𝑥

𝑛 𝑛−1 𝑛−1 𝑢𝑖+1 2ℎ𝑖 𝑢𝑖

𝑛 𝑛 𝑛 − 𝑢𝑖−1 𝑛−1 2 ℎ𝑖+1 − 𝑢𝑖−1 + (𝑢𝑖 ) 2∆𝑥 2∆𝑥

dan untuk bentuk kedua 1 𝜕 2 1 𝜕 𝑔 ℎ (𝑥, 𝑡) = 2 𝑔ℎ(𝑥, 𝑡) ℎ(𝑥, 𝑡) 2 𝜕𝑥 2 𝜕𝑥 ≈ 𝑔ℎ𝑖𝑛−1

𝑛 𝑛 ℎ𝑖+1 − ℎ𝑖−1 . 2∆𝑥

Jadi, persamaan beda dari (3.12) dan (3.13) menjadi ℎ𝑖𝑛 − ℎ𝑖𝑛−1 ∆𝑡 𝑛 𝑛−1 𝑛 𝑛−1 𝑛−1 𝑢𝑖 − 𝑢𝑖 𝑛−1 ℎ𝑖 − ℎ𝑖 ℎ + 𝑢 𝑖 [ 𝑖 ] ∆𝑡 ∆𝑡 𝑛 𝑛−1 ℎ𝑖+1 𝑢𝑖

𝑛 𝑛 𝑛 − ℎ𝑖−1 𝑛−1 𝑢𝑖+1 − 𝑢𝑖−1 + ℎ𝑖 0 2∆𝑥 2∆𝑥 + =[ ] 𝑛 𝑛 𝑛 𝑛 𝑛 𝑛 ℎ𝑖+1 − 𝑢𝑖−1 ℎ𝑖+1 − ℎ𝑖−1 0 𝑛−1 𝑛−1 𝑢𝑖+1 − 𝑢𝑖−1 + (𝑢𝑖𝑛−1 )2 + 𝑔ℎ𝑖𝑛−1 [2ℎ𝑖 𝑢𝑖 ] 2∆𝑥 2∆𝑥 2∆𝑥

atau dapat ditulis berdasarkan persamaan konservasi massa dalam bentuk diskrit


𝑛 𝑛 𝑛 𝑛 ℎ𝑖𝑛 − ℎ𝑖𝑛−1 ℎ𝑖+1 − ℎ𝑖−1 𝑢𝑖+1 − 𝑢𝑖−1 + 𝑢𝑖𝑛−1 + ℎ𝑖𝑛−1 =0 ∆𝑡 2∆𝑥 2∆𝑥

(3.16)

dan konservasi momentum dalam bentuk diskrit 𝑛 𝑛−1 𝑢𝑖 ℎ𝑖

𝑛 𝑛−1 𝑛 𝑛 − 𝑢𝑖𝑛−1 𝑛−1 ℎ𝑖 − ℎ𝑖 𝑛−1 𝑛−1 𝑢𝑖+1 − 𝑢𝑖−1 + 𝑢𝑖 + 2ℎ𝑖 𝑢𝑖 ∆𝑡 ∆𝑡 2∆𝑥

+ (𝑢𝑖𝑛−1 )2

(3.17)

𝑛 𝑛 𝑛 𝑛 ℎ𝑖+1 − 𝑢𝑖−1 ℎ𝑖+1 − ℎ𝑖−1 + 𝑔ℎ𝑖𝑛−1 =0 2∆𝑥 2∆𝑥

Persamaan (3.16) dikali dengan 𝑢𝑖𝑛−1 diperoleh 𝑛 𝑛−1 ℎ𝑖 𝑢𝑖

𝑛 𝑛 𝑛 𝑛 − ℎ𝑖𝑛−1 𝑛−1 2 ℎ𝑖+1 − ℎ𝑖−1 𝑛−1 𝑛−1 𝑢𝑖+1 − 𝑢𝑖−1 + (𝑢𝑖 ) + 𝑢𝑖 ℎ𝑖 =0 ∆𝑡 2∆𝑥 2∆𝑥

(3.18)

Persamaan (3.17) dikurang (3.18) diperoleh 𝑛 𝑛−1 𝑢𝑖 ℎ𝑖

𝑛 𝑛 𝑛 𝑛 − 𝑢𝑖𝑛−1 𝑛−1 𝑛−1 𝑢𝑖+1 − 𝑢𝑖−1 𝑛−1 ℎ𝑖+1 − ℎ𝑖−1 + ℎ𝑖 𝑢𝑖 + 𝑔ℎ𝑖 =0 ∆𝑡 2∆𝑥 2∆𝑥

Sehingga diperoleh dua persamaan

𝑢𝑖𝑛−1

𝑛 𝑛 𝑛 𝑛 ℎ𝑖𝑛 − ℎ𝑖𝑛−1 ℎ𝑖+1 − ℎ𝑖−1 𝑢𝑖+1 − 𝑢𝑖−1 + (𝑢𝑖𝑛−1 )2 + 𝑢𝑖𝑛−1 ℎ𝑖𝑛−1 ∆𝑡 2∆𝑥 2∆𝑥

(3.19)

=0 ℎ𝑖𝑛−1

𝑛 𝑛 𝑛 𝑛 𝑢𝑖𝑛 − 𝑢𝑖𝑛−1 𝑢𝑖+1 − 𝑢𝑖−1 ℎ𝑖+1 − ℎ𝑖−1 + ℎ𝑖𝑛−1 𝑢𝑖𝑛−1 + 𝑔ℎ𝑖𝑛−1 =0 ∆𝑡 2∆𝑥 2∆𝑥

Persamaan konservasi massa pada persamaan (3.19) dibagi dengan 𝑢𝑖𝑛−1 , diperoleh 𝑛 𝑛 𝑛 𝑛 ℎ𝑖𝑛 − ℎ𝑖𝑛−1 𝑛−1 ℎ𝑖+1 − ℎ𝑖−1 𝑛−1 𝑢𝑖+1 − 𝑢𝑖−1 + 𝑢𝑖 + ℎ𝑖 =0 ∆𝑡 2∆𝑥 2∆𝑥

ℎ𝑖𝑛−1

𝑛 𝑛 𝑛 𝑛 𝑢𝑖𝑛 − 𝑢𝑖𝑛−1 𝑢𝑖+1 − 𝑢𝑖−1 ℎ𝑖+1 − ℎ𝑖−1 + ℎ𝑖𝑛−1 𝑢𝑖𝑛−1 + 𝑔ℎ𝑖𝑛−1 ∆𝑡 2∆𝑥 2∆𝑥

=0

(3.20)


Persamaan (3.20) dikali dengan 2∆𝑥∆𝑡 diperoleh 𝑛 𝑛 𝑛 𝑛 ) + ∆𝑡ℎ𝑖𝑛−1 (𝑢𝑖+1 ) 2∆𝑥(ℎ𝑖𝑛 − ℎ𝑖𝑛−1 ) + ∆𝑡𝑢𝑖𝑛−1 (ℎ𝑖+1 − ℎ𝑖−1 − 𝑢𝑖−1

(3.21)

=0 𝑛 𝑛 ) 2∆𝑥ℎ𝑖𝑛−1 (𝑢𝑖𝑛 − 𝑢𝑖𝑛−1 ) + ∆𝑡ℎ𝑖𝑛−1 𝑢𝑖𝑛−1 (𝑢𝑖+1 − 𝑢𝑖−1 𝑛 𝑛 ) = 0. + ∆𝑡𝑔ℎ𝑖𝑛−1 (ℎ𝑖+1 − ℎ𝑖−1

Dengan menggunakan sifat distributif pada perkalian, persamaan (3.21) menjadi 𝑛 𝑛 𝑛 𝑛 2∆𝑥ℎ𝑖𝑛 − 2∆𝑥ℎ𝑖𝑛−1 + ∆𝑡𝑢𝑖𝑛−1 ℎ𝑖+1 − ∆𝑡𝑢𝑖𝑛−1 ℎ𝑖−1 + ∆𝑡ℎ𝑖𝑛−1 𝑢𝑖+1 − ∆𝑡ℎ𝑖𝑛−1 𝑢𝑖−1

=0 𝑛 𝑛 2∆𝑥ℎ𝑖𝑛−1 𝑢𝑖𝑛 − 2∆𝑥ℎ𝑖𝑛−1 𝑢𝑖𝑛−1 + ∆𝑡ℎ𝑖𝑛−1 𝑢𝑖𝑛−1 𝑢𝑖+1 − ∆𝑡ℎ𝑖𝑛−1 𝑢𝑖𝑛−1 𝑢𝑖−1 𝑛 𝑛 + ∆𝑡𝑔ℎ𝑖𝑛−1 ℎ𝑖+1 − ∆𝑡𝑔ℎ𝑖𝑛−1 ℎ𝑖−1 =0

Sehingga persamaan (3.16) menjadi 𝑛 𝑛 𝑛 𝑛 2∆𝑥ℎ𝑖𝑛 + ∆𝑡𝑢𝑖𝑛−1 ℎ𝑖+1 − ∆𝑡𝑢𝑖𝑛−1 ℎ𝑖−1 + ∆𝑡ℎ𝑖𝑛−1 𝑢𝑖+1 − ∆𝑡ℎ𝑖𝑛−1 𝑢𝑖−1

(3.22)

= 2∆𝑥ℎ𝑖𝑛−1 dan (3.17) menjadi 𝑛 𝑛 𝑛 2∆𝑥ℎ𝑖𝑛−1 𝑢𝑖𝑛 + ∆𝑡ℎ𝑖𝑛−1 𝑢𝑖𝑛−1 𝑢𝑖+1 − ∆𝑡ℎ𝑖𝑛−1 𝑢𝑖𝑛−1 𝑢𝑖−1 + ∆𝑡𝑔ℎ𝑖𝑛−1 ℎ𝑖+1

(3.23)

𝑛 − ∆𝑡𝑔ℎ𝑖𝑛−1 ℎ𝑖−1 = 2∆𝑥ℎ𝑖𝑛−1 𝑢𝑖𝑛−1

Diambil nilai 𝑁 = 4 , kemudian dicari solusi awal untuk persamaan linear diskritisasi sebagai berikut: Untuk 𝑖 = 1


2∆𝑥ℎ1𝑛 + ∆𝑡𝑢1𝑛−1 ℎ2𝑛 − ∆𝑡𝑢1𝑛−1 ℎ0𝑛 + ∆𝑡ℎ1𝑛−1 𝑢2𝑛 − ∆𝑡ℎ1𝑛−1 𝑢0𝑛

(3.24)

= 2∆𝑥ℎ1𝑛−1 dan 2∆𝑥ℎ1𝑛−1 𝑢1𝑛 + ∆𝑡ℎ1𝑛−1 𝑢1𝑛−1 𝑢2𝑛 − ∆𝑡ℎ1𝑛−1 𝑢1𝑛−1 𝑢0𝑛 + ∆𝑡𝑔ℎ1𝑛−1 ℎ2𝑛

(3.25)

− ∆𝑡𝑔ℎ1𝑛−1 ℎ0𝑛 = 2∆𝑥ℎ1𝑛−1 𝑢1𝑛−1 Untuk 𝑖 = 2 2∆𝑥ℎ2𝑛 + ∆𝑡𝑢2𝑛−1 ℎ3𝑛 − ∆𝑡𝑢2𝑛−1 ℎ1𝑛 + ∆𝑡ℎ2𝑛−1 𝑢3𝑛 − ∆𝑡ℎ2𝑛−1 𝑢1𝑛

(3.26)


(3.27)

− ∆𝑡𝑔ℎ2𝑛−1 ℎ1𝑛 = 2∆𝑥ℎ2𝑛−1 𝑢2𝑛−1 Untuk 𝑖 = 3 2∆𝑥ℎ3𝑛 + ∆𝑡𝑢3𝑛−1 ℎ4𝑛 − ∆𝑡𝑢3𝑛−1 ℎ2𝑛 + ∆𝑡ℎ3𝑛−1 𝑢4𝑛 − ∆𝑡ℎ3𝑛−1 𝑢2𝑛

(3.28)


(3.29)

− ∆𝑡𝑔ℎ3𝑛−1 ℎ2𝑛 = 2∆𝑥ℎ3𝑛−1 𝑢3𝑛−1 Variabel 𝑢1𝑛 , 𝑢2𝑛 , 𝑢3𝑛 , ℎ1𝑛 , ℎ2𝑛 , ℎ3𝑛 tidak diketahui dan enam persamaan berikut dapat disederhanakan menjadi: Untuk 𝑖 = 1


2∆𝑥ℎ1𝑛 + ∆𝑡𝑢1𝑛−1 ℎ2𝑛 + ∆𝑡ℎ1𝑛−1 𝑢2𝑛 − ∆𝑡ℎ1𝑛−1 𝑢0𝑛 = 2∆𝑥ℎ1𝑛−1 +

(3.30)

∆𝑡𝑢1𝑛−1 ℎ0𝑛 +−∆𝑡ℎ1𝑛−1 𝑢0𝑛 2∆𝑥ℎ1𝑛−1 𝑢1𝑛 + ∆𝑡ℎ1𝑛−1 𝑢1𝑛−1 𝑢2𝑛 + ∆𝑡𝑔ℎ1𝑛−1 ℎ2𝑛

(3.31)

= 2∆𝑥ℎ1𝑛−1 𝑢1𝑛−1 + ∆𝑡𝑔ℎ1𝑛−1 ℎ0𝑛 + ∆𝑡ℎ1𝑛−1 𝑢1𝑛−1 𝑢0𝑛 Untuk 𝑖 = 2 2∆𝑥ℎ2𝑛 + ∆𝑡𝑢2𝑛−1 ℎ3𝑛 − ∆𝑡𝑢2𝑛−1 ℎ1𝑛 + ∆𝑡ℎ2𝑛−1 𝑢3𝑛 − ∆𝑡ℎ2𝑛−1 𝑢1𝑛

(3.32)


(3.33)

− ∆𝑡𝑔ℎ2𝑛−1 ℎ1𝑛 = 2∆𝑥ℎ2𝑛−1 𝑢2𝑛−1 Untuk 𝑖 = 3 2∆𝑥ℎ3𝑛 − ∆𝑡𝑢3𝑛−1 ℎ2𝑛 − ∆𝑡ℎ3𝑛−1 𝑢2𝑛

(3.34)

= 2∆𝑥ℎ3𝑛−1 − ∆𝑡𝑢3𝑛−1 ℎ4𝑛 − ∆𝑡ℎ3𝑛−1 𝑢4𝑛 dan 2∆𝑥ℎ3𝑛−1 𝑢3𝑛 − ∆𝑡ℎ3𝑛−1 𝑢3𝑛−1 𝑢2𝑛 − ∆𝑡𝑔ℎ3𝑛−1 ℎ2𝑛

(3.35)

= 2∆𝑥ℎ3𝑛−1 𝑢3𝑛−1 − ∆𝑡𝑔ℎ3𝑛−1 ℎ4𝑛 − ∆𝑡ℎ3𝑛−1 𝑢3𝑛−1 𝑢4𝑛 Jika diketahui syarat batas kecepatan awal 𝑢(𝑥𝑎 , 𝑡) = 𝑢(𝑥𝑏 , 𝑡) = 0, dan ketinggian awal ℎ(𝑥𝑎 , 𝑡) = 1, ℎ(𝑥𝑏 , 𝑡) = 1 maka diperoleh matriks 0

2∆𝑥ℎ𝑛−1 1 −∆𝑡ℎ𝑛−1 2 𝑛−1 𝑛−1 −∆𝑡ℎ2 𝑢2 0 [

0

∆𝑡ℎ𝑛−1 1 𝑛−1 ∆𝑡ℎ𝑛−1 1 𝑢1 0

2∆𝑥ℎ𝑛−1 2 −∆𝑡𝑢𝑛−1 3 𝑛−1 −∆𝑡ℎ𝑛−1 3 𝑢3

0

2∆𝑥

0

0

∆𝑡ℎ𝑛−1 2 𝑛−1 𝑛−1 ∆𝑡ℎ2 𝑢2

−∆𝑡𝑢𝑛−1 2 −∆𝑡𝑔ℎ𝑛−1 2

0

0

2∆𝑥ℎ𝑛−1 3

0

∆𝑡𝑢𝑛−1 1 ∆𝑡𝑔ℎ𝑛−1 1

0

2∆𝑥

∆𝑡𝑢𝑛−1 2 ∆𝑡𝑔ℎ𝑛−1 2

0

−∆𝑡ℎ𝑛−1 3 −∆𝑡𝑔ℎ𝑛−1 3

0

2∆𝑥 0

]


2∆𝑥ℎ1𝑛−1 + ∆𝑡𝑢1𝑛−1 𝑢1𝑛 2∆𝑥ℎ1𝑛−1 𝑢1𝑛−1 + ∆𝑡𝑔ℎ1𝑛−1 𝑢2𝑛 𝑢3𝑛 2∆𝑥ℎ2𝑛−1 = ℎ1𝑛 2∆𝑥ℎ2𝑛−1 𝑢2𝑛−1 𝑛 ℎ2 2∆𝑥ℎ3𝑛−1 − ∆𝑡𝑢3𝑛−1 𝑛 [ℎ3 ] [2∆𝑥ℎ3𝑛−1 𝑢3𝑛−1 − ∆𝑡𝑔ℎ3𝑛−1 ] Misalkan 𝐴 0

=

2∆𝑥ℎ𝑛−1 1 −∆𝑡ℎ𝑛−1 2 𝑛−1 𝑛−1 −∆𝑡ℎ2 𝑢2 0 (

0

∆𝑡ℎ𝑛−1 1 𝑛−1 𝑛−1 ∆𝑡ℎ1 𝑢1 0

2∆𝑥ℎ𝑛−1 2 −∆𝑡𝑢𝑛−1 3 𝑛−1 𝑛−1 −∆𝑡ℎ3 𝑢3

0

2∆𝑥

0

0

∆𝑡ℎ𝑛−1 2 𝑛−1 𝑛−1 ∆𝑡ℎ2 𝑢2

−∆𝑡𝑢𝑛−1 2 −∆𝑡𝑔ℎ𝑛−1 2

0

0

2∆𝑥ℎ𝑛−1 3

0

∆𝑡𝑢𝑛−1 1 ∆𝑡𝑔ℎ𝑛−1 1

0

2∆𝑥

∆𝑡𝑢𝑛−1 2 ∆𝑡𝑔ℎ𝑛−1 2

0

−∆𝑡ℎ𝑛−1 3 −∆𝑡𝑔ℎ𝑛−1 3

0

2∆𝑥 0

𝑢1𝑛 𝑢2𝑛 𝑢3𝑛 𝑋= , ℎ1𝑛 ℎ2𝑛 (ℎ3𝑛 ) 2∆𝑥ℎ1𝑛−1 + ∆𝑡𝑢1𝑛−1 2∆𝑥ℎ1𝑛−1 𝑢1𝑛−1 + ∆𝑡𝑔ℎ1𝑛−1 2∆𝑥ℎ2𝑛−1 𝑏= , 2∆𝑥ℎ2𝑛−1 𝑢2𝑛−1 2∆𝑥ℎ3𝑛−1 − ∆𝑡𝑢3𝑛−1 (2∆𝑥ℎ3𝑛−1 𝑢3𝑛−1 − ∆𝑡𝑔ℎ3𝑛−1 ) Sekarang enam persamaan tersebut memiliki bentuk linear 𝐴𝑋 = 𝑏 yang dapat diselesaikan. Dalam program MATLAB digunakan syarat batas untuk kedalaman awal air ℎ1 = 10, ℎ0 = 4. Hasil simulasi penyelesaian model gelombang air dangkal dengan

)


metode beda hingga grid kolokasi dengan menggunakan program MATLAB ditunjukkan dalam Gambar 3.4. Hasil simulasi program berhenti pada waktu 𝑡 = 0.2.

Gambar 3.4: Hasil simulasi penyelesaian model gelombang air dangkal dengan metode beda hingga grid kolokasi.

F. Metode Beda Hingga Grid Selang-Seling Sekarang persamaan air dangkal akan diselesaikan dengan metode beda hingga grid selang-seling kemudian solusi numerisnya dibandingkan dengan metode volume hingga dan metode beda hingga grid kolokasi yang sudah dilakukan di atas. Dipandang persamaan air dangkal


ℎ𝑢 ℎ [ ] + [ 2 1 2 ] = 0. ℎ𝑢 + 𝑔ℎ ℎ𝑢 𝑡 2 𝑥 Perhitungan domain 0 < 𝑥 < 𝑋, 𝑡 > 0 dengan grid selang-seling pada domain ruang 𝑥1 = 0, 𝑥1 , … , 𝑥𝑖−1 , 𝑥𝑖 , 𝑥𝑖+1 , … , 𝑥𝑁 2

2

1 𝑥 +2

2

= 𝑋. Metode beda hingga selang-

seling adalah pendekatan dua persamaan (3.12) dan (3.13) pada grid berbeda. Pada rumus selang-seling, nilai variabel ℎ dihitung pada grid bilangan bulat 𝑥𝑖 , 𝑖 = 1, … , 𝑁𝑥 dan 𝑢 dihitung pada grid pecahan 𝑥𝑖+1 , 𝑖 = 1, … , 𝑁𝑥 . Selanjutnya, hukum 2

kekekalan massa pada persamaan (3.12) dihitung dengan pendekatan pada grid [𝑥𝑖−1 , 𝑥𝑖+1 ] dan hukum kekekalan momentum pada persamaan (3.13) dihitung 2

2

dengan pendekatan pada grid [𝑥𝑗 , 𝑥𝑗+1 ]. Solusi Numeris Metode Beda Hingga Grid Selang-Seling Pendekatan konsisten terhadap ruang dari (3.12) dan (3.13) adalah ∆ℎ𝑖 =− ∆𝑡 ∆(ℎ𝑖+1 𝑢𝑖+1 ) 2

2

∆𝑡

𝑞𝑖+1 − 𝑞𝑖−1 2

(3.36)

2

∆𝑥

2 1 ℎ𝑖+1 − ℎ𝑖2 𝑞𝑖+1 𝑢̂𝑖+1 − 𝑞𝑖 𝑢̂𝑖 =− 𝑔 − 2 ∆𝑥 ∆𝑥

(3.37)

dengan 𝑞𝑖+1 = ℎ̂𝑖+1 𝑢𝑖+1 , ℎ𝑖+1 = 2

2

2

2

1 1 (ℎ𝑖 + ℎ𝑖+1 ), 𝑞𝑖 = (𝑞 1 + 𝑞 1 ) 𝑖− 2 2 𝑖+2 2

(3.38)

Dengan menggunakan metode upwind orde satu, pendekatan nilai dari ℎ𝑖 ,

ℎ̂𝑖+1 = { ℎ𝑖+1 , 2 dan

jika 𝑢𝑖+1 ≥ 0 2

jika 𝑢𝑖+1 < 0 2

(3.39)


𝑢̂𝑖 = {

𝑢𝑖−1 ,

jika 𝑞𝑖 ≥ 0

𝑢𝑖+1 ,

jika 𝑞𝑖 < 0

2 2

(3.40)

Dalam melakukan perhitungan, variabel ℎ dihitung pada grid ℎ𝑖 , 𝑖 = 1,2, … , 𝑁𝑥 dan variabel 𝑢 dihitung pada grid 𝑢𝑖+1 , 𝑖 = 0,1, … , 𝑁𝑥 . Dengan menerapkan beda maju 2

untuk persamaan (3.36) dan beda mundur untuk persamaan (3.37), persamaan (3.36) dan (3.37) menjadi

ℎ𝑖𝑛+1

− ∆𝑡

𝑛+1 𝑛 𝑛 ℎ𝑛+1 1 𝑢 1 − ℎ 1𝑢 𝑖+

2

𝑖+

𝑖+

2

∆𝑡

2

𝑖+

1 2

ℎ𝑖𝑛

=−

𝑞𝑖+1 − 𝑞𝑖−1 2

(3.41)

2

∆𝑥

2 1 ℎ𝑖+1 − ℎ𝑖2 𝑞𝑖+1 𝑢̂𝑖+1 − 𝑞𝑖 𝑢̂𝑖 =− 𝑔 − 2 ∆𝑥 ∆𝑥

(3.42)

Persamaan (3.41) dan (3.42) dihitung dengan menggunakan program MATLAB. Hasil simulasi gelombang air dangkal dengan metode beda hingga grid selangseling ditunjukkan dalam Gambar 3.5. Pada simulasi penyelesaian persamaan gelombang air dangkal berikut, program berhenti pada saat 𝑡 = 0.2 dengan kedalaman awal air adalah ℎ1 = 10 dan ℎ0 = 4.


Gambar 3.5: Hasil simulasi model gelombang air dangkal dengan metode beda hingga grid selang-seling.


BAB IV PERBANDINGAN BEBERAPA METODE NUMERIS DALAM PENYELESAIAN MODEL GELOMBANG AIR DANGKAL

Pada bagian ini akan dibahas hasil-hasil simulasi numeris untuk tiga metode, yaitu metode volume hingga, metode beda hingga grid kolokasi dan metode beda hingga grid selang-seling. Simulasi numeris untuk masing-masing metode dilakukan menggunakan program MATLAB dengan konsidi awal ditunjukkan pada Gambar 4.1.

Permukaan air

Bendungan air

ℎ1 = 10

Permukaan air ℎ0 = 4 𝑥

Gambar 4.1: Bendungan air dengan kondisi awal adalah air yang tenang. Kedalaman air di sebelah kiri bendungan adalah ℎ1 = 10 dan di sebelah kanan bendungan adalah ℎ0 = 4.

57


Galat solusi perhitungan 𝑢 dan ℎ dihitung dengan menggunakan rumus 𝑁

1 𝐺𝑎𝑙𝑎𝑡 𝑢 = ∑ |𝑢(𝑖) − 𝑢_𝑒𝑥(𝑖)| 𝑁 𝑖=1

dan 𝑁

1 𝐺𝑎𝑙𝑎𝑡 ℎ = ∑ |ℎ(𝑖) − ℎ_𝑒𝑥(𝑖)|, 𝑁 𝑖=1

di mana 𝑢 dan ℎ adalah nilai numeris di titik 𝑥𝑖 , 𝑢_𝑒𝑥 dan ℎ_𝑒𝑥 adalah nilai eksak di titik 𝑥𝑖 , dan 𝑁 adalah banyaknya sel yang digunakan untuk mendiskretkan domain ruang 𝑥.

A.

Metode Volume Hingga Lax-Friedrichs Pada Bab III dibahas mengenai solusi numeris dari persamaan gelombang air

dangkal dengan metode volume hingga. Dalam Bab IV ini dibahas mengenai hasil simulasi solusi numeris persamaan gelombang air dangkal dengan menggunakan MATLAB. Simulasi ini dilakukan untuk beberapa nilai 𝑁, yaitu 100, 200, 400, 800, 1600, dan 3200. Berikut merupakan hasil simulasi numeris untuk persamaan gelombang air dangkal menggunakan metode volume hingga Lax-Friedrichs. Grafik simulasi numeris ditunjukkan pada Gambar 4.2.


Gambar 4.2: Grafik simulasi numeris dengan metode volume hingga LaxFriedrichs untuk 𝑁 = 1000, 𝑡 = 0.2 detik, ∆𝑥 = 0.01, dan ∆𝑡 = 0.05 ∆𝑥. Gambar 4.2 memberikan ilustrasi geometris mengenai simulasi yang telah dilakukan. Terlihat pada gambar bahwa selisih solusi numeris dan eksaknya kecil. Dengan kata lain, solusi numerisnya hampir mendekati solusi eksak. Simulasi ini menggunakan nilai 𝑁 = 1000, ∆𝑥 = 0.01, dan ∆𝑡 = 0.05 ∆𝑥. Simulasi berhenti pada saat 𝑡 = 0.2. Tabel 4.1. Hasil simulasi numeris menggunakan metode volume hingga LaxFriedrichs pada 𝑡 = 0.2

N 100 200 400 800 1600

Galat pada Galat u 0.2650 0.1684 0.1028 0.0601 0.0344

Galat h 0.2112 0.1362 0.0842 0.0498 0.0287


3200 Rata-rata galat

0.0198

0.0166

0.1084

0.0877

Dari Tabel 4.1 tampak bahwa semakin kecil 𝑁 yang diambil, semakin besar galat (error) yang dihasilkan. Ketika diambil nilai 𝑁 yang cukup besar maka galat simulasi akan semakin kecil karena banyaknya partisi pada ruang di sumbu 𝑥 yang semakin banyak dan mendekati solusi eksaknya. Artinya ketika 𝑁 besar solusi yang diperoleh akan lebih akurat. Dapat dilihat juga pada grafik bahwa solusi eksak dan numerisnya berhimpit. Namun pada ujung-ujung grafik, kedua grafik terlihat tidak berhimpit. Artinya masih ada galat perhitungan pada penggunaan metode ini. Ilustrasi galat secara geometris ditunjukkan pada Gambar 4.3.

Gambar 4.3: Grafik galat kecepatan aliran air (𝑢) dan kedalaman air (ℎ) menggunakan metode volume hingga Lax-Friedrichs


Gambar 4.4: Grafik log galat kecepatan aliran air (𝑢) dan log galat kedalaman aliran air (ℎ) menggunakan metode volume hingga Lax-Friedrichs Gambar 4.4 menunjukkan nilai log galat 𝑢 𝐿1, ℎ 𝐿1, 𝑢 𝐿2, ℎ 𝐿2, 𝑢 𝐿∞ dan ℎ 𝐿∞ . Rumus galat 𝐿∞ tidak tepat digunakan sebagai alat untuk mengukur keakuratan metode volume hingga Lax-Friedrichs untuk persamaan gelombang air dangkal. Rumus galat 𝐿∞ tidak dapat menunjukkan bahwa semakin banyak titik, komputasi errornya akan semakin kecil. Sedangkan rumus galat 𝐿1 dan 𝐿2 merupakan rumus yang tepat untuk digunakan sebagai alat untuk mengukur keakuratan metode ini untuk persamaan gelombang air dangkal karena 𝐿1 dan 𝐿2 dapat menunjukkan bahwa semakin banyak titik, komputasi errornya akan semakin kecil.


Kelebihan dari metode ini adalah skema numeriknya sederhana dan perhitungan komputasinya juga sangat mudah. Kekurangan metode ini adalah ketika diambil ∆𝑡 cukup besar metode volume hingga Lax-Friedrichs tidak stabil.

B.

Metode Beda Hingga Grid Kolokasi Selanjutnya akan dipaparkan hasil simulasi numeris penyelesaian gelombang

air dangkal menggunakan metode beda hingga grid koloksi. Pada Bab III telah dibahas solusi numeris untuk model gelombang air dangkal dengan metode beda hingga grid kolokasi. Solusi numeris yang diperoleh disimulasikan dengan menggunakan MATLAB. Kondisi awal yang digunakan adalah sama dengan pada simulasi sebelumnya. Grafik simulasi numeris ditunjukkan pada Gambar 4.5. Grafik yang dihasilkan tidak mulus karena banyak getaran semu pada grafik. Simulasi ini menggunakan nilai 𝑁 = 1000, ∆𝑥 = 0.01, dan ∆𝑡 = 0.05 ∆𝑥. Simulasi berhenti pada saat 𝑡 = 0.2.


Gambar 4.5: Grafik simulasi numeris dengan metode beda hingga grid kolokasi untuk 𝑁 = 1000, 𝑡 = 0.2 detik, ∆𝑥 = 0.01, dan ∆𝑡 = 0.01.

Tampak pada gambar bahwa selisih solusi numeris dengan solusi eksak sangat besar. Hal ini mengakibatkan solusi numeris yang dihasilkan tidak akurat. Berikut merupakan hasil simulasi numeris untuk persamaan gelombang air dangkal menggunakan metode beda hingga grid kolokasi. Tabel 4.2. Hasil simulasi numeris metode beda hingga grid kolokasi pada 𝑡 = 0.2

𝑁 100 200 400 800 1600 3200

Galat pada Galat 𝑢 1.2075 1.2168 1.2122 1.2099 1.2110 1.2105

Galat ℎ 0.9629 0.9751 0.9746 0.9744 0.9759 0.9758


Rata-rata galat

1.2113

0.9731

Dari Tabel 4.2 terlihat bahwa tidak ada perbedaan galat yang signifikan. Misalnya pada 𝑁 = 200 dan 𝑁 = 1600, di mana galat yang dihasilkan hampir sama. Nilai galat 𝑢 cukup besar yaitu berkisar 1.2 sedangkan nilai galat ℎ berkisar 0.9, hal ini tidak sesuai dengan yang diharapkan. Ilustrasi galat secara geometris ditunjukkan pada gambar 4.6.

Gambar 4.6: Grafik galat kecepatan aliran air (𝑢) dan kedalaman air (ℎ) menggunakan metode beda hingga grid kolokasi


Gambar 4.7: Grafik galat kecepatan aliran air (𝑢) dan kedalaman air (ℎ) menggunakan metode beda hingga gris kolokasi Gambar 4.7 menunjukkan nilai log galat 𝑢 𝐿1, ℎ 𝐿1, 𝑢 𝐿2, ℎ 𝐿2, 𝑢 𝐿∞ dan ℎ 𝐿∞ . Rumus galat 𝐿1, 𝐿2, dan 𝐿∞ tidak tepat digunakan sebagai alat untuk mengukur keakuratan metode beda hingga grid kolokasi untuk persamaan gelombang air dangkal. Rumus galat 𝐿1, 𝐿2, dan 𝐿∞ tidak dapat menunjukkan bahwa semakin banyak titik, komputasi errornya akan semakin kecil. Metode ini memiliki kelebihan yaitu ketika diambil ∆𝑡 cukup besar metode ini tetap stabil. Adapun kekurangan dari metode ini yaitu memiliki skema perhitungan numerik yang cukup rumit jika dibandingkan dengan skema eksplisit metode volume hingga Lax-Friedrichs. Metode ini memiliki sifat stabil tanpa syarat. Selain itu, kekurangan pada model ini adalah memiliki galat yang besar pada perhitungan numeriknya karena hanya dipandang satu grid saja untuk menghitung variabel yang tidak diketahui. Dengan demikian metode beda hingga grid kolokasi


bukan pilihan yang tepat untuk menyelesaikan persamaan gelombang air dangkal terkait dengan masalah bendungan bobol.

C.

Metode Beda Hingga Grid Selang-Seling Pada bagian sebelumnya telah dilihat hasil simulasi menggunakan metode

volume hingga dan metode beda hingga grid kolokasi. Pada metode beda hingga grid kolokasi masih terdapat getaran semu dalam hasil grafik yang dapat dilihat pada Gambar 4.8. Sehingga diperlukan model grid selang-seling untuk memperbaiki model grid kolokasi agar tidak ada getaran semu dalam hasil simulasi persamaan gelombang air dangkal. Hasil simulasi numeris menggunakan model grid selang-seling ditunjukkan pada Gambar 4.8.

Gambar 4.8: Grafik simulasi numeris dengan metode beda hingga grid selangseling untuk 𝑁 = 1000, 𝑡 = 0.2 detik, ∆𝑥 = 0.01, dan ∆𝑡 = 0.05 ∆𝑥.


Tampak pada Gambar 4.8 bahwa grafik kedalaman dan kecepatan aliran air terlihat lebih halus dari simulasi-simulasi sebelumnya. Grafik solusi numeris terlihat berhimpit dengan grafik solusi eksak. Namun pada ujung kiri grafik kecepatan aliran air tampak tidak berhimpit. Artinya, galat yang dihasilkan sangat kecil. Dibandingkan dengan kedua metode yang sudah dibahas diatas, metode beda hingga grid selang-seling memberikan hasil yang lebih akurat. Secara grafis terlihat bahwa solusi numeris tampak berhimpit dengan solusi eksaknya. Dengan kata lain, solusi numeris sudah sangat dekat dengan solusi eksaknya. Sekarang akan dilihat hasil simulasi numeris menggunakan metode beda hingga grid selang-seling. Tabel 4.3. Hasil simulasi numeris menggunakan metode beda hingga grid selangseling pada 𝑡 = 0.2 𝑁 100 200 400 800 1600 3200 Rata-rata galat

Galat pada Galat 𝑢 0.1804 0.0897 0.0415 0.0285 0.0144 0.0089

Galat ℎ 0.1239 0.0621 0.0416 0.0209 0.0140 0.0067

0.0605

0.0448

Tabel 4.3 menunjukkan bahwa semakin besar nilai 𝑁 yang diambil, semakin kecil galat yang dihasilkan. Galat yang dihasilkan dengan metode beda hingga selang-seling yaitu setengah kali dari galat dengan nilai 𝑁 yang diambil


sebelumnya. Misalnya pada 𝑁 = 100 diperoleh galat 𝑢 = 0.1804 dan galat ℎ = 0.1239. Kemudian diuji untuk nilai 𝑁 = 200 diperoleh galat 𝑢 = 0.0897 dan galat ℎ = 0.0621. Artinya, galat pada 𝑁 = 200 adalah setengah kalinya galat pada 𝑁 = 100. Ilustrasi galat secara geometris ditunjukkan pada Gambar 4.9.

Gambar 4.9: Grafik galat kecepatan aliran air (𝑢) dan kedalaman air (ℎ) menggunakan metode beda hingga grid selang-seling Galat yang dihasilkan pada metode beda hingga grid selang-seling lebih kecil jika dibandingkan dengan metode beda volume hingga dan metode beda hingga kolokasi. Metode beda hingga grid kolokasi memiliki galat yang lebih besar dari kedua metode yang lain.


Gambar 4.10: Grafik galat kecepatan aliran air (𝑢) dan kedalaman air (ℎ) menggunakan beda hingga grid selang-seling Gambar 4.10 menunjukkan nilai log galat 𝑢 𝐿1, ℎ 𝐿1, 𝑢 𝐿2, ℎ 𝐿2, 𝑢 𝐿∞ dan ℎ 𝐿∞ . Rumus galat 𝐿∞ tidak tepat digunakan sebagai alat untuk mengukur keakuratan metode beda hingga grid selang-seling untuk persamaan gelombang air dangkal. Rumus galat 𝐿∞ tidak dapat menunjukkan bahwa semakin banyak titik, komputasi errornya akan semakin kecil. Sedangkan rumus galat 𝐿1 dan 𝐿2 merupakan rumus yang tepat untuk digunakan sebagai alat untuk mengukur keakuratan metode ini untuk persamaan gelombang air dangkal karena 𝐿1 dan 𝐿2 dapat menunjukkan bahwa semakin banyak titik, komputasi errornya akan semakin kecil. Model grid selang-seling memiliki kelebihan dan kekurangan. Kelebihan dari model ini adalah perhitungan numeriknya lebih akurat dibandingkan dengan


metode beda hingga grid kolokasi dan metode volume hingga Lax-Friedrichs. Kekurangan metode ini yaitu skema untuk menyelesaikan persamaan lebih rumit dibandingkan dengan skema grid kolokasi karena dalam perhitungan dipandang dua grid untuk menghitung setiap variabel yang berbeda. Jadi, setiap variabel yang tidak diketahui dihitung pada grid yang berbeda. Sehingga menghasilkan galat yang kecil. Dengan demikian, untuk menyelesaikan persamaan gelombang air dangkal secara numeris terkait dengan masalah bendungan bobol, penulis menyarankan penggunaan skema beda hingga grid selang-seling. Hasil-hasil pada bab ini telah diseminarkan pada The 2016 International Conference on Information Systems and Applied Mathematics (Mungkasi dan Ilga Purnama Sari, 2016).


BAB V PENUTUP Pada bab ini diberikan kesimpulan atas pembahasan bab-bab sebelumnya serta saran untuk penelitian selanjutnya. A.

Kesimpulan Telah dilihat bahwa untuk menyelesaikan model gelombang air dangkal

secara numeris dapat diselesaikan dengan berbagai metode numeris. Pada bab-bab sebelumnya telah diuji beberapa metode yaitu metode volume hingga, metode beda hingga grid kolokasi dan metode beda hingga grid selang-seling. Berdasarkan hasil simulasi numeris didapatkan penyelesaian yang terbaik dan sesuai yang diharapkan. Nilai rata-rata galat dari ketiga metode yang sudah dibahas pada Bab IV adalah sebagai berikut Metode Numeris Rata-rata galat 𝑢 Metode volume hingga 0.1084 Lax-Friedrichs

Rata-rata galat ℎ 0.0877

Metode beda hingga grid kolokasi

1.2113

0.9731

0.0605

0.0448

Metode beda hingga grid selang-seling

Dari ketiga metode tersebut dapat dilihat bahwa metode beda hingga grid selangseling memiliki rata-rata galat 𝑢 dan ℎ yang paling kecil dibandingkan dengan kedua metode lainnya. Metode beda hingga grid selang-seling adalah metode yang terbaik berdasarkan hasil simulasi numeris yang dilakukan, apabila dibandingkan dengan metode volume hingga dan metode beda hingga grid kolokasi. 71


Dengan demikian, untuk menyelesaikan model gelombang air dangkal terkait dengan masalah bendungan bobol, penulis menyarankan untuk menggunakan metode beda hingga grid selang-seling.

B.

Saran Penulis sadar bahwa dalam penulisan skripsi ini masih banyak kekurangan.

Oleh sebab itu, penulis mengharapkan kelak akan ada yang melanjutkan penelitian ini. Tulisan ini hanya membahas model gelombang air dangkal satu dimensi. Penulis berharap kelak akan ada yang melanjutkan penelitian ini di ruang dimensi yang lebih tinggi dan jika dimungkinkan menggunakan model yang lain yang mungkin memberikan hasil yang lebih baik lagi.


DAFTAR PUSTAKA

Anton, H. (2010). Elementary Linear Algebra: Applications version. Hoboken: John Wiley & Sons. Crowhurst, P. (2013). Numerical Solutions of One-Dimensional Shallow Water Equations. Melbourne: Australian Mathematical Sciences Institute. Crowhurst P. dan Li Z. Proceedings of the 2013 UKSim 15th International Conference on Computer Modelling and Simulation (IEEE, 2013), pp. 55-60. David. B dan George. C. (1995). Basic Partial Differential Equations, Honolulu: Hawaii. Doyen, D. dan Gunawan, P. H. (2014). An Explicit Staggered Finite Volume Scheme for the Shallow Water Equations. Finite Volumes for Complex Apllications VII-Methods and Theoritical Aspects, 227-235. LeVeque, R. J. (1992). Numerical Methods for Conservation Laws, Basel: Birkhauser. LeVeque, R. J. (2004). Finite-Volume Method for Hyperbolic Problems. Cambridge: Cambridge University Press. Mungkasi, S dan Ilga Purnama Sari. (2016). Numerical Solution to the Shallow Water Equations Using Explicit and Implicit Schemes. International Conference on Information System and Applied Mathematics. Submitted to AIP Conference Proceedings. Mungkasi, S. (2008). Finite Volume Methods for the One-Dimensional Shallow Water Equations. Masters Thesis, Canberra: Australian National University. Munir, R. (2008). Metode Numerik, Bandung: Informatika Bandung. Stelling, G. S. dan Duinmeijer, S. P. A. (2003). A Staggered Conservative Scheme for every Froude Number in Rapidly Varied Shallow Water Flows. International Journal for Numerical Methods in Fluids, 43, 1329-1354. Thomas, G. B. (2010). Thomas’ Calculus Early Transcendentals. Boston: Pearson Education.

73


Lampiran Berikut ini adalah code program MATLAB untuk masing-masing metode yang digunakan dalam simulasi numeris untuk persamaan gelombang air dangkal satu dimensi. 1. Metode Volume Hingga Lax-Friedrichs a. Code untuk metode volume hingga Lax-Friedrichs close all; clear all; clc; N = 1000; %banyaknya langkah pada ruang x1 L = 5; dx = 2*L/N; %ukuran langkah pada ruang x x1 = -L:dx:L; %membuat langkah pada ruang dengan dx adalah jarak dari titik a ke b x = x1(1:N)+dx/2; %membuat langkah pada ruang x1 sehingga jarak langkah semakin kecil dt = 0.05*dx; %ukuran langkah waktu t % initial conditions u = zeros(N,1); %membuat ruang untuk kecepatan aliran air u h = zeros(N,1); %membuat ruang untuk kedalaman aliran air h uh = zeros(N,1); %membuat ruang untuk debit air Nm = N/2+1; %kedalaman air di atas kedalaman air disebelah kanan dinding hL=10; hR=4; h(1:Nm) = h(1:Nm)+ hL; %membuat ruang untuk kedalaman air di sebelah kiri dinding h(Nm+1:N) = h(Nm+1:N)+ hR; %membuat ruang untuk kedalaman air di sebelah kanan dinding h1 = h; u1 = u; uh1 = uh; g = 9.81;

%percepatan gravitasi bumi

tStop = 0.2; %waktu akhir iterasi t = 0; %waktu awal j=0;

74


uL=0; uR=0; xmin=-L; xmax=L; delta=(L-(-L))/N; %for shock wave exact computation c_L = sqrt(g*hL); % wave speed at left side c_R = sqrt(g*hR); % wave speed at right side Smin=-101.0; % initial bracketing for bisection method for shock speed S Smax=0.0; % terminal bracketing for bisection method for shock speed S for i=1:100 S=(Smin+Smax)/2.0; % bisecting u2=S-c_R*c_R/4.0/S*(1.0+sqrt(1.0+8.0*S*S/c_R/c_R)); % implicit function for velocity behind shock c2=c_R*sqrt(0.5*(sqrt(1.0+8.0*S*S/c_R/c_R)-1.0)); % implicit function for wave speed behind shock func = -2.0*c_L -u2 +2.0*c2; % characteristic relation between velocity and wave speed if (func > 0.0) % bisection procedure Smin=S; else % bisection procedure Smax=S; end end

while t < tStop j=j+1; for i=2:N-1 F2_right=0.5*(u(i+1)^2*h(i+1) + 0.5*g*h(i+1)^2 + u(i)^2*h(i)+0.5*g*h(i)^2) - dx/(2*dt)*(u(i+1)*h(i+1)-u(i)*h(i)); %menghitung flux 2 kanan F2_left=0.5*(u(i)^2*h(i) + 0.5*g*h(i)^2 + u(i-1)^2*h(i1)+0.5*g*h(i-1)^2) - dx/(2*dt)*(u(i)*h(i)-u(i-1)*h(i-1)); %menghitung flux 2 kiri F1_right=0.5*(u(i+1) *h(i+1) + u(i)*h(i)) dx/(2*dt)*(h(i+1)-h(i)); %menghitung flux 1 kanan F1_left=0.5*(u(i) *h(i) + u(i-1)*h(i-1)) dx/(2*dt)*(h(i)-h(i-1)); %menghitung flux 1 kiri h1(i)=h(i)-dt*(F1_right- F1_left)/dx; %menghitung kedalaman aliran air uh1(i)=uh(i)-dt*(F2_right- F2_left)/dx; %menghitung debit air yang mengalir end t = t + dt;

%update waktu


h = h1; %memperbarui dan memanggil hasil akhir perhitungan kedalaman aliran air uh = uh1; %memperbarui dan memanggil hasil akhir debit air yang mengalir u = uh1./h1; %memperbarui dan memanggil hasil akhir perhitungan kecepatan aliran air %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Exact solution u_ex=xmin+dx/2:dx:xmax-dx/2; % initial storage for exact velocity u_ex h_ex=xmin+dx/2:dx:xmax-dx/2; % initial storage for exact height h_ex % the following iteration computes the exact solution at every spatial point xC=x; for i=1:length(xC) if xC(i) <= -c_L*t u_ex(i) = 0.0; h_ex(i) = hL; elseif xC(i) <= -(u2+c2)*t u_ex(i) = 2.0/3*(xC(i)/t + sqrt(g*hL)); h_ex(i) = 4.0/(9*g)*(-xC(i)/(2.0*t)+sqrt(g*hL))^2; elseif xC(i) < -S*t u_ex(i) = -u2; h_ex(i) = c2^2/g; else u_ex(i) = 0.0; h_ex(i) = hR; end end Q_ex = h_ex.*u_ex; % the exact momentum or discharge Q of water % Galat solusi eksak dengan solusi numeris %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% E_u = (1/N)*sum(abs(u(1:N)'-u_ex(1:N))) E_h = (1/N)*sum(abs(h(1:N)'-h_ex(1:N))) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% subplot(2,1,1) %memesan tempat untuk gambar grafik kedalaman aliran air plot(x,h,'b', xC,h_ex,'r--'); %menggambar kedalaman air legend('Numerical Solution','Exact Solution'); title(sprintf('Depth at t=%4.3f',t)) %membuat judul grafik xlim([-L L]); %membuat ukuran sumbu x ylim([0 11]); %membuat ukuran sumbu y subplot(2,1,2) %memesan tempat untuk gambar grafik debit air yang mengalir plot(x,u,'b', xC,u_ex,'r--'); %menggambar kecepatan air yang mengalir legend('Numerical Solution','Exact Solution'); title(sprintf('Velocity at t=%4.3f',t)) %membuat judul grafik xlim([-L L]); %membuat ukuran sumbu x ylim([0 5]); %membuat ukuran sumbu y


pause(0.001);

%program akan berhenti setiap t detik

end max(h) min(h)

%hasil akhir kedalaman maksimum aliran air %hasil akhir kedalaman minimum aliran air

b. Code untuk grafik galat 𝑢 dan ℎ clc N = [100,200,400,800,1600,3200]; Galat_u = [0.265, 0.1684, 0.1028, 0.0601, 0.0344, 0.0198]; Galat_h = [0,2112, 0,1362, 0,0842, 0,0498, 0,0287, 0,0166]; plot(N,Galat_u,'--r*',N,Galat_h,'--b*') legend('Galat u','Galat h'); xlabel('Banyaknya langkah pada domain ruang (N)') ylabel('Galat') title('Galat Kecepatan dan Kedalaman Aliran Air')

c. Code untuk grafik galat 𝐿1, 𝐿2, dan 𝐿∞ clear clc N1 = [100 200 400 800 1600 3200]; E_u_L1 = zeros(6,1) E_h_L1 = zeros(6,1) E_u_L2 = zeros(6,1) E_h_L2 = zeros(6,1) E_u_Linf = zeros(6,1) E_h_Linf = zeros(6,1) E_u_L1a = E_h_L1a = E_u_L2a = E_h_L2a = E_u_Linfa E_h_Linfa

E_u_L1 E_h_L1 E_u_L2 E_h_L2 = E_u_Linf = E_h_Linf

for kk = 1:length(N1) N = N1(kk) for k=1:length(N) L = 5; dx = 2*L/N; %ukuran langkah pada ruang x x1 = -L:dx:L; %membuat langkah pada ruang dengan dx adalah jarak dari titik a ke b x = x1(1:N)+dx/2; %membuat langkah pada ruang x1 sehingga jarak langkah semakin kecil dt = 0.05*dx; %ukuran langkah waktu t


% initial conditions u = zeros(N,1); %membuat ruang untuk kecepatan aliran air u h

= zeros(N,1);

%membuat ruang untuk kedalaman aliran

air h uh = zeros(N,1);

%membuat ruang untuk debit air

Nm = N/2+1; %kedalaman air di atas kedalaman air disebelah kanan dinding hL=10; hR=4; h(1:Nm) = h(1:Nm)+ hL; %membuat ruang untuk kedalaman air di sebelah kiri dinding h(Nm+1:N) = h(Nm+1:N)+ hR; %membuat ruang untuk kedalaman air di sebelah kanan dinding h1 = h; u1 = u; uh1 = uh; g = 9.81;

%percepatan gravitasi bumi

tStop = 0.2; %waktu akhir iterasi t = 0; %waktu awal j=0; uL=0; uR=0; xmin=-L; xmax=L; delta=(L-(-L))/N; %for shock wave exact computation c_L = sqrt(g*hL); % wave speed at left side c_R = sqrt(g*hR); % wave speed at right side Smin=-101.0; % initial bracketing for bisection method for shock speed S Smax=0.0; % terminal bracketing for bisection method for shock speed S for i=1:100 S=(Smin+Smax)/2.0; % bisecting u2=S-c_R*c_R/4.0/S*(1.0+sqrt(1.0+8.0*S*S/c_R/c_R)); % implicit function for velocity behind shock c2=c_R*sqrt(0.5*(sqrt(1.0+8.0*S*S/c_R/c_R)-1.0)); % implicit function for wave speed behind shock func = -2.0*c_L -u2 +2.0*c2; % characteristic relation between velocity and wave speed if (func > 0.0) % bisection procedure Smin=S; else % bisection procedure Smax=S; end end


while t < tStop j=j+1; for i=2:N-1 F2_right=0.5*(u(i+1)^2*h(i+1) + 0.5*g*h(i+1)^2 + u(i)^2*h(i)+0.5*g*h(i)^2) - dx/(2*dt)*(u(i+1)*h(i+1)-u(i)*h(i)); %menghitung flux 2 kanan F2_left=0.5*(u(i)^2*h(i) + 0.5*g*h(i)^2 + u(i1)^2*h(i-1)+0.5*g*h(i-1)^2) - dx/(2*dt)*(u(i)*h(i)-u(i-1)*h(i-1)); %menghitung flux 2 kiri F1_right=0.5*(u(i+1) *h(i+1) + u(i)*h(i)) dx/(2*dt)*(h(i+1)-h(i)); %menghitung flux 1 kanan F1_left=0.5*(u(i) *h(i) + u(i-1)*h(i-1)) dx/(2*dt)*(h(i)-h(i-1)); %menghitung flux 1 kiri h1(i)=h(i)-dt*(F1_right- F1_left)/dx; %menghitung kedalaman aliran air uh1(i)=uh(i)-dt*(F2_right- F2_left)/dx; %menghitung debit air yang mengalir end t = t + dt; %update waktu h = h1; %memperbarui dan memanggil hasil akhir perhitungan kedalaman aliran air uh = uh1; %memperbarui dan memanggil hasil akhir debit air yang mengalir u = uh1./h1; %memperbarui dan memanggil hasil akhir perhitungan kecepatan aliran air %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Exact solution u_ex=xmin+dx/2:dx:xmax-dx/2; % initial storage for exact velocity u_ex h_ex=xmin+dx/2:dx:xmax-dx/2; % initial storage for exact height h_ex % the following iteration computes the exact solution at every spatial point xC=x; for i=1:length(xC) if xC(i) <= -c_L*t u_ex(i) = 0.0; h_ex(i) = hL; elseif xC(i) <= -(u2+c2)*t u_ex(i) = 2.0/3*(xC(i)/t + sqrt(g*hL)); h_ex(i) = 4.0/(9*g)*(xC(i)/(2.0*t)+sqrt(g*hL))^2; elseif xC(i) < -S*t u_ex(i) = -u2; h_ex(i) = c2^2/g; else u_ex(i) = 0.0; h_ex(i) = hR;


end end Q_ex = h_ex.*u_ex; % the exact momentum or discharge Q of water subplot(2,1,1) %memesan tempat kedalaman aliran air plot(x,h,'b', xC,h_ex,'r--'); air legend('Solusi Numerik','Solusi title(sprintf('Kedalaman air at judul grafik xlim([-L L]); %membuat ukuran ylim([0 11]); %membuat ukuran

untuk gambar grafik %menggambar kedalaman Eksak'); t=%4.3f',t))

%membuat

sumbu x sumbu y

subplot(2,1,2) %memesan tempat untuk gambar grafik debit air yang mengalir plot(x,u,'b', xC,u_ex,'r--'); %menggambar kecepatan air yang mengalir legend('Solusi Numerik','Solusi Eksak'); title(sprintf('Kecepatan aliran air at t=%4.3f',t)) %membuat judul grafik xlim([-L L]); %membuat ukuran sumbu x ylim([0 5]); %membuat ukuran sumbu y pause(0.001); %program akan berhenti setiap t detik end end

% Galat L1, L2, dan L inf E_u_L1a(kk,1) =(1/N)*sum(abs(u(1:N)'-u_ex(1:N))) E_h_L1a(kk,1) =(1/N)*sum(abs(h(1:N)'-h_ex(1:N))) E_u_L2a(kk,1) = (1/N)*sum((u(1:N)'-u_ex(1:N)).^2) E_h_L2a(kk,1) = (1/N)*sum((h(1:N)'-h_ex(1:N)).^2) E_u_Linfa(kk,1) = E_h_Linfa(kk,1) =

max(abs(u(1:N)'-u_ex(1:N))) max(abs(h(1:N)'-h_ex(1:N)))

end E_u_L1 = E_h_L1 = E_u_L2 = E_h_L2 = E_u_Linf E_h_Linf

E_u_L1a E_h_L1a E_u_L2a E_h_L2a = E_u_Linfa = E_h_Linfa

subplot(2,1,1) plot(log(N1),log(E_u_L1),'--r*',log(N1),log(E_u_L2),'--b*', log(N1), log(E_u_Linf),'--g*') legend('Galat u L1','Galat u L2', 'Galat u Linf'); xlabel('log(N)') ylabel('log(Galat u)') title('Galat Kecepatan Aliran Air')


subplot(2,1,2) plot(log(N1),log(E_h_L1),'--r*',log(N1),log(E_h_L2),'--b*', log(N1), log(E_h_Linf),'--g*') legend('Galat h L1','Galat h L2', 'Galat h Linf'); xlabel('log(N)') ylabel('log(Galat h)') title('Galat Kedalaman Aliran Air')

2. Metode Beda Hingga Grid Kolokasi a. Code untuk metode beda hingga grid kolokasi clear all; clear figure; N=1000; % jumlah langkah pada ruang zeta=0.01; % ukuran langkah waktu tStop=0.2; M=tStop/zeta; % jumlah langkah pada waktu L=5; hL=10; hR=4; uL=0; uR=0; xmin=-L; xmax=L; g=9.81; % konstanta percepatan gravitasi delta=(L-(-L))/N; dx=delta; u(:,:)=zeros(N+1,M+1); % membuat ruang untuk kecepatan h(:,:)=zeros(N+1,M+1); % membuat ruang untuk ketinggian u(1,:)=uL; u(N+1,:)=uR; h(1,:)=hL; h(N+1,:)=hR; x=-L:delta:L; % langkah ruang dan jarak t = 0; %for shock wave exact computation c_L = sqrt(g*hL); % wave speed at left side c_R = sqrt(g*hR); % wave speed at right side Smin=-101.0; % initial bracketing for bisection method for shock speed S Smax=0.0; % terminal bracketing for bisection method for shock speed S for i=1:100 S=(Smin+Smax)/2.0; % bisecting u2=S-c_R*c_R/4.0/S*(1.0+sqrt(1.0+8.0*S*S/c_R/c_R)); % implicit function for velocity behind shock c2=c_R*sqrt(0.5*(sqrt(1.0+8.0*S*S/c_R/c_R)-1.0)); % implicit function for wave speed behind shock


func = -2.0*c_L -u2 characteristic relation if (func > 0.0) % Smin=S; else % Smax=S; end end

+2.0*c2; % between velocity and wave speed bisection procedure bisection procedure

for k=2:N if k < N/2 h(k,1)= hL; % kedalaman aliran air di kiri dinding else h(k,1)= hR; % kedalaman aliran air di kanan dinding end end % Penyelesaian Matriks A=zeros(2*(N-1),2*(N-1)); %memesan tempat untuk matriks A b=zeros(2*(N-1),1); %memesan tempat untuk matriks b for j=2:M+1 A(1,N)=2*delta; A(1,N+1)=zeta*u(2,j-1); A(1,2)=zeta*h(2,j-1); % Persamaan pertama untuk i=1 b(1,1)=2*delta*h(2,j-1)+zeta*u(2,j-1)*h(1,j)+zeta*h(2,j1)*u(1,j); A(2,1)=2*delta*h(2,j-1); A(2,2)=zeta*u(2,j-1)*h(2,j-1); % persamaan kedua untuk i=1 A(2,N+1)=zeta*g*h(2,j-1); b(2,1)=2*delta*h(2,j-1)*u(2,j-1)+zeta*u(2,j-1)*h(2,j1)*u(1,j)+zeta*g*h(2,j-1)*h(1,j); A(2*N-3,2*N-2)=2*delta; A(2*N-3,N-2+N-1)=-zeta*u(N,j-1); A(2*N-3,N-1-1)=-zeta*h(N,j-1); % Persamaan pertama untuk i=3 b(2*N-3,1)=2*delta*h(N,j-1)-zeta*u(N,j-1)*h(N+1,j)-zeta*h(N,j1)*u(N+1,j); A(2*N-2,N-1)=2*delta*h(N,j-1); A(2*N-2,N-1-1)=-zeta*u(N,j-1)*h(N,j-1); A(2*N-2,N-2+N-1)=-zeta*g*h(N,j-1);

% Persamaan kedua untuk i=3 b(2*N-2,1)=2*delta*h(N,j-1)*u(N,j-1)-zeta*u(N,j-1)*h(N,j1)*u(N+1,j)-zeta*g*h(N,j-1)*h(N+1,j); for i=3:N-1 A(2*i-3,N-2+i)=2*delta; A(2*i-3,N-2+i+1)=zeta*u(i,j-1);


A(2*i-3,N-2+i-1)=-zeta*u(i,j-1); A(2*i-3,i+1-1)=zeta*h(i,j-1); A(2*i-3,i-1-1)=-zeta*h(i,j-1); % persamaan pertama untuk i=2 b(2*i-3,1)=2*delta*h(i,j-1); A(2*i-2,i-1)=2*delta*h(i,j-1); A(2*i-2,i+1-1)=zeta*u(i,j-1)*h(i,j-1); A(2*i-2,i-1-1)=-zeta*u(i,j-1)*h(i,j-1); A(2*i-2,N-2+i+1)=zeta*g*h(i,j-1); A(2*i-2,N-2+i-1)=-zeta*g*h(i,j-1); % persamaan kedua untuk i=2 b(2*i-2,1)=2*delta*h(i,j-1)*u(i,j-1); end % menyelesaikan matriks A y=A\b;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % SOLUSI NUMERIS % solusi untuk ruang kecepatan dan kedalaman aliran air u(2:N,j)=y(1:N-1); h(2:N,j)=y(N:2*N-2); % Syarat batas u(1:5,j)=0; u(end-5:end,j)=0; h(1:5,j)=10; h(end-5:end,j)=4; t = t+zeta; %%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % SOLUSI EKSAK u_ex=xmin+dx/2:dx:xmax-dx/2; % initial storage for exact velocity u_ex h_ex=xmin+dx/2:dx:xmax-dx/2; % initial storage for exact height h_ex % the following iteration computes the exact solution at every spatial point xC=x; for i=1:length(xC) if xC(i) <= -c_L*t u_ex(i) = 0.0; h_ex(i) = hL; elseif xC(i) <= -(u2+c2)*t u_ex(i) = 2.0/3*(xC(i)/t + sqrt(g*hL)); h_ex(i) = 4.0/(9*g)*(-xC(i)/(2.0*t)+sqrt(g*hL))^2; elseif xC(i) < -S*t u_ex(i) = -u2; h_ex(i) = c2^2/g; else u_ex(i) = 0.0;


h_ex(i) = hR; end end Q_ex = h_ex.*u_ex; % the exact momentum or discharge Q of water

% Galat solusi eksak dengan solusi numeris %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% E_u = (1/N)*sum(abs(u(1:N)-u_ex(1:N))) E_h = (1/N)*sum(abs(h(1:N)-h_ex(1:N))) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % membuat gambar grafik kedalaman dan kecepatan air subplot(2,1,1); plot(x,h(:,j),'b-', xC,h_ex, 'r--'); legend('Numerical Solution','Exact Solution'); title(sprintf('Depth at t=%4.3f',t)) axis([-L L 0 11]); subplot(2,1,2); plot(x,u(:,j),'b-',xC,u_ex,'r--'); legend('Numerical Solution','Exact Solution'); title(sprintf('Velocity at t=%4.3f',t)) axis([-L L 0 5]); pause(0.1) end

b. Code untuk grafik galat 𝑢 dan ℎ clc N = [100,200,400,800,1600,3200]; Galat_u = [1.2075, 1.2168, 1.2122, 1.2099, 1.2110, 1.2105]; Galat_h = [0.9629, 0.9751, 0.9746, 0.9744, 0.9759, 0.9758]; plot(N,Galat_u,'--r*',N,Galat_h,'--b*') legend('Galat u','Galat h'); xlabel('Banyaknya langkah pada domain ruang (N)') ylabel('Galat') title('Galat Kecepatan dan Kedalaman Aliran Air')

c. Code untuk grafik galat 𝐿1, 𝐿2, dan 𝐿∞ clear all; clear figure; N1 = [100 200 400 800 1600 3200]; E_u_L1=zeros(6,1) E_h_L1=zeros(6,1) E_u_L2=zeros(6,1) E_h_L2=zeros(6,1)


E_u_Linf=zeros(6,1) E_h_Linf=zeros(6,1) E_u_L1a=E_u_L1 E_h_L1a=E_h_L1 E_u_L2a=E_u_L2 E_h_L2a=E_h_L2 E_u_Linfa=E_u_Linf E_h_Linfa=E_h_Linf for kk = 1:length(N1) N = N1(kk) for jj=1:length(N) zeta=0.01; % ukuran langkah waktu tStop=0.2; M=tStop/zeta; % jumlah langkah pada waktu L=5; hL=10; hR=4; uL=0; uR=0; xmin=-L; xmax=L; g=9.81; % konstanta percepatan gravitasi delta=(L-(-L))/N; dx=delta; u=zeros(N+1,M+1); % membuat ruang untuk kecepatan h=zeros(N+1,M+1); % membuat ruang untuk kedalaman u(1,:)=uL; u(N+1,:)=uR; h(1,:)=hL; h(N+1,:)=hR; x=-L:delta:L; % langkah ruang dan jarak t = 0; %for shock wave exact computation c_L = sqrt(g*hL); % wave speed at left side c_R = sqrt(g*hR); % wave speed at right side Smin=-101.0; % initial bracketing for bisection method for shock speed S Smax=0.0; % terminal bracketing for bisection method for shock speed S for i=1:100 S=(Smin+Smax)/2.0; % bisecting u2=S-c_R*c_R/4.0/S*(1.0+sqrt(1.0+8.0*S*S/c_R/c_R)); % implicit function for velocity behind shock c2=c_R*sqrt(0.5*(sqrt(1.0+8.0*S*S/c_R/c_R)-1.0)); % implicit function for wave speed behind shock func = -2.0*c_L -u2 +2.0*c2; % characteristic relation between velocity and wave speed if (func > 0.0) % bisection procedure Smin=S; else % bisection procedure Smax=S; end end


for k=2:N if k < N/2 h(k,1)= hL; % kedalaman aliran air di kiri dinding else h(k,1)= hR; % kedalaman aliran air di kanan dinding end end % Penyelesaian Matriks A=zeros(2*(N-1),2*(N-1)); %memesan tempat untuk matriks A b=zeros(2*(N-1),1);

%memesan tempat untuk matriks b

for j=2:M+1 A(1,N)=2*delta; A(1,N+1)=zeta*u(2,j-1); A(1,2)=zeta*h(2,j-1); % Persamaan pertama untuk i=1 b(1,1)=2*delta*h(2,j-1)+zeta*u(2,j1)*h(1,j)+zeta*h(2,j-1)*u(1,j); A(2,1)=2*delta*h(2,j-1); A(2,2)=zeta*u(2,j-1)*h(2,j-1); % persamaan kedua untuk i=1 A(2,N+1)=zeta*g*h(2,j-1); b(2,1)=2*delta*h(2,j-1)*u(2,j-1)+zeta*u(2,j-1)*h(2,j1)*u(1,j)+zeta*g*h(2,j-1)*h(1,j); A(2*N-3,2*N-2)=2*delta; A(2*N-3,N-2+N-1)=-zeta*u(N,j-1); A(2*N-3,N-1-1)=-zeta*h(N,j-1); % Persamaan pertama untuk i=3 b(2*N-3,1)=2*delta*h(N,j-1)-zeta*u(N,j-1)*h(N+1,j)zeta*h(N,j-1)*u(N+1,j); A(2*N-2,N-1)=2*delta*h(N,j-1); A(2*N-2,N-1-1)=-zeta*u(N,j-1)*h(N,j-1); A(2*N-2,N-2+N-1)=-zeta*g*h(N,j-1);

% Persamaan kedua untuk i=3 b(2*N-2,1)=2*delta*h(N,j-1)*u(N,j-1)-zeta*u(N,j1)*h(N,j-1)*u(N+1,j)-zeta*g*h(N,j-1)*h(N+1,j); for i=3:N-1 A(2*i-3,N-2+i)=2*delta; A(2*i-3,N-2+i+1)=zeta*u(i,j-1); A(2*i-3,N-2+i-1)=-zeta*u(i,j-1); A(2*i-3,i+1-1)=zeta*h(i,j-1); A(2*i-3,i-1-1)=-zeta*h(i,j-1); % persamaan pertama untuk i=2


b(2*i-3,1)=2*delta*h(i,j-1); A(2*i-2,i-1)=2*delta*h(i,j-1); A(2*i-2,i+1-1)=zeta*u(i,j-1)*h(i,j-1); A(2*i-2,i-1-1)=-zeta*u(i,j-1)*h(i,j-1); A(2*i-2,N-2+i+1)=zeta*g*h(i,j-1); A(2*i-2,N-2+i-1)=-zeta*g*h(i,j-1); % persamaan kedua untuk i=2 b(2*i-2,1)=2*delta*h(i,j-1)*u(i,j-1); end % menyelesaikan matriks A y=A\b;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % SOLUSI NUMERIS % solusi untuk ruang kecepatan dan kedalaman aliran air u(2:N,j)=y(1:N-1); h(2:N,j)=y(N:2*N-2); % Syarat batas u(1:5,j)=0; u(end-5:end,j)=0; h(1:5,j)=10; h(end-5:end,j)=4; t = t+zeta; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % SOLUSI EKSAK u_ex=xmin+dx/2:dx:xmax-dx/2; % initial storage for exact velocity u_ex h_ex=xmin+dx/2:dx:xmax-dx/2; % initial storage for exact height h_ex % the following iteration computes the exact solution at every spatial point xC=x; for i=1:length(xC) if xC(i) <= -c_L*t u_ex(i) = 0.0; h_ex(i) = hL; elseif xC(i) <= -(u2+c2)*t u_ex(i) = 2.0/3*(xC(i)/t + sqrt(g*hL)); h_ex(i) = 4.0/(9*g)*(xC(i)/(2.0*t)+sqrt(g*hL))^2; elseif xC(i) < -S*t u_ex(i) = -u2; h_ex(i) = c2^2/g; else u_ex(i) = 0.0; h_ex(i) = hR;


end end Q_ex = h_ex.*u_ex; % the exact momentum or discharge Q of water

% Galat solusi eksak dengan solusi numeris %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% E_u = (1/N)*sum(abs(u(1:N)-u_ex(1:N))) E_h = (1/N)*sum(abs(h(1:N)-h_ex(1:N))) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% membuat gambar grafik kedalaman dan kecepatan air subplot(2,1,1); plot(x,h(:,j),'b-', xC,h_ex, 'r--'); legend('Numerical Solution','Exact Solution'); title(sprintf('Depth at t=%4.3f',t)) axis([-L L 0 11]); subplot(2,1,2); plot(x,u(:,j),'b-',xC,u_ex,'r--'); legend('Numerical Solution','Exact Solution'); title(sprintf('Velocity at t=%4.3f',t)) axis([-L L 0 5]); pause(0.1) end end % Galat L1, L2, dan L inf E_u_L1a(kk,1) =(1/N)*sum(abs(u(1:N)-u_ex(1:N))) E_h_L1a(kk,1) =(1/N)*sum(abs(h(1:N)-h_ex(1:N))) E_u_L2a(kk,1) = (1/N)*sum((u(1:N)-u_ex(1:N)).^2) E_h_L2a(kk,1) = (1/N)*sum((h(1:N)-h_ex(1:N)).^2) E_u_Linfa(kk,1) = E_h_Linfa(kk,1) = end E_u_L1 = E_h_L1 = E_u_L2 = E_h_L2 = E_u_Linf E_h_Linf

max(abs(u(1:N)-u_ex(1:N))) max(abs(h(1:N)-h_ex(1:N)))


subplot(2,1,1) plot(log(N1),log(E_u_L1),'--r*',log(N1),log(E_u_L2),'--b*', log(N1), log(E_u_Linf),'--g*')


legend('Galat u L1','Galat u L2', 'Galat u Linf'); xlabel('log(N)') ylabel('log(Galat u)') title('Galat Kecepatan Aliran Air') subplot(2,1,2) plot(log(N1),log(E_h_L1),'--r*',log(N1),log(E_h_L2),'--b*', log(N1), log(E_h_Linf),'--g*') legend('Galat h L1','Galat h L2', 'Galat h Linf'); xlabel('log(N)') ylabel('log(Galat h)') title('Galat Kedalaman Aliran Air')

3. Metode Beda Hinggs Grid Selang-Seling a. Code untuk metode beda hingga grid selang-seling clc clear all g = 9.81; % percepatan gravitasi N = 100; % banyaknya langkah untuk mendiskritkan domain ruang L = 5; % ujung-ujung interval xmin=-L; xmax=L; dx = (L-(-L))/N; % jarak grid a ke grid b hL = 10; % ketinggian awal di sebelah kiri dinding, dinding berada pada x=0 hR = 4; % ketinggian awal di sebelah kanan dinding dt = 0.05*dx; % langkah waktu t tstop =0.2; % program ini berjalan sampai waktu tStop detik %for shock wave exact computation c_L = sqrt(g*hL); % wave speed at left side c_R = sqrt(g*hR); % wave speed at right side Smin=-101.0; % initial bracketing for bisection method for shock speed S Smax=0.0; % terminal bracketing for bisection method for shock speed S for i=1:100 S=(Smin+Smax)/2.0; % bisecting u2=S-c_R*c_R/4.0/S*(1.0+sqrt(1.0+8.0*S*S/c_R/c_R)); % implicit function for velocity behind shock c2=c_R*sqrt(0.5*(sqrt(1.0+8.0*S*S/c_R/c_R)-1.0)); % implicit function for wave speed behind shock func = -2.0*c_L -u2 +2.0*c2; % characteristic relation between velocity and wave speed if (func > 0.0) % bisection procedure Smin=S; else % bisection procedure Smax=S; end end


% make a grid x = (-L+dx/2 : dx : L+dx/2); % membuat langkah pada ruang x dengan ujung-ujung interval -L+dx/2 dan L+dx/2 h = zeros(1,N); % membuat ruang untuk kedalaman aliran air u = zeros(1,N); % membuat ruang untuk kecepatan aliran air q = zeros(1,N); % membuat ruang untuk debit aliran air xodd = zeros(1,N/2); % membuat ruang untuk menghitung kecepatan aliran air uodd = zeros(1,N/2); % membuat ruang untuk perhitungan kecepatan aliran air pada grid ganjil qodd = zeros(1,N/2); % membuat ruang untuk perhitungan debit aliran air pada grid ganjil xeven = zeros(1,N/2); % membuat ruang untuk menghitung kedalaman aliran air heven = zeros(1,N/2); % membuat ruang untuk perhitungan kedalaman aliran air pada grid genap % menghitung kecepatan dan kedalaman aliran air. Jika hasil modulonya nol % maka hasil perhitungannya masuk ke grid genap. Jika hasil modulonya tidak nol % maka hasil perhitungannya masuk ke grid ganjil. for i = 1:N if mod(i,2) == 0 xeven(i/2) = x(i); else xodd((i+1)/2) = x(i); end end %initial condition h(1,1:N/2) = hL; h(1,N/2:N) = hR; %now let start the evolution figure(1) t=0; %initial condition h(1) = hL; h(end) = hR; u(1) = 0; u(end) = 0; while (t

hhtimh = hh(i-2); uuimh = uu(i-1); qimh = hhtimh*uuimh; hnew(i) = hh(i) - dt/(2*dx)*(qiph - qimh); end end %boundary condition hnew(1:3) = hL; hnew(end-3:end) = hR; unew(1:3) = 0; unew(end-3:end) = 0; % perhitungan dengan metode beda hingga grid selang-seling for i=4:N-4 if mod(i,2)~=0; %u grids i.e. odd grids uutip1 = uu(i); uuti = uu(i-2); uuip3h = uu(i+2); hhtip3h = hh(i+1); qip3h = hhtip3h*uuip3h; uuiph = uu(i); hhtiph = hh(i-1); qiph = hhtiph*uuiph; qip1 = 0.5*(qiph+qip3h); uuimh = uu(i-2); hhtimh = hh(i-3); qimh = hhtimh*uuimh; qi = 0.5*(qimh+qiph); hhip1 = hh(i+1); hhi = hh(i-1); hhiph = 0.5*(hhi+hhip1); hip1 = hnew(i+1); hi = hnew(i-1); hiph = 0.5*(hi+hip1); unew(i) = (1/hiph)*(hhiph*uuiph dt/(2*dx)*(qip1*uutip1 - qi*uuti + 0.5*g*(hip1^2 - hi^2) ) ); end end %boundary condition hnew(1:4) = hL; hnew(end-4:end) = hR; unew(1:4) = 0; unew(end-4:end) = 0; h = hnew; u = unew;

% memperbarui hasil perhitungan h % memperbarui hasil perhitungan u


t=t+dt; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Exact solution u_ex=xmin+dx/2:dx:xmax-dx/2; % initial storage for exact velocity u_ex h_ex=xmin+dx/2:dx:xmax-dx/2; % initial storage for exact height h_ex % the following iteration computes the exact solution at every spatial point xC=x; for i=1:length(xC) if xC(i) <= -c_L*t u_ex(i) = 0.0; h_ex(i) = hL; elseif xC(i) <= -(u2+c2)*t u_ex(i) = 2.0/3*(xC(i)/t + sqrt(g*hL)); h_ex(i) = 4.0/(9*g)*(-xC(i)/(2.0*t)+sqrt(g*hL))^2; elseif xC(i) < -S*t u_ex(i) = -u2; h_ex(i) = c2^2/g; else u_ex(i) = 0.0; h_ex(i) = hR; end end Q_ex = h_ex.*u_ex; % the exact momentum or discharge Q of water % Menggambar hasil perhitungan h dan u hevenex = 0*heven; uoddex = 0*uodd; for i = 1:N if mod(i,2) == 0 heven(i/2) = h(i); hevenex(i/2) = h_ex(i); else uodd((i+1)/2) = u(i); uoddex((i+1)/2) = u_ex(i); end end % Galat solusi eksak dengan solusi numeris %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% E_u = 1/(N/2)*sum(abs(uodd-uoddex)) E_h = 1/(N/2)*sum(abs(heven-hevenex)) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure(1) subplot(2,1,1) plot(xeven,heven,'b-', xeven,hevenex,'r-') legend('solusi numeris','solusi eksak')


title(sprintf('Kedalaman air pada waktu t=%4.3f',t)) xlim([-L L]) ylim([0 11]) subplot(2,1,2) plot(xodd,uodd,'b-', xodd,uoddex,'r-') legend('solusi numeris','solusi eksak') title(sprintf('Kecepatan air pada waktu t=%4.3f',t)) xlim([-L L]) ylim([0 5]) pause(0.1) end

b. Code untuk grafik galat 𝑢 dan ℎ clc N = [100,200,400,800,1600,3200]; Galat_u = [0.1804, 0.0897, 0.0415, 0.0285, 0.0144, 0.0089]; Galat_h = [0.1239, 0.0621, 0.0416, 0.0209, 0.0140, 0.0067]; plot(N,Galat_u,'--r*',N,Galat_h,'--b*') legend('Galat u','Galat h'); xlabel('Banyaknya langkah pada domain ruang (N)') ylabel('Galat') title('Galat Kecepatan dan Kedalaman Aliran Air')

c. Code untuk grafik galat 𝐿1, 𝐿2, dan 𝐿∞ clc clear all N1 = [100 200 400 800 1600 3200]; E_u_L1=zeros(6,1) E_h_L1=zeros(6,1) E_u_L2=zeros(6,1) E_h_L2=zeros(6,1) E_u_Linf=zeros(6,1) E_h_Linf=zeros(6,1) E_u_L1a=E_u_L1 E_h_L1a=E_h_L1 E_u_L2a=E_u_L2 E_h_L2a=E_h_L2 E_u_Linfa=E_u_Linf E_h_Linfa=E_h_Linf for kk = 1:length(N1) N = N1(kk) for k=1:length(N) g = 9.81; % percepatan gravitasi L = 5; % ujung-ujung interval


xmin=-L; xmax=L; dx = (L-(-L))/N; % jarak grid a ke grid b hL = 10; % ketinggian awal di sebelah kiri dinding, dinding berada pada x=0 hR = 4; % ketinggian awal di sebelah kanan dinding dt = 0.05*dx; % langkah waktu t tstop =0.2; % program ini berjalan sampai waktu tStop detik %for shock wave exact computation c_L = sqrt(g*hL); % wave speed at left side c_R = sqrt(g*hR); % wave speed at right side Smin=-101.0; % initial bracketing for bisection method for shock speed S Smax=0.0; % terminal bracketing for bisection method for shock speed S for i=1:100 S=(Smin+Smax)/2.0; % bisecting u2=S-c_R*c_R/4.0/S*(1.0+sqrt(1.0+8.0*S*S/c_R/c_R)); % implicit function for velocity behind shock c2=c_R*sqrt(0.5*(sqrt(1.0+8.0*S*S/c_R/c_R)-1.0)); % implicit function for wave speed behind shock func = -2.0*c_L -u2 +2.0*c2; % characteristic relation between velocity and wave speed if (func > 0.0) % bisection procedure Smin=S; else % bisection procedure Smax=S; end end % make a grid x = (-L+dx/2 : dx : L+dx/2); % membuat langkah pada ruang x dengan ujung-ujung interval -L+dx/2 dan L+dx/2 h = zeros(1,N); % membuat ruang untuk kedalaman aliran air u = zeros(1,N); % membuat ruang untuk kecepatan aliran air q = zeros(1,N); % membuat ruang untuk debit aliran air xodd = zeros(1,N/2); % membuat ruang untuk menghitung kecepatan aliran air uodd = zeros(1,N/2); % membuat ruang untuk perhitungan kecepatan aliran air pada grid ganjil qodd = zeros(1,N/2); % membuat ruang untuk perhitungan debit aliran air pada grid ganjil xeven = zeros(1,N/2); % membuat ruang untuk menghitung kedalaman aliran air heven = zeros(1,N/2); % membuat ruang untuk perhitungan kedalaman aliran air pada grid genap % menghitung kecepatan dan kedalaman aliran air. Jika hasil modulonya nol % maka hasil perhitungannya masuk ke grid genap. Jika hasil modulonya tidak nol % maka hasil perhitungannya masuk ke grid ganjil. for i = 1:N if mod(i,2) == 0


xeven(i/2) = x(i); else xodd((i+1)/2) = x(i); end end %initial condition h(1,1:N/2) = hL; h(1,N/2:N) = hR; %now let start the evolution figure(1) t=0; %initial condition h(1) = hL; h(end) = hR; u(1) = 0; u(end) = 0; while (t

hhtiph = hh(i-1); qiph = hhtiph*uuiph; qip1 = 0.5*(qiph+qip3h); uuimh = uu(i-2); hhtimh = hh(i-3); qimh = hhtimh*uuimh; qi = 0.5*(qimh+qiph); hhip1 = hh(i+1); hhi = hh(i-1); hhiph = 0.5*(hhi+hhip1); hip1 = hnew(i+1); hi = hnew(i-1); hiph = 0.5*(hi+hip1); unew(i) = (1/hiph)*(hhiph*uuiph dt/(2*dx)*(qip1*uutip1 - qi*uuti + 0.5*g*(hip1^2 - hi^2) ) ); end end %boundary condition hnew(1:4) = hL; hnew(end-4:end) = hR; unew(1:4) = 0; unew(end-4:end) = 0; h = hnew; u = unew;

% memperbarui hasil perhitungan h % memperbarui hasil perhitungan u

t=t+dt; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Exact solution u_ex=xmin+dx/2:dx:xmax-dx/2; % initial storage for exact velocity u_ex h_ex=xmin+dx/2:dx:xmax-dx/2; % initial storage for exact height h_ex % the following iteration computes the exact solution at every spatial point xC=x; for i=1:length(xC) if xC(i) <= -c_L*t u_ex(i) = 0.0; h_ex(i) = hL; elseif xC(i) <= -(u2+c2)*t u_ex(i) = 2.0/3*(xC(i)/t + sqrt(g*hL)); h_ex(i) = 4.0/(9*g)*(xC(i)/(2.0*t)+sqrt(g*hL))^2; elseif xC(i) < -S*t u_ex(i) = -u2; h_ex(i) = c2^2/g; else u_ex(i) = 0.0; h_ex(i) = hR;


end end Q_ex = h_ex.*u_ex; % the exact momentum or discharge Q of water % Menggambar hasil perhitungan h dan u hevenex = 0*heven; uoddex = 0*uodd; for i = 1:N if mod(i,2) == 0 heven(i/2) = h(i); hevenex(i/2) = h_ex(i); else uodd((i+1)/2) = u(i); uoddex((i+1)/2) = u_ex(i); end end % Galat solusi eksak dengan solusi numeris %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% E_u = 1/(N/2)*sum(abs(uodd-uoddex)) E_h = 1/(N/2)*sum(abs(heven-hevenex)) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure(1) subplot(2,1,1) plot(xeven,heven,'b-', xeven,hevenex,'r-') legend('Numerical solution','Exact Solution') title(sprintf('Depth at t=%4.3f',t)) xlim([-L L]) ylim([0 11]) subplot(2,1,2) plot(xodd,uodd,'b-', xodd,uoddex,'r-') legend('Numerical solution','Exact Solution') title(sprintf('Velocity at t=%4.3f',t)) xlim([-L L]) ylim([0 5]) pause(0.1) end end % Galat solusi eksak dengan solusi numeris E_u_L1a(kk,1) =(1/N)*sum(abs(uodd-uoddex)) E_h_L1a(kk,1) =(1/N)*sum(abs(heven-hevenex)) E_u_L2a(kk,1) = (1/N)*sum((uodd-uoddex).^2) E_h_L2a(kk,1) = (1/N)*sum((heven-hevenex).^2) E_u_Linfa(kk,1) = E_h_Linfa(kk,1) = end

max(abs(uodd-uoddex)) max(abs(heven-hevenex))


E_u_L1 = E_h_L1 = E_u_L2 = E_h_L2 = E_u_Linf E_h_Linf


subplot(2,1,1) plot(log(N1),log(E_u_L1),'--r*',log(N1),log(E_u_L2),'--b*', log(N1), log(E_u_Linf),'--g*') legend('Galat u L1','Galat u L2', 'Galat u Linf'); xlabel('log(N)') ylabel('log(Galat u)') title('Galat Kecepatan Aliran Air') subplot(2,1,2) plot(log(N1),log(E_h_L1),'--r*',log(N1),log(E_h_L2),'--b*', log(N1), log(E_h_Linf),'--g*') legend('Galat h L1','Galat h L2', 'Galat h Linf'); xlabel('log(N)') ylabel('log(Galat h)') title('Galat Kedalaman Aliran Air')

PENYELESAIAN PERSAMAAN GELOMBANG AIR DANGKAL DENGAN BEBERAPA METODE NUMERIS

Recommend Documents