1 Pendahuluan pdp 2. 4 Persamaan Difusi Prinsip Maksimum Fungsi Green Metoda separasi variable, recall

1 Pendahuluan pdp 2 Persamaan Type Hiperbolik 2.1 Persamaan Transport . . . . 2.1.1 Metoda karakteristik . 2.1.2 Koefisien tak konstan 2.2 Persamaan Gelombang . . . . 2.2.1 Energi . . . . . . . . . 2.3 Well-posed problem . . . . .

SRP '12

Contents 2

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

3 Metoda Beda Hingga pada Persamaan Tipe Hiperbolik 3.1 Persamaan Transport . . . . . . . . . . . . . . . . . . . . 3.1.1 Metode Courant-Isaacson-Rees (FTBS, upwind) . 3.1.2 Metode Richardson . . . . . . . . . . . . . . . . . . 3.1.3 Metode Lax . . . . . . . . . . . . . . . . . . . . . . 3.1.4 Metode Lax - Wendroff One Step . . . . . . . . . . 3.2 Persamaan Gelombang . . . . . . . . . . . . . . . . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

6 . 6 . 7 . 9 . 9 . 12 . 12

. . . . . .

14 14 15 17 17 19 22

. . . . . .

4 Persamaan Difusi 26 4.1 Prinsip Maksimum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4.2 Fungsi Green . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 4.3 Metoda separasi variable, recall . . . . . . . . . . . . . . . . . . . . . . . . . 30 5 Metode beda hingga bagi persamaan difusi 5.1 Metode FTCS (Forward Time Center Space) . 5.2 Metode Implisit BTCS (Backward Time Center 5.3 Metode Crank - Nicholson . . . . . . . . . . . . 5.4 Metode - θ . . . . . . . . . . . . . . . . . . . . 5.5 Soal-soal . . . . . . . . . . . . . . . . . . . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

31 31 32 34 35 36

6 Persamaan Laplace & Poisson 6.1 Prinsip maksimum dan ketunggalan solusi persamaan Laplace 6.2 Persamaan Laplace pada domain persegi panjang . . . . . . . 6.3 Rumus Poisson . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Formula Green bagi solusi persamaan Laplace . . . . . . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

38 39 41 42 45

. . . . Space) . . . . . . . . . . . .

. . . . .

. . . . .

. . . . .

7 Metode beda hingga bagi persamaan Laplace 2D 46 7.1 Metode langsung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 7.2 Metode Iterasi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

1

Bab 1

SRP '12

Pendahuluan pdp Persamaan diferensial parsial adalah persamaan diferensial bagi fungsi peubah banyak u(x, y, · · · ). Orde dari pdp adalah turunan tertinggi yang muncul pada pdp tersebut. Untuk fungsi dua peubah u(x, y) bentuk umum pdp orde satu adalah: F (x, y, u(x, y), ux (x, y), uy (x, y)) = F (x, y, u, ux , uy ) = 0, sedangkan bentuk umum pdp orde dua adalah: F (x, y, u, ux , uy , uxx , uxy , uyy ) = 0.

(1.0.1)

(1.0.2)

Solusi dari pdp adalah fungsi u(x, y) yang memenuhi persamaan diferensial tersebut, untuk suatu daerah di bidang−xy.

Pada pencarian solusi persamaan diferensial biasa, kadang kala variabel bebas dan variable 3 tak bebas boleh ditukar, misalnya mencari solusi du dx = u . Untuk persamaan diferensial parsial hal ini tidak diperbolehkan. Peran variabel bebas dan variabel tak bebas tak dapat ditukar. Berbagai persamaan diferensial parsial yang penting 1. ux + ut = 0 (persamaan transport)

2. ux +uut = 0 (persamaan Burgers, merupakan bentuk khusus pers. Buckley-Leverett) 3. ut = kuxx (persamaan panas, persamaan difusi) 4. utt − c2 uxx = 0 (persamaan gelombang) 5. uxx + uyy = 0 ≡ △u = 0 (persamaan Laplace)

Semua pdp pada contoh di atas adalah pdp linier, kecuali pdp no.3.

Klasifikasi pdp Perhatikan bentuk umum persamaan diferensial parsial orde-2 berikut Auxx + 2Buxy + Cuyy = F (x, y, u, ux , uy ). Persamaan diferensial dikatakan bertipe 1. eliptik jika AC − B 2 > 0, 2

(1.0.3)

BAB 1. PENDAHULUAN PDP

3

2. parabolik jika AC − B 2 = 0, 3. hiperbolik jika AC − B 2 = 0. Tunjukkan bahwa persamaan difusi ut = kuxx bertipe parabolik, persamaan gelombang utt = c2 uxx bertipe hiperbolik. persamaan Laplace uxx + uyy = 0 bertipe eliptik,

Definisi 1.0.1. Suatu operator L dikatakan linier jika L(u + v) = Lu + Lv,

and L(cu) = cLu,

untuk setiap fungsi u, v dan untuk setiap bilangan real c. Suatu pdp berbentuk Lu = 0 dikatakan linier jika L operator linier.

SRP '12

Persamaan diferensial yang bertipe sama mempunyai perilaku solusi yang serupa. Oleh karena itu pada kuliah ini kita akan mempelajari perilaku solusi persamaan bertipe parabolik melalui bentuk kanoniknya, yaitu persamaan difusi. Begitu pula dengan perilaku solusi persamaan tipe hiperbolik melalui persamaan gelombang, dan perilaku solusi persamaan tipe eliptik melalui persamaan Laplace.

(1.0.4)

Contoh 1.0.2. Akan dibuktikan bahwa persamaan difusi adalah pdp linier. Persamaan difusi dapat dituliskan dalam bentuk Lu = 0, dengan operator L = ∂t − k∂xx . Selanjutnya L(u + v) = (∂t − k∂xx )(u + v) = ∂t u − k∂xx u + ∂t v − k∂xx v = Lu + Lv.

Coba buktikan bahwa persamaan transport dan persamaan gelombang adalah juga pdp linier.

Semua pdp pada contoh di atas dikatakan homogen, karena dapat dituliskan dalam bentuk Lu = 0. Sedangkan ux +ut = x adalah pdp tak homogen. Bentuk umum pdp tak homogen: Lu = g,

dengan g ̸= 0.

Contoh 1.0.3. 1. Buktikan prinsip superposisi solusi pdp linier. Untuk suatu pdp linier Lu = 0, dengan L operator linier, maka berlaku prinsip superposisi solusi. Jika u1 dan u2 masing-masing solusi, maka αu1 + βu2 , dengan α, β ∈ ℜ juga merupakan solusi.

2. Jika u1 solusi pdp linier tak homogen Lu = g untuk suatu g, sedangkan u0 adalah solusi pdp homogennya Lu = 0, maka u1 + u0 solusi pdp yang mana? Buktikan jawab Anda. 3. Perhatikan bahwa pdp berikut dapat dicari solusinya dengan cara persamaan diferensial biasa. (a) ux = x

BAB 1. PENDAHULUAN PDP

4

(b) uxx = 0 (c) uxx + u = 0 (d) uxy = 0 4. (a) Cari semua polinom derajat satu berbentuk p(x, y) = ax+by+c yang memenuhi uxx + uyy = 0. (b) Cari semua polinom homogen derajat dua berbentuk p(x, y) = ax2 + bxy + cy 2 yang memenuhi uxx + uyy = 0. Sketsakan solusi polinom derajat dua yang relatif sederhana.

SRP '12

Dari contoh-contoh di atas dapat disimpulkan bahwa solusi dari pdp memuat bukan hanya kon stanta sebarang melainkan fungsi sebarang. 5. Buktikan bahwa u(x, t) = f (bx−at), untuk sebarang fungsi f , memenuhi persamaan transport aux + but = 0. Misalkan diketahui a = 2, b = 1, dan { 1, 0 < x < 3 u(x, 0) = f (x) = 0, untuk x lainnya Sketsakan u(x, 0), u(x, 1), u(x, 2) pada satu sumbu koordinat. Amati apa yang terjadi. Cari solusi dari △u = 0 pada persegi satuan [0, 1]×[0, 1] yang memenuhi syarat batas u(x, 1) = sin πx dan u = 0 pada ketiga sisi lainnya. (Petunjuk: Mengingat syarat batasnya, cobalah solusi berbentuk u(x, y) = a(y) sin πx.) Catatan: Jika disubstitusikan akan diperoleh solusi u(x, y) = 6. (sinh πy/ sinh π) sin πx. Gambar di samping menunjukkan permukaan solusi persamaan Laplace u(x, y). Plot solusi persamaan Laplace dengan syarat batas tipe Dirichlet dapat dibayangkan berupa permukaan membran yang tepinya dibengkokkan mengikuti syarat batas Dirichlet yang diberikan. Diketahui persamaan Laplace beserta syarat batasnya uxx + uyy = 0, ux (0, y) = 0,

pada pers. panj. [0, 2] × [0, 1]

ux (2, y) = −2,

uy (x, 0) = 0,

uy (x, 1) = 1,

7. Tentukan solusi u(x, y) berbentuk polinom berderajat-2 yang memenuhi persamaan Laplace beserta syarat batasnya. Catatan: Plot solusi masalah ini diberikan pada gambar di samping. Yakinilah bahwa permukaan pada gambar di samping memenuhi keempat syarat batas yang diberikan. Soal Latihan 1.0.4. linier?

1. Manakah diantara operator berikut yang merupakan operator

(a) Lu = ux + xuy (b) Lu = ux + uuy (c) Lu = ux + u2y (d) Lu = ux + uy + 1

BAB 1. PENDAHULUAN PDP √ (e) Lu = 1 + x2 (cos y)ux + uyxy − [arctan(x/y)]u

5

2. Buktikan bahwa u(x, y) = f (x)g(y) adalah solusi dari pdp uuxy = ux uy untuk sebarang fungsi f dan g yang diferensiabel. 3. Buktikan melalui substitusi langsung bahwa un (x, y) = sin nx sinh ny adalah solusi dari uxx + uyy = 0 untuk setiap n > 0.

ux (0, y) = 0,

ux (1, y) = −1

uy (x, 0) = 0,

SRP '12

4. Perhatikan persamaan Laplace uxx + uyy = 0 dengan syarat batas uy (x, 1) = 1.

Carilah solusi eksak masalah di atas yang berbentuk polinom derajat dua

5. (Pembuktian solusi pdp dengan substitusi langsung) Soal-soal Kreyszig 11.1 no 2 13, 14(b). 6. (Pdp yang dapat diselesaikan sebagai ode) Soal-soal Kreyszig 11.1 no 15 - 22.

Bab 2

2.1

Persamaan Transport

SRP '12

Persamaan Type Hiperbolik

Persamaan transport, dikenal juga sebagai persamaan konveksi, seringkali muncul pada masalah transport dari berbagai substansi, misalnya polutan, gas, atau fluida lainnya. Bentuk umum persamaan transport adalah sbb: aηt + bηx = 0.

(2.1.1)

Berikut ini akan diturunkan model persamaan transport melalui konservasi massa. Perhatikan fluida yang mengalir dengan kecepatan u pada sebuah kanal dengan penampang konstan. Misalkan η(x, t) menyatakan ketinggian fluida pada posisi x saat t. Perhatikan massa fluida pada domain dengan batas kiri x dan batas kanan x + ∆x. Dalam selang waktu ∆t [akumulasi] = [lajumasuk] − [lajukeluar] (η|t+∆t − η|t )∆xb = ((uη)|x − (uη)|x+∆x ) b∆t

h ( x, t )

u

u

b

Jika dibagi dengan b∆x∆t, dan diambil limit ∆x → 0, ∆t → 0, akan diperoleh

x + Dx

ηt + (uη)x = 0.

(2.1.2)

Persamaan (2.1.2) dikenal sebagai persamaan konservasi massa. Latihan: Coba turunkan persamaan konservasi massa jika dasar kanal tak rata dan topografinya h(x) (lebar kanal tetap b).

Jika fluks u konstan, maka persamaan (2.1.2) berubah menjadi ηt + uηx = 0, suatu persamaan transport dalam variabel η(x, t) dengan koefisien konstan u. Metoda integral Sebagai alternatif, persamaan transport dapat juga diturunkan melalui formulasi integral. Massa total fluida pada domain pengamatan [x0 , x1 ] adalah ∫ x1 M= η dx. x0

Misalkan u(x, t) menyatakan kecepatan fluida di posisi x saat t. Laju aliran fluida keluar, atau fluks di posisi xi adalah η(xi , t)u(xi , t), dengan i = 0, 1. Dengan demikian laju perubahan massa fluida pada domain pengamatan [x0 , x1 ] adalah laju perubahan =fluks masuk - fluks keluar 6

BAB 2. PERSAMAAN TYPE HIPERBOLIK ∫ d x1 η(x, t) dx = η(x1 , t)u(x1 , t) − η(x0 , t)u(x0 , t). dt x0

7 (2.1.3)

Persamaan (2.1.3) merupakan persamaan konservasi massa dalam bentuk integral. Untuk memperoleh bentuk lain dari konservasi massa, maka persamaan di atas diintegralkan terhadap waktu dari t0 ke t1 dan menghasilkan ∫ x1 ∫ x1 ∫ t1 ∫ t1 η(x, t1 )dx − η(x, t0 )dx = η(x0 , t)u(x0 , t)dt − η(x1 , t)u(x1 , t)dt. (2.1.4) x0

x0

t0

t0

dan

∫ η(x1 , t)u(x1 , t) − η(x0 , t)u(x0 , t) =

x1

x0

∂ (η(x, t)u(x, t)) dx. ∂x

Selanjutnya (2.1.4) dapat dituliskan sebagai } ∫ t1 ∫ x1 { ∂ ∂ η(x, t) + (η(x, t)u(x, t)) dxdt = 0. ∂t ∂x t0 x0

SRP '12

Assumsikan η(x, t) dan u(x, t) merupakan fungsi-fungsi yang terdiferensialkan, maka ∫ t1 ∂ η(x, t1 ) − η(x, t0 ) = η(x1 , t)dt ∂t t0

Mengingat persamaan di atas berlaku untuk sebarang selang [x0 , x1 ] dan sebarang selang waktu [t0 , t1 ], maka integrannya haruslah nol, atau ηt + (uη)x = 0.

(2.1.5)

Persamaan (2.1.5) merupakan konservasi massa dalam bentuk diferensial. Dalam hal laju aliran u(x, t) sebanding dengan η(x, t), maka persamaan konservasi massa menjadi ηt + ηηx = 0

(2.1.6)

Persaman Burger terkenal akan solusinya berupa gelombang kejut (shock wave). Persamaan Burger merupakan bentuk khusus dari persamaan Buckley-Leverett, suatu model bagi aliran fluida dua phasa, misal air dan minyak. Namun kedua persamaan di atas tidak akan dibahas dalam kuliah ini. Di kuliah ini kita akan membahas kasus bila u(x, t) = c konstan, dan persamaan konservasi massa menjadi ηt + c ηx = 0.

2.1.1

Metoda karakteristik

Untuk pembahasan selanjutnya kita menggunakan persamaan transport (koefisien konstan) dengan bentuk umum berikut (perhatikan adanya perubahan notasi η → u) aux + but = 0,

(2.1.7)

dengan konstanta a, b bilangan real sebarang. Perhatikan bahwa persamaan transport di atas dapat dituliskan sebagai turunan berarah Dv u = v · ∇u = 0, dengan vektor arah ( ) a v= . b

BAB 2. PERSAMAAN TYPE HIPERBOLIK

t

Ini berarti u(x, t) solusi (3.1.1), jika dan hanya jika turunan berarah dari u(x, t) pada arah v = (a b)T sama dengan nol! Perhatikan garis-garis yang sejajar dengan vektor v pada gambar berikut. Persamaan garis-garis sejajar itu adalah t = ab x − C ≡ bx − at = C. Garis-garis ini disebut sebagai garis karakteristik.

8

bx-at=C

x

Perhatikan tiga pernyataan equivalen berikut.

2. Dv u = 0 dengan v = (a, b)T

SRP '12

1. u(x, t) solusi aux + but = 0.

3. Nilai fungsi u(x, t) konstan selama x, t terletak pada satu garis karakteristik.

Point 3 akan dielaborasi. Nilai fungsi u(x, t) konstan selama x, t terletak pada satu garis karakteristik bx − at = C. Jadi u(x, t) hanya bergantung pada nilai C pada persamaan garis karakteristik tersebut, atau u(x, t) hanya bergantung pada bx − at, sehingga solusi pers transport aux + but = 0 adalah u(x, t) = f (bx − at), dengan f suatu fungsi sebarang. Fungsi sebarang f dapat diperoleh jika diberikan syarat tambahan. Syarat tersebut bisa berupa syarat awal atau syarat batas.

Soal Latihan 2.1.1. 1. Buktikan bahwa u(x, t) = f (bx − at), untuk sebarang fungsi f , memenuhi persamaan transport aux + but = 0. Misalkan diketahui a = 2, b = 1, dan { 1, 0 < x < 3 u(x, 0) = f (x) = 0, untuk x lainnya Sketsakan u(x, 0), u(x, 1), u(x, 2) pada satu sumbu koordinat. Amati apa yang terjadi. 2. Diberikan suatu persamaan transport 3ux − 4ut = 0 dengan syarat awal { 4 − x2 , |x| < 2 u(x, 0) = 0, |x| > 2 (a) Tentukan persamaan garis karakteristiknya lalu gambarkan.

(b) Sketsakan u(x, 0), u(x, 1), u(x, 2), u(x, 3) pada satu sumbu koordinat. 3. Tentukan solusi persamaan transport 2ux − ut = 0, dengan syarat awal { sin πx 0 ≤ x ≤ 1 u(x, 0) = 0 untuk x lainnya.

Sketsakan u(x, 0), u(x, 1), u(x, 2) pada satu sumbu koordinat. Tentukan rumusan solusi u(x, t). Gunakan Maple untuk memplot u(x, 0), u(x, 1), u(x, 2) pada satu sumbu koordinat (Pastikan hasilnya sesuai). 4. Tentukan solusi persamaan transport ux + ut = 0, dengan syarat batas { 1 t≥0 u(0, t) = 0 untuk t lainnya. Sketsakan u(x, 0), u(x, 1), u(x, 2) pada satu sumbu koordinat.


9

5. Perhatikan pdp 2ut + ux = x. Persamaan di atas dapat ditulis sebagai Lu = x, dengan operator diferensial L = 2∂t + ∂x . (a) Buktikan bahwa operator L adalah operator linier, sehingga pdp yang bersesuaian juga merupakan pdp linier. (b) Tentukan satu solusi pdp tak homogen di atas. Selanjutnya mengingat bagi pdp linier berlaku sifat: solusi umum pdp tak homogen = solusi umum pdp homogen + satu solusi pdp tak homogennya, tentukan solusi umum pdp di atas.

v = ax + bt,

z = bx − at.

SRP '12

6. Perhatikan persamaan transport aux + but = 0. Terapkan transformasi koordinat, yang mengubah variabel bebas (x, t) menjadi variabel bebas baru (v, z) berikut

Gunakan aturan rantai untuk menyatakan ux dalam turunan-turunan u terhadap variabel bebas baru: uv dan uz . Lakukan hal yang sama terhadap ut . Substitusikan kedua hasilnya ke dalam persamaan transport, dan tunjukkan bahwa persamaan transport dalam variabel baru menjadi uv = 0. Selanjutnya tunjukkan bahwa solusinya adalah u(v, z) = f (z) atau u(x, t) = f (bx − at), untuk sebarang fungsi f . Metoda pencarian solusi dengan cara di atas, dikenal dengan nama metoda koordinat, dan bentuk pdp yang mana setelah transformasi koordinat dapat diselesaikan dengan pengintegralan langsung tsb dikenal sebagai bentuk normal. Untuk soal-soal solusi pdp melalui bentuk normal, lihat Kreyszig 8ed Subbab 11.4 no 11 - 17.

2.1.2

Koefisien tak konstan

Soal Latihan 2.1.2. 1. Tentukan solusi ux + tut = 0, u(0, t) = t3 . Periksa jawabnya dengan substitusi langsung. 2. Tentukan solusi ux + 2xt2 ut = 0.

3. Tentukan solusi (1 + x2 )ux + ut = 0. Sketsakan beberapa kurva karakteristiknya. √ 4. Tentukan solusi 1 − x2 ux + ut = 0, u(0, t) = t. 5. Tentukan solusi aux + buy + cu = 0. 6. Tentukan solusi ux + uy + u = ex+2y , u(x, 0) = 0.

7. Gunakan metoda koordinat untuk menentukan solusi persamaan differensial ux + 2uy + (2x − y)u = 2x2 + 3xy − 2y 2 .

2.2

Persamaan Gelombang

Getaran senar gitar mengikuti persamaan gelombang (bentuk yang paling sederhana) berikut utt = c2 uxx ,

(2.2.1)

dengan c2 = Tρ , dimana T gaya tegang senar, dan ρ rapat massa tali. Di sini u menyatakan simpangan senar dari kondisi setimbang, lihat Gambar 2.2.1. Penurunan persamaan gelombang tali dapat dilihat antara lain pada Strauss 1.3 hal. 11. Selain itu,


10

u x Gambar 2.2.1: Fungsi u(x, t) menyatakan simpangan di posisi x saat t diukur dari kondisi setimbang.

SRP '12

√ persamaan (3.2.1) juga merupakan model perambatan gelombang air, dengan c = gH, dan H menyatakan kedalaman air. Di sini c menyatakan cepat rambat gelombang air. Solusi persamaan (3.2.1) dapat diperoleh dengan beberapa cara sesuai dengan domain keberlakuan persamaan gelombang tersebut. Solusi persamaan (3.2.1) pada domain R dapat diperoleh dengan metode d’Alembert. Misalkan u merupakan suatu fungsi yang bergantung pada variabel baru ξ dan η, dengan ξ dan η dinyatakan sebagai berikut, ξ = x + ct dan η = x − ct, sehingga dapat diperoleh ∂x = ∂ξ + ∂η ; ∂t = c(∂ξ − ∂η ), utt = c2 (uξξ − 2uξη + uηη ), uxx = uξξ + 2uξη + uηη . Dengan demikian, persamaan (3.2.1) menjadi −4uξη = 0.

(2.2.2)

Solusi dari (2.2.2) adalah u(ξ, η) = F (ξ) + G(η), u(x, t) = F (x + ct) + G(x − ct),

(2.2.3)

dengan F dan G adalah sebarang fungsi differensial. Misal ditambahkan syarat awal sebagai berikut, { u(x, 0) = ϕ(x) ut (x, 0) = ψ(x). Substitusikan (2.2.3) ke (2.2.4) diperoleh : { F (x) + G(x) = ϕ(x) c(F ′ (x) − G′ (x)) = ψ(x)

(2.2.4)

{ ⇒

F ′ + G ′ = ϕ′ F ′ − G′ = 1c ψ,

sehingga akan didapatkan, 1 F = 2 ′

( ) ( ) ψ 1 ψ ′ ′ ′ ϕ + ; G = ϕ − . c 2 c

Jadi, 1 1 F (s) = ϕ(s) + 2 2c

∫

s

ψ(r)dr + k1 , 0

BAB 2. PERSAMAAN TYPE HIPERBOLIK ∫ 1 1 s G(s) = ϕ(s) − ψ(r)dr + k2 , 2 2c 0

11

sehingga, 1 1 u(x, t) = {ϕ(x + ct) + ϕ(x − ct)} + 2 2c

∫

x+ct

ψ(r)dr.

(2.2.5)

x−ct

Persamaan (2.2.5) dikenal sebagai rumus d’Alembert.

c. u(x, 0) = φ(x), dan ut (x, 0) = ψ(x). Amati maxx u(x, t).

SRP '12

Soal: Plot kurva-kurva solusi persamaan gelombang utt − c2 uxx = 0 dengan tiga tipe pilihan syarat awal berikut. { 1 − x2 , untuk |x| ≤ 1 a. u(x, 0) = φ(x) = dan ut (x, 0) = 0. 0, untuk x lainnya { 1, untuk |x| ≤ 1 b. u(x, 0) = 0, dan ut (x, 0) = ψ(x) = . 0, untuk x lainnya

Gambar 2.2.2: Garis-garis karakteristik persamaan gelombang.

t

t (x,t) x+ct=xo

x-ct=xo

x x+ct

x-ct (a)

x

xo (b)

Gambar 2.2.3: (a). Daerah kebergantungan (domain of dependence) dari titik (x, t), (b). Daerah pengaruh (domain of influence) dari titik (x0 , 0).

Perhatikan bahwa rumus solusi d’Alembert berlaku untuk persamaan gelombang pada domain tak hingga. Apabila domain keberlakuan persamaan gelombang hanya berupa setengah selang, misalnya (0, ∞) atau (−∞, 0) maka rumus d’Alembert dapat dimodifikasi. Diskusi lengkap mengenai hal ini ada pada Bab 3.2 Strauss. Apabila domain keberlakuan persamaan gelombang adalah berupa domain berhingga, maka cara penyelesaiannya adalah menggunakan metoda separasi variabel, Bab 4 Strauss.


12

Jika dipunyai syarat batas u(0, t) = 0 (tipe Dirichlet) maka syarat awal u(x, 0) = ϕ(x) pada (0, ∞) diperluas menjadi fungsi ganjil. Selanjutnya karena pada persamaan gelombang, gelombang menjalar mengikuti garis karakteristik, kita dapatkan sketsa gelombang untuk t = t0 , lihat Gambar 1.4. Hasil tersebut sesuai dengan kondisi gelombang yang merambat pada tali dengan ujung kiri terikat.

2.2.1

Energi

SRP '12

Jika dipunyai syarat batas ux (0, t) = 0 (tipe Neumann) maka syarat awal u(x, 0) = ϕ(x) pada (0, ∞) diperluas menjadi fungsi genap. Selanjutnya karena pada persamaan gelombang, gelombang menjalar mengikuti garis karakteristik, kita dapatkan sketsa gelombang untuk t = t0 , lihat Gambar 1.5. Hasil tersebut sesuai dengan kondisi gelombang yang merambat pada tali dengan ujung kiri bebas.

Berikut ini akan ditunjukkan bahwa solusi (3.2.1) memenuhi kekekalan energi. Perhatikan bahwa energi kinetik sistem ∫ 1 KE = ρ u2t dx, 2 sedangkan energi potensialnya = T × pertambahan panjang, sehingga ∫ 1 P E = T × (ds − dx) = T u2x dx. 2 sehingga energi totalnya adalah 1 E = KE + P E = 2 Kemudian, dE dt

=

∫ (ρu2t + T u2x )dx.

∫

(ρut utt + T ux uxt )dx ∫ 2 = ∫ (c ρut uxx + T ux uxt )dx = (−T utx ux + T ux uxt )dx

dE = 0. dt yang berarti sepanjang evolusi energi E bernilai konstan.

(Strauss, Section 2.2 no. 5) Periksa perubahan energi E(u) untuk u(x, t) solusi persamaan gelombang yang memperhitungkan faktor redaman utt = c2 uxx − rut , r > 0.

2.3

Well-posed problem

Suatu persamaan diferensial parsial dengan syarat awal dan syarat batas (atau syaratsyarat lainnya) dikatakan well-posed jika mempunyai solusi tunggal (existence and uniqueness) dan solusi tersebut stabil terhadap perubahan nilai awal (stability). Contoh 2.3.1. Perhatikan persamaan gelombang utt − c2 uxx = 0, untuk x ∈ R, t > 0 dengan syarat awal u(x, 0) = Φ(x), ut (x, 0) = Ψ(x) bersifat well-posed. Penjelasan, mengingat adanya rumus solusi d’Alembert dari persamaan gelombang, maka jelas memenuhi existence &


13

uniqueness. Berikut adalah sekilas mengenai kestabilan solusi terhadap syarat awal u(x, 0) = Φ(x), ut (x, 0) = 0. Misalkan ui (x, t), i = 1, 2 berturut-turut adalah solusi persamaan gelombang dengan syarat awal u(x, 0) = Φi (x), i = 1, 2. 1 1 (u2 − u1 )(x, t) = (Φ2 − Φ1 )(x − ct) + (Φ2 − Φ1 )(x + ct) 2 2

||Φ2 − Φ1 || < δ =⇒ ||u2 − u2 || < ε. Perhatikan sistim persamaan linier Am×n Un×1 = bm×1 .

SRP '12

Kestabilan dipenuhi jika untuk setiap ε > 0, terdapat δ > 0 sehingga

• Jika m > n yang berarti jumlah baris lebih besar dari jumlah kolom, maka SPL di atas bersifat overdetermined, jadi SPL tak mempunyai solusi.

• Jika m < n yang berarti jumlah baris lebih kecil dari jumlah kolom, maka SPL di atas bersifat underdetermined, jadi SPL tak mempunyai tak hingga banyak solusi.

• Jika m = n dan matriks A bersifat singular (det A = 0), maka SPL tak punya solusi.

• Jika m = n dan nilai eigen dari matriks A kecil sekali, maka SPL tidak memenuhi konsep kestabilan. Untuk setiap kasus yang diuraikan di atas, SPL bersifat ill-posed.

Sebagai catatan: masalah difusi untuk t ≥ 0 bersifat well-posed, namun backward-diffusion tidak.

Bab 3

SRP '12

Metoda Beda Hingga pada Persamaan Tipe Hiperbolik

Sebuah persamaan diferensial apabila didiskritisasi dengan metode beda hingga akan menjadi sebuah persamaan beda. Jika persamaan diferensial parsial mempunyai solusi eksak u(x, t), maka persamaan beda akan mempunyai solusi hampiran u(xj , tn ). Kaitan antara pdp dengan persamaan beda, dan solusi u(x, t) dan solusi hampiran u(xj , tn ) tak lain adalah konsep kekonsistenan, kestabilan, dan kekonvergenan suatu persamaan beda, lihat Gambar 3.0.1. kekonsistenan

kekonvergenan

Gambar 3.0.1: Skema hubungan antara pdp dan persamaan beda serta solusi-solusinya.

3.1

Persamaan Transport

Sebagai langkah awal pengenalan metoda beda hingga, akan digunakan pdp yang paling sederhana yaitu persamaan transport ut + dux = 0, untuk (x, t) ∈ [0, L] × [0, T ]

(3.1.1)

beserta syarat awalnya u(x, 0) = f (x), untuk 0 ≤ x ≤ L

(3.1.2)

yang telah kita ketahui mempunyai solusi eksak u(x, t) = f (x − dt). Persamaan beda hingga dikatakan stabil jika persamaan beda menghasilkan solusi unj yang berhingga. Persamaan beda dikatakan konsisten terhadap pdpnya jika selisih antara persamaan beda dengan pdpnya (suku-suku truncation error ) menuju nol jika lebar grid menuju nol, ∆x → 0, ∆t → 0. Persamaan beda hingga dikatakan konvergen jika solusi 14

BAB 3. METODA BEDA HINGGA PADA PERSAMAAN TIPE HIPERBOLIK

15

persamaan beda mendekati solusi pdp jika lebar grid menuju nol, ∆x → 0, ∆t → 0. Teorema Equivalensi Lax: Untuk suatu masalah nilai awal yang well-(properly-)posed, jika suatu persamaan beda konsisten dan stabil, maka persamaan beda tersebut pastilah konvergen.

3.1.1

Metode Courant-Isaacson-Rees (FTBS, upwind)

un+1 − unj j ∆t

+d

unj − unj−1 = 0. ∆x

atau = (1 − C)unj + Cunj−1 , un+1 j dengan C ≡

d△t △x

SRP '12

Perhatikan selang [0, L] yang dipartisi dengan ukuran △x dengan titik-titik partisi xj = (j − 1)△x, untuk j = 1, 2, · · · , N x. Selang [0, T ] dipartisi dengan ukuran △t dengan titik-titik partisi tn = (n − 1)△t, untuk n = 1, 2, · · · , N t. Akan dicari u(xj , tn ) untuk j = 1, 2, · · · , N x, dan n = 1, 2, · · · , N t. Akan digunakan notasi unj ≡ u(xj , tn ). Metode Courant-Isaacson-Rees tak lain adalah Metode FTBS (Forward Time Backward Space) dikenal juga sebagai metode upwind. sehingga persamaan beda untuk persamaan transport (3.1.1) adalah: (3.1.3)

(3.1.4)

suatu bilangan Courant. Metode ini memiliki akurasi O(△t, △x).

Gambar 3.1.1: Stencil metoda upwind untuk persamaan transport.

Perhatikan domain perhitungan [0, L] × [0, T ] yang telah dipartisi. Syarat awal (3.1.2) telah diketahui sehingga u1j , untuk j = 1, 2, · · · , N x telah diketahui. Untuk satu time-step berikutnya dengan 3.1.4 kita dapat menghitung u2j , untuk j = 2, · · · , N x, namun u21 tak dapat dihitung. Ini berarti untuk pendekatan numerik ini kita membutuhkan syarat batas kiri. Jika digunakan syarat batas kiri u(0, t) = 0 atau un1 = 0 untuk n = 1, 2, · · · , N t, maka FTBS dapat digunakan untuk menentukan unj , untuk setiap j = 1, 2, · · · , N x, n = 1, 2, · · · , N t. Perhatikan bahwa FTBS pada persamaan transport tidak membutuhkan syarat batas kanan. Catatan: di sini syarat batas kiri kita tambahkan melulu untuk keperluan perhitungan numerik. Pilihan untuk syarat batas kiri tersebut tentu berragam. Jika kita ingin mensimulasikan situasi tertentu, pilihan syarat batas ini menjadi hal yang penting untuk dicermati. Lihat soal-soal pada akhir bab ini.


16

Kestabilan: Metode berikut ini disebut metode kestabilan von Neumann1 Substitusikan unj = ρn eiaj ke dalam (3.1.4). Setelah dibagi dengan unj , akan diperoleh ρ = 1 − C(1 − e−ia ). atau ρ = 1 − C(1 − cos a) − iC sin a. Persamaan beda stabil jika dan hanya jika |ρ| < 1 atau (

)2 1 − 2C sin2 a2 + (C sin a)2 2 sin4 a + 4C 2 sin2 a cos2 a 1 − 4C(sin2 a2 + 4C 2 2 2) ) ( 4C sin2 a2 −1 + C sin2 a2 + C cos2 a2 4C(C − 1) sin2 a2

≤ 1, ≤ 1, ≤ 0, ≤ 0.

SRP '12

|ρ|2 =

d∆t ≤ 1. Ketaksamaan terakhir dipenuhi untuk setiap a ∈ R jika dan hanya jika 0 ≤ C = ∆x d∆t Jadi, syarat kestabilan metode FTBS adalah 0 ≤ ≤ 1. ∆x Konsistensi: Perhatikan ekspansi Taylor dari un+1 dan unj−1 masing-masing di sekitar j n uj berikut. un+1 j unj−1

=

unj

n ∂u n 1 2 ∂ 2 u + ∆t + ∆t + .... ∂t j 2 ∂t2 j

=

unj

2 n ∂u n 1 2 ∂ u + ∆x − ∆x + .... ∂x j 2 ∂x2 j

(3.1.5)

(3.1.6)

Substitusikan (3.1.5) dan (3.1.6) ke dalam persamaan beda (3.1.4), dan setelah menggunakan hubungan ( ) ( ) 2 ∂2u ∂ ∂u ∂ ∂u 2∂ u −d = −d = d = . ∂t2 ∂t ∂x ∂x ∂t ∂x2 akan diperoleh modified differential equation sebagai berikut: ( ) n ∂u n ∂u n 1 ∂ 2 u +d + d(d∆t − ∆x) + .... = 0. ∂t j ∂x j 2 ∂x2 j

(3.1.7)

Suku pertama pada (3.1.7) tak lain adalah persamaan transport yang akan diselesaikan. Suku kedua dan seterusnya pada (3.1.7) tak lain adalah suku tambahan yang kita dapatkan saat kita bekerja dengan persamaan beda (3.1.4), dan disebut truncation term. n 1 ∂ 2 u truncation term = d(d∆t − ∆x) 2 ∂x2 j Perhatikan bahwa jika ∆t → 0 dan ∆x → 0, maka truncation term→ 0. Jadi metode FTBS konsisten terhadap persamaan transport. Catatan: Modified differential equation adalah persamaan yang diperoleh dari persamaan beda, yang mana beberapa sukunya diekspansikan ke dalam deret Taylor, sehingga akan ditemui suku-suku pada pdp semula. Suku-suku pada modified differential equation yang 1

Terdapat tiga metode untuk menguji kestabilan suatu persamaan beda, yaitu von Neumann stability analisis, metode matriks, dan discrete perturbation method.


17

tidak terdapat pada pdp semula disebut truncation term. Pengamatan lebih lanjut, saat kita bekerja dengan persamaan beda (3.1.4) sebenarnya kita bukan menyelesaikan persamaan transport (3.1.1) melainkan menyelesaikan modified differential equation (3.1.7). Suku pertama dari truncation term akan nol jika dan hanya jika d∆t = ∆x (C = 1). Suku kedua dari truncation term adalah ( ) 1 3 3 1 2 1 3 − d ∆t + d ∆x∆t − d∆x ) uxxx |nj 3 2 6

3.1.2

Metode Richardson

SRP '12

dan akan bernilai nol jika C = 1. Ternyata hal ini juga berlaku untuk suku ketiga dan seterusnya pada truncation term. Jadi jika C = 1 persamaan beda (3.1.4) ekuivalen dengan (3.1.1). Dengan kata lain metode FTBS pada persamaan transport dengan C = 1 akan menghasilkan solusi eksak, dan uraian di atas adalah buktinya. Jika C < 1, maka ∂2u , dan akan menimbulkan error nusuku pertama truncation term berupa suku difusi ∂x2 merik berupa damping.

Metode Richardson biasa disebut sebagai Metode FTCS (Forward Time Center Space), akurasim metode ini O(∆t, ∆x2 ). Persamaan beda hingga untuk persamaan transport (3.1.1) adalah sebagai berikut: un+1 − unj j ∆t

+d

unj+1 − unj−1 = 0. 2∆x

(3.1.8)

Gambarkan stencil metode Richardson.

Kestabilan: Substitusikan unj = ρn eiaj ke dalam (3.1.8), kemudian dibagi dengan unj , maka persamaan (3.1.8) menjadi: eia − e−ia ρ−1 +d = 0, ∆t 2∆x sehingga diperoleh

∆t i sin a, ∆x Perhatikan bahwa |ρ| > 1 untuk setiap a ∈ R, ini berarti metode FTCS / metode Richardson selalu tidak stabil. ρ=1−d

3.1.3

Metode Lax

Metode ini merupakan perbaikan dari metode Richardson, yang mana nilai dari unj diganti menjadi rata-rata dari unj+1 dan unj−1 , sehingga n+1 ∂u n uj − = ∂t j

1 2

(

unj+1 + unj−1

)

∆t

∂u n unj+1 − unj−1 = + O(∆x2 ). ∂x j 2∆x

+ O(∆t),

(3.1.9)

(3.1.10)


18

Substitusikan (3.1.9) dan (3.1.10) ke (3.1.1) sehingga diperoleh persamaan beda sebagai berikut: ( ) 1 n n u + u un+1 − unj+1 − unj−1 j+1 j−1 j 2 +d = 0. (3.1.11) ∆t 2∆x Akurasi metode ini O(∆t, ∆x2 ). Gambarkan stencil metode Lax. Kestabilan: Substitusikan unj = ρn eia∆xj ke dalam (3.1.11), sehingga diperoleh:

dengan C =

) 1 ( ) 1 ( ia∆x e + e−ia∆x − C eia∆x − e−ia∆x , 2 2

d∆t . Kemudian juga dapat diperoleh, ∆x ρ = cos √ a∆x − iC sin a∆x, |ρ| = cos2 a∆x + C 2 sin2 a∆x,

Jadi, jika C =

SRP '12

ρ=

d∆t ≤ 1 maka |ρ| ≤ 1 untuk setiap a, jadi skema stabil. ∆x

Konsistensi: Perhatikan ekspansi Taylor dari un+1 dan unj±1 masing-masing di sekitar j n uj berikut. un+1 j unj±1

=

unj

n ∂u n 1 2 ∂ 2 u + ∆t + ∆t + .... ∂t j 2 ∂t2 j

=

unj

2 n ∂u n 1 2 ∂ u ± ∆x + ∆x + .... ∂x j 2 ∂x2 j

(3.1.12)

(3.1.13)

Substitusikan (3.1.12) dan (3.1.13) ke dalam 3.1.11 akan menghasilkan unj

n 2 n ∂u n 1 2 ∂ 2 u 1 2 ∂ u n + ∆t + ∆t + ... = uj + ∆x + ... 2 ∂t j 2 ∂t2 j 2 ∂x j ( ) 3 n 1 ∂u n 3 ∂ u + (∆x) + ... −C ∆x ∂x j 3! ∂x3 j n 2 n 3 n ∂u n 1 2 ∂ 2 u ∂u n 1 1 2 ∂ u 2 ∂ u ∆t + ∆t + ... = ∆x + ... − d∆t − d∆t∆x ∂t j 2 ∂t2 j 2 ∂x2 j ∂x j 3! ∂x3 j (3.1.14)

Perhatikan bahwa di sini berlaku pers. transport (3.1.1), maka n 2 n ∂ 2 u 2 ∂ u =d . ∂t2 j ∂x2 j

Sehingga diperoleh modified difference equation berikut ( ) ( 2 ) ∂ 2 u n ∂u ∂u 1 2 + O(△x3 , △t3 ) ∆t +d + ∆x C − 1 ∂t ∂x 2 ∂x2 j Jelas bahwa suku truncation term akan menuju nol jika ∆x → 0 dan ∆t → 0. Jadi persamaan beda (3.1.11) konsisten.


19

Lebih jauh dapat dibuktikan bahwa, suku pertama truncation term yang berorde-2: O(∆x2 ) akan bernilai nol jika dan hanya jika C = 1. Begitu juga dengan suku-suku truncation term berorde lebih tinggi (cek!). Jika C < 1, suku berorde-2 tersebut akan menimbulkan efek difusi. Selain itu, suku-suku truncation term yang merupakan turunan orde genap akan memberikan efek error numerik berupa difusi atau damping, disipasi, sedangkan sukusuku yang merupakan turunan orde ganjil akan memberikan efek error numerik berupa efek dispersi.

3.1.4

Metode Lax - Wendroff One Step

Pada metode ini digunakan hampiran sebagai berikut: 1 un+1 = unj + ∆t ut |nj + ∆t2 utt |nj + O(∆t3 ). j 2

SRP '12

Terdapat sebuah aturan umum yang mengatakan bahwa approksimasi orde-1 bagi ∂t akan menghasilkan error numerik berupa difusi sedangkan aproksimasi orde-2 bagi ∂t akan menghasilkan error numerik berupa dispersi. Pada uraian selanjutnya akan dibahas metode beda hingga berorde-2.

(3.1.15)

Kemudian dari persamaan transport dapat diperoleh bahwa ut = −dux dan utt = d2 uxx , sehingga persamaan (3.1.15) menjadi n 1 un+1 = unj − ∆t dux |nj − ∆t2 d2 uxx j + O(∆t3 ). j 2

(3.1.16)

Aproksimasi ux dan uxx dengan center difference sehingga diperoleh persamaan beda dari metode Lax-Wendroff One Step bagi persamaan transport sebagai berikut: un+1 = unj − j dengan C =

) C2 ( n ) C( n uj+1 − unj−1 + uj+1 − 2unj + unj−1 + O(∆t3 , ∆x3 ), 2 2

(3.1.17)

d∆t . Akurasi dari metode Lax - Wendroff One Step: O(∆t2 , ∆x2 ). ∆x

Gambar 3.1.2: Stencil metode Lax-Wendroff One Step untuk persamaan transport.

Soal Latihan 3.1.1. 1. Implementasikan metoda upwind untuk persamaan ut + ux = 0 dengan syarat awal u(x, 0) = sin 8πx. Gunakan ∆x = 2.5×10−2 dan ∆t = 2×10−2 . Gambarkan solusi numerik bersama-sama dengan solusi eksaknya saat t = 1 dan hasilkan Gambar 3.1.3. Perhatikan bahwa metoda upwind memuat error berupa damping.

20

SRP '12


Gambar 3.1.3: Solusi numerik beserta solusi eksak dari metode upwind saat t = 1.

d∆t ≤ 1. 2. Buktikan syarat kestabilan metoda Lax-Wendroff: C = ∆x Tentukan suku pertama truncation term dari metoda Lax-Wendroff. Kerjakan Soal di atas dengan menggunakan metoda Lax-Wendroff. Gambarkan solusi numerik bersama-sama dengan solusi eksaknya saat t = 1 dan hasilkan Gambar 3.1.4. Perhatikan bahwa error damping berkurang, namun solusi numerik memiliki error phase yang cukup nyata.

Gambar 3.1.4: Solusi numerik beserta solusi eksak dari metode Lax-Wendroff.

3. Terapkan metoda upwind dan metode Lax-Wendroff untuk mensimulasikan perambatan gelombang awal berupa gelombang persegi, dan hasilkan gambar seperti Gambar 3.1.5.

4. Perhatikan persamaan transport ut +bux = 0, dengan b konstanta real sebarang yang tidak diketahui tandanya, dan syarat awal u(x, 0) = f (x). Perhatikan persamaan beda dengan akurasi orde satu berikut 1 n+1 b+ n b− n (ui − unj ) = − (ui − uni−1 ) − (u − uni ) ∆t ∆x ∆x i+1 dengan b+ ≡ max(b, 0) dan b− ≡ min(b, 0). Implementasikan persamaan beda di atas. Pilih sendiri syarat awal f (x).

21

SRP '12


Gambar 3.1.5: Perbandingan antara solusi eksak dan numerik dari metoda upwind (kiri) dan metode Lax-Wendroff (kanan) sesudah 25, 50 dan 75 timestep.

(a) Jika b > 0, syarat batas apa yang diperlukan? pilih sendiri. Hasil apa yang Anda lihat? (b) Pertanyaan yang sama seperti di atas untuk b < 0. Berikut ini adalah metode-metode lainnya.

Soal Latihan 3.1.2. 1. Metode Leapfrog memiliki akurasi O(∆t2 , ∆x2 ). Persamaan beda untuk persamaan transport (3.1.1) dengan metode ini adalah un+1 − un−1 j j

unj+1 − unj−1 +d = 0, 2∆t 2∆x atau dapat ditulis sebagai berikut: un+1 = un−1 − C(unj+1 − unj−1 ), j j

(3.1.18)

d∆t . Gambarkan stencil metoda Leapfrog. Buktikan bahwa syarat ∆x d∆t kestabilan Metode Leapfrog adalah C = ≤ 1. Buktikan bahwa suku pertama ∆x truncation term berupa dispersi (memuat suku uxxx ). Detailnya dapat dilihat pada Hoffmann hal. 588.

dengan C =

2. (Metode Lax-Richtmyer two step) memiliki keakuratan O(∆t2 , ∆x2 ). Pada metode ini digunakan dua hampiran sebagai berikut: n+1 un+2 = unj − C(un+1 j j+1 − uj−1 ),

(3.1.19)


22

C 1 (3.1.20) un+1 = (unj+1 + unj−1 ) − (unj+1 − unj−1 ), j 2 2 d∆t dengan C = . Substitusikan (3.1.20) ke dalam (3.1.19) untuk memperoleh per∆x samaan beda berikut ini, C n C2 n (uj+2 − unj−2 ) + (u − 2unj + unj−2 ). 2 2 j+2

un+2 = unj − j

(3.1.21)

3.2

Persamaan Gelombang

Perhatikan persamaan gelombang utt − c2 uxx = 0,

T ρ

c2 =

(3.2.1)

dengan syarat awal u(x, 0) = ϕ(x),

SRP '12

Perhatikan bahwa persamaan beda (3.1.21) tepat sama dengan persamaan beda dari metoda Lax-Wendroff untuk lebar selang 2∆t dan 2∆x, yaitu (3.1.17). Dengan demikian, pembahasan mengenai kestabilan, kekonsistenan, dll. pada metode ini sama persis dengan pembahasan pada metode Lax-Wendroff one step. Gambarkan stencil metode Lax - Richtmyer two step.

ut (x, 0) = ψ(x).

(3.2.2)

Terapkan beda pusat pada kedua sukunya menghasilkan persamaan beda dengan akurasi O(∆t2 , ∆x2 ) berikut un+1 − 2unj + un−1 j j ∆t2

un 2 j+1

=c

− 2unj + unj−1 , ∆x2

(3.2.3)

atau dapat ditulis sebagai berikut: un+1 = S(unj+1 + unj−1 ) + 2(1 − S)unj − un−1 , j j dengan S =

(3.2.4)

c2 (∆t)2 . (∆x)2

S

S 2-2S

-1

Gambar 3.2.1: Stencil metode beda hingga untuk persamaan gelombang.

Perhatikan bahwa persamaan beda (3.2.4) membutuhkan dua baris syarat awal, sementara permasalahan kita hanya mempunyai syarat awal (3.2.2), yang berarti u0j , untuk j = 1, · · · , N x. Bagaimana mendapatkan u1j ? Perhatikan yang berikut, terapkan beda pusat dengan akurasi O(∆t2 ) pada ut |0j : u1j − u−1 j 2∆t

= ψj .

(3.2.5)


23

Persamaan (3.2.4) untuk n = 0 menghasilkan u1j = S(u0j+1 + u0j−1 ) + 2(1 − S)u0j − u−1 j , 1 dan karena −u−1 j = 2∆tψj − uj , maka dapat diperoleh:

u1j =

S 0 (ϕ + ϕ0j−1 ) + (1 − S)ϕ0j + ∆tψj . 2 j+1

(3.2.6)

Jadi, persamaan beda (3.2.3) dapat diterapkan dengan u0j = ϕj .

SRP '12

(3.2.7)

dan (3.2.6) sebagai dua baris nilai awal. Hampiran orde rendah bagi ut |0j , yaitu u1j − u0j = ψj + O(∆t), ∆t sehingga u1j = ϕj + ∆tψj + O(∆t2 ),

dapat digunakan sebagai nilai awal, namun akurasi O(∆t) tidak memenuhi syarat karena akurasi O(∆t) yang rendah ini akan masuk dan ’tersebar’ ke dalam domain perhitungan sehingga akurasi dari keseluruhan perhitungan hanya dapat dikatakan berorde-1 saja, dan akurasi orde-2 dari persamaan beda (3.2.4) menjadi sia-sia. Kestabilan: Substitusikan unj = ρn η j dimana η = eia∆x ke dalam (3.2.4) sehingga diperoleh 1 1 ρ − 2 + = S(η − 2 + ) = 2S(cos a∆x − 1) = 2p, ρ η yang dapat disederhanakan menjadi ρ2 − 2(1 + p)ρ + 1 = 0. Dengan demikian, diperoleh akar-akar berikut ini √ ρ1,2 = (1 + p) ± p2 + 2p , p ≤ 0. Nilai dari ρ1,2 terbagi menjadi dua kasus yaitu:

1. Jika p2 + 2p > 0 , p < −2, diperoleh bahwa ρ1,2 bernilai real dan salah satu diantaranya bernilai < −1. Jadi, skema tidak stabil. √ 2. Jika p2 + 2p < 0 , − 2 < p ≤ 0, diperoleh bahwa ρ1,2 = (1 + p) ± i −p2 − 2p merupakan bilangan kompleks dengan |ρ1,2 | = 1. Jadi ρ1,2 = cos θ + i sin θ. 3. Jika p = −2, maka akan dipeoleh ρ = −1.

Dengan demikian, skema beda hingga stabil jika p ∈ [−2, 0], ∀a. Hal ini akan terjadi jika dan hanya jika: −2 ≤ S(cos a∆x − 1) ≤ 0 , ∀a dan karena −2 ≤ cos a∆x − 1 ≤ 0 , ∀a, maka batasan untuk S adalah −2 ≤ −2S ≤ 0 ⇔ 0 < S ≤ 1.


24

Jadi syarat kestabilan untuk skema ini adalah ∆t 2 ) ≤ 1. (3.2.8) ∆x Ilustrasi garis-garis karakteristik pada persamaan gelombang dapat dilihat pada Gambar 3.2.2. Pada gambar tersebut dapat dilihat bahwa jika untuk S = 1, yang berarti ∆t = c∆x, maka garis karakteristik melewati titik-titik grid. Sedangkan jika S < 1, yang berarti ∆t < c∆x, maka garis karakteristik lebih curam yang berarti kecepatan gelombang hasil numerik lebih kecil daripada kecepatan sebenarnya. Sehingga syarat kestabilan persamaan gelombang (3.2.8) mengartikan bahwa numerical propagation speed tidak bisa lebih besar dari continuous propagation speed. (Strauss hal 204.)

(a)

(b)

SRP '12

S = c2 (

Gambar 3.2.2: Garis karakteristik persamaan gelombang untuk (a) S = 1 ⇔ ∆t = c∆x, (b) S < 1 ⇔ ∆t < c∆x.

Konsistensi: Perhatikan dua uraian Taylor berikut:

1 1 un±1 = unj ± ∆t ut |nj + ∆t2 utt |nj ± O(∆t3 ) + ∆t4 u4t |nj + .... j 2 2 1 1 unj±1 = unj ± ∆x ux |nj + ∆x2 uxx |nj ± O(∆x3 ) + ∆x4 u4x |nj + .... 2 2 Dari persamaan (3.2.9) dan (3.2.10) diperoleh: ( ) 1 1 n n n n n 2 4 uj+1 + uj−1 = 2 uj + ∆x uxx |j + ∆x u4x |j + .... 2 4! ) ( 1 4 1 2 n n n+1 n−1 n uj + uj = 2 uj + ∆t utt |j + ∆t u4t |j + .... 2 4!

(3.2.9)

(3.2.10)

(3.2.11) (3.2.12)

Kemudian substitusikan (3.2.11) dan (3.2.12) ke dalam persamaan beda (3.2.4) dan dengan sedikit manipulasi aljabar akan didapatkan 2 4 2 ∆t u4t |nj − 2Sunj − S∆x2 uxx |nj − S ∆x4 u4x |nj = 2(1 − S)unj , 4! 4! yang dapat disederhanakan menjadi 2unj + ∆t2 utt |nj +

2 ((∆t)4 u4t − S(∆x)4 u4x ) |nj . 4! Jelas bahwa skema konsisten. Lebih jauh lagi suku pertama truncation terms-nya adalah (∆t)2 (utt − c2 uxx ) |nj +

2 2 2 ((∆t)4 u4t − S(∆x)4 u4x ) = ((∆t)4 c2 u4x − S(∆x)4 u4x ) = (∆x)4 (S 2 − S)u4x , 4! 4! 4! yang berupa suku difusi, dan akan bernilai nol jika dan hanya jika S 2 − S = 0 atau c2 (∆t)2 = 1. S= (∆x)2


25

Soal Latihan 3.2.1. 1. (Standing wave) Untuk k tertentu, tentukan ω sehingga fungsi u(x, t) = cos kx cos ωt solusi bagi persamaan gelombang utt − c2 uxx = 0. Gunakan Maple untuk mensimulasikan solusi standing wave tersebut. 2. Perhatikan persamaan gelombang utt = uxx , untuk 0 ≤ x ≤ 10, dengan syarat batas ux (0, t) = 0 dan u(10, t) = 0, dan syarat awal ut (x, 0) = 0 dan { 2 2 2 16 (x − 3) (x − 7) , 3 ≤ x ≤ 7 u(x, 0) = 0, untuk x lainnya

SRP '12

(a) Terapkan skema beda hingga orde-2 untuk persamaan gelombang. Gunakan pula hampiran orde-2 untuk menaksir u(x, ∆t), juga untuk syarat batas kiri. (b) Jika dipilih ∆x = 1, diperoleh suatu hampiran bagi simpangan awal berikut u(x, 0) =

0

0

0

0

1

2

1

0

0

0

0

Hitung u(x, t) untuk beberapa selang waktu (hand-calculation), pilih ∆t = 1. Apa yang Anda lihat pada batas?

(c) Implementasikan dan simulasikan perpecahan gundukan awal sampai gelombang pecahannya menabrak batas kanan dan kiri, kemudian berbalik. Amatilah! (d) Kerjakan soal (c) namun simpangan awal nol, dan kecepatan awal { 1, |x − 5| ≤ 1 u(x, 0) = 0, x lainnya.

Bab 4

SRP '12

Persamaan Difusi Persamaan difusi adalah persamaan tipe hiperbolik. Dalam uraian berikut ini persamaan difusi akan diformulasikan berdasarkan prinsip kesetimbangan energi. Perhatikan masalah distribusi panas pada suatu batang penghantar yang panjangnya L. Misalkan u(x, t) menyatakan suhu pada posisi x saat t. Perhatikan elemen batang dengan batas kiri x dan batas kanan x + ∆x, internal energi pada elemen tersebut adalah ρc u(¯ x, t)∆x, dengan c, ρ berturut-turut menyatakan specific heat, rapat massa batang. Jika dimisalkan fluks aliran panas sebagai fungsi f (u), maka selama selang waktu ∆t, perubahan internal energi pada elemen batang adalah ρc(u|t+∆t − u|t )∆x = (f (u)|x − f (u)|x+∆x )∆t

f (u )

f (u )

Jika kedua ruas dibagi dengan ∆t∆x, dan setelah diambil limitnya, diperoleh:

ρcut + (f (u))x = 0. (4.0.1) x + Dx Persamaan 4.0.1 dengan ρc dianggap 1 atau ut +(f (u))x = 0 merupakan bentuk umum persamaan konservasi bagi masalah-masalah konservasi lainnya, misalnya konservasi massa. Tentunya setiap proses penurunan persamaan konservasi tertentu memerlukan formulasi yang cermat. Jika selanjutnya diterapkan hukum Fick’s: fluks energi sebanding dengan gradien temperatur: fluks= −κux , dengan κ menyatakan heat conductivities. maka persamaan (4.0.1) menjadi persamaan difusi ut = kuxx ,

(4.0.2)

κ . Konstanta difusi k > 0 menyatakan koefisien penghantar panas, dan dengan k = ρc 2 berdimensi [L ]/[T ]. Tanda negatif pada hukum Fick’s berarti arah aliran searah dengan arah berkurangnya suhu.

Alternatif lain, persamaan difusi dapat ∫ x+∆xjuga diperoleh melalui formulasi integral. Internal energi dari elemen batang adalah x ρcu dx, sedangkan laju perubahannya bergantung pada fluks masuk dikurangi fluks ke luar d dt

∫

x+∆x

ρcu dx = f (u)|x − f (u)|x+∆x .

x

26

BAB 4. PERSAMAAN DIFUSI

27

Perhatikan bahwa di sini batas kiri dan kanan x dan x + ∆x tidak bergerak sehingga d dt ∫

∫

∫

x+∆x

ρcu dx = x x+∆x

x+∆x

ρcut dx = x

∫ x+∆x x

(−f (u))x dx

ρcut + (f (u))x dx = 0. x

ut = kuxx ,

SRP '12

Karena relasi di atas berlaku untuk setiap x dan ∆x, maka berarti haruslah integran bernilai nol, atau u(x, t) memenuhi persamaan (4.0.1). Jika selanjutnya dimisalkan berlaku hukum Ficks: f (u) = −κux , dengan κ > 0, maka u(x, t) memenuhi persamaan difusi (4.0.3)

κ dengan k = ρc . Selain sebagai persamaan pengatur masalah distribusi panas, persamaan difusi juga merupakan persamaan pengatur proses melarutnya tinta dalam zat pelarut tertentu, dan konstanta difusi k menyatakan seberapa cepat tinta dapat terlarut.

Terdapat tiga tipe syarat batas:

1. Syarat batas tipe Dirichlet: jika nilai u diketahui. Pikirkan situasi yang bersesuaian dengan syarat batas kiri Dirichlet homogen: u(0, t) = 0 (masalah panas & masalah tinta)

2. Syarat batas tipe Neumann: jika turunan dari u diketahui. Pikirkan situasi yang bersesuaian dengan syarat batas kanan Neumann homogen: ux (L, t) = 0 (masalah panas & masalah tinta) 3. Syarat batas tipe Robin: kombinasi dari Dirichlet dan Neumann. Pikirkan situasi yang bersesuaian dengan syarat batas kanan Robin: ux (L, t) = u(L, t). Uraian lanjut mengenai arti fisis syarat batas Robin dapat dilihat pada Strauss Section 4.3. [Strauss 1.4 no 5] Dua batang penghantar berpenampang sama...

4.1

Prinsip Maksimum

Teorema 4.1.1. Prinsip Maksimum: Solusi persamaan difusi (5.0.1) pada ’persegi panjang’ [0, L] × [0, ∞) memenuhi u(x, t) ≤ max u pada batas ’persegi panjang’,

Dengan kata lain, solusi persamaan difusi u(x, t) mencapai ekstrim pada saat awal t = 0 atau pada batas kiri x = 0 atau kanan x = L: u(x, t) ≤ max {u(x, 0), u(0, t), u(L, t)} .

Prinsip Maksimum versi kuat: Solusi persamaan difusi (5.0.1) pada persegi panjang [0, L]× [0, ∞) memenuhi u(x, t) < max {u(x, 0), u(0, t), u(L, t)} , atau u(x, t) = konstan. Pikirkan kesesuaian prinsip maksimum pada masalah distribusi tinta. Bukti prinsip maksimum: Strauss hal 41.


28

Coba buktikan prinsip minimum berikut. Jika u(x, t) solusi persamaan difusi pada persegi panjang [0, L] × [0, ∞) maka berlaku min {u(x, 0), u(0, t), u(L, t)} ≤ u(x, t). Gabungan antara versi kuat prinsip maksimum dan mininum adalah: solusi u(x, t) dari persamaan difusi memenuhi

SRP '12

min {u(x, 0), u(0, t), u(L, t)} < u(x, t) < max {u(x, 0), u(0, t), u(L, t)} , atau u(x, t) = konstan. Ketunggalan Selanjutnya kita akan membuktikan ketunggalan solusi persamaan difusi dengan syarat batas tipe Dirichlet berikut:   ut − kuxx = f (x, t), 0 < x < L, t > 0 u(x, 0) = ϕ(x) (4.1.1)  u(0, t) = g(t), u(L, t) = h(t)

Bukti: Misalkan u1 (x, t) dan u2 (x, t) merupakan dua solusi (4.1.1) akan ditunjukkan bahwa u1 (x, t) = u2 (x, t). Tuliskan w(x, t) ≡ u1 (x, t) − u2 (x, t), maka w(x, t) memenuhi   wt − kwxx = 0, 0 < x < L, t > 0 u(x, 0) = 0 (4.1.2)  w(0, t) = 0, w(L, t) = 0 Terapkan prinsip maksimum dan minimum pada persamaan difusi dengan syarat batas dan syarat awal di atas akan menghasilkan w(x, t) = 0 atau u1 (x, t) = u2 (x, t).

Alternatif: metoda energi untuk membuktikan ketunggalan solusi masalah difusi. Misalkan u1 (x, t) dan u2 (x, t) merupakan dua solusi (4.1.1), misalkan w(x, t) ≡ u1 (x, t) − u2 (x, t), maka w(x, t) memenuhi (4.1.2). 0 = 0 · w = (wt − kwxx )w = (1/2w2 )t + (−kwx w)x + kwx2 . Integralkan terhadap x dengan batas 0 dan L, maka ∫ L ∫ L 2 x=L 0= (1/2w )t dx − kwx w|x=0 + k wx2 dx 0

0

atau ∫ ∫ L d L1 2 w dx = −k wx2 dx. dt 0 2 0 ∫L 2 Berarti 0 w dx turun untuk setiap t ≥ 0 ∫

L

∫ w dx ≤ 2

0

Sedangkan

L

w(x, 0)2 dx, 0

∫L 0

w(x, 0)2 dx = 0, jadi w ≡ 0 atau u1 (x, t) = u2 (x, t).

Interpretasi syarat batas Robin bagi persamaan difusi (soal Strauss 2.3 no 8). Perhatikan persamaan difusi ut = kuxx pada domain 0 < x < L, t > 0 dengan syarat batas Robin


29

ux (0, t) − a0 u(0, t) = 0, dan ux (L, t) − aL u(L, t) = 0. Jika a0 > 0 dan aL > 0, gunakan metoda ∫ L 2 energi untuk membuktikan bahwa kedua titik ujung berperan dalam berkurangnya 0 u (x, t) dx. Dengan kata lain, sebagian dari energi hilang melalui batas, dan syarat batas ini dikenal sebagai syarat batas ’radiating’ atau ’dissipative’. Soal: Buktikan ketunggalan solusi persamaan difusi non-linear berikut ut = uxx − u3 , x ∈ [0, L], t > 0 u(0, t) = u(L, t) = 0 u(x, 0) = f (x)

SRP '12

Kestabilan Misalkan ui (x, t) untuk i = 1, 2 masing-masing nerupakan solusi (4.1.1) dengan syarat awal u(x, 0) = ϕi (x, t) untuk i = 1, 2. Akan ditunjukkan bahwa untuk setiap ε > 0, terdapat δ > 0 sehingga |ϕ1 − ϕ2 | < δ berakibat |u1 − u2 | < ε. Bukti: Tuliskan w(x, t) ≡ u1 (x, t) − u2 (x, t), maka w(x, t) memenuhi persamaan difusi   ut − kuxx = 0, 0 < x < L, t > 0 u(x, 0) = ϕ1 (x) − ϕ2 (x)  u(0, t) = 0, u(L, t) = 0 Penerapan prinsip maksimum menghasilkan w(x, t) ≤ max{ϕ1 − ϕ2 , 0} ≤ maxx |ϕ1 − ϕ2 |, sedangkan penerapan prinsip minimum menghasilkan min{ϕ1 − ϕ2 , 0} ≤ w(x, t). Lebih jauh lagi dapat diperoleh hubungan

−max|ϕ1 − ϕ2 | ≤ −max{ϕ2 − ϕ1 , 0} = min{ϕ1 − ϕ2 , 0} ≤ w(x, t). Sehingga akhirnya diperoleh −maxx |ϕ1 − ϕ2 | ≤ w ≤ maxx |ϕ1 − ϕ2 | atau |w| = |u1 − u2 | ≤ maxx |ϕ1 − ϕ2 |.

Tugas: Buktikan kestabilan solusi persamaan difusi menggunakan metoda energi.

4.2

Fungsi Green

Solusi persamaan difusi (5.0.1) untuk x ∈ R dengan syarat awal u(x, 0) = ϕ(x) dapat direpresentasikan secara eksplisit: ∫ ∞ u(x, t) = S(x − y, t)ϕ(y)dy, (4.2.1) −∞

dengan fungsi Green S(x, t) adalah 1 2 S(x, t) = √ e−x /4kt , t > 0. 4πkt

(4.2.2)

Fungsi S(x, t) dikenal dengan nama fungsi Green. Nama lainnya: Gaussian, propagator, fundamental solution. Tampak dari Gambar 1.6 bahwa kurva S(x, t) menurun menuju nol seiring dengan bertambahnya


30 S

t kecil

t sedang t besar

Catatan:

SRP '12

Dengan demikian kita memiliki rumusan eksplisit bagi solusi persamaan difusi dengan syarat awal u(x, 0) = ϕ(x), yaitu ∫ ∞ 1 2 e−(x−y) /4kt ϕ(y)dy. (4.2.3) u(x, t) = √ 4πkt −∞

• S(x, t) merupakan solusi persamaan (5.0.1) dengan syarat awal ϕ(x) = δ(x) fungsi delta Dirac. ∫∞ • −∞ S(x, t)dx = 1, • limt→0 S(x, t) = 0. Soal:

1. (Smoothing processes in the solution of diffusion equation) Tentukan { solusi per0, x < 0 samaan difusi ut = kuxx , x ∈ R, t > 0 dengan syarat awal u(x, 0) = 1, x > 0 Simulasikan menggunakan Maple solusi u(x, t) untuk menunjukkan proses smoothing yang berasal dari persamaan difusi. Cermati laju dari proses smoothing tersebut untuk berbagai pilihan k, misalnya 1/10, 1/2, 1, 10. 2. Perhatikan bahwa ∂u ∂x merupakan solusi persamaan difusi dengan syarat awal δ(x). Solusi inilah yang dikenal sebagai Fungsi Green: S(x, t) ≡

4.3

∂u 1 2 = √ e−x /4kt ∂x 2 πkt

Metoda separasi variable, recall

Solusi persamaan difusi (5.0.1) pada selang berhingga [0, L] dapat diperoleh melalui metode separasi variabel. Pelajari metoda separasi variabel persamaan difusi. Solusinya berupa jumlahan suku-suku eiax e−ka t = e−ka t (cos ax + i sin ax). 2

2

Suku e−ka t menjelaskan bahwa efek damping membesar seiring dengan bertambahnya waktu. Selain itu, untuk bilangan gelombang a yang lebih besar faktor damping exp(−ka2 t) juga lebih besar. Dengan demikian jelas bahwa persamaan difusi merupakan persamaan pengatur bagi masalah-masalah difusi atau damping. 2

Metode beda hingga bagi persamaan difusi

SRP '12

Bab 5

Pada subbab ini akan diuraikan berbagai metode beda hingga untuk persamaan difusi ut = duxx ,

(5.0.1)

dengan d suatu konstanta.

5.1

Metode FTCS (Forward Time Center Space)

Metode FTCS sering disebut dengan metode eksplisit bagi persamaan difusi. Pada metode ini, forward time diterapkan pada ut dengan akurasi O(∆t) dan metode beda pusat yang diterapkan pada uxx dengan akurasi O(∆x2 ), sehingga diperoleh persamaan beda berikut: un+1 = unj + d j

∆t (un − 2unj + unj−1 ), (∆x)2 j+1

(5.1.1)

atau un+1 = (1 − 2S)unj + S(unj+1 + unj−1 ), dengan S = d j

∆t . (∆x)2

j, n+1

j-1, n

j, n

j+1, n

Gambar 5.1.1: Stencil FTCS untuk persamaan difusi. Kestabilan Substitusikan unj = ρn eia∆xj ke dalam (5.1.2), sehingga diperoleh:

ρ = (1 − 2S) + S(eia∆x + e−ia∆x ) = 1 + 2S(cos a∆x − 1). Agar skema stabil, maka |ρ| ≤ 1, yaitu |ρ| = |1 + 2S(cos a∆x − 1)| −1 ≤ 1 + 2S(cos a∆x − 1) −1 ≤ S(cos a∆x − 1) 0≤ S(1 − cos ( a∆x) ) 0≤ 2S sin2 a∆x 2 31

≤ 1, ≤ 1, ≤ 0, ≤ 1, ≤ 1,

(5.1.2)

BAB 5. METODE BEDA HINGGA BAGI PERSAMAAN DIFUSI sehingga 2S ≤ 1. Jadi, skema akan stabil jika S = d

32

1 ∆t ≤ . 2 (∆x) 2

Konsistensi Diberikan dua hampiran berikut: 1 1 un+1 = unj + ∆t ut |nj + ∆t2 utt |nj + O(∆t3 ) + ∆t4 u4t |nj + .... j 2 2

(5.1.4)

SRP '12

1 1 unj±1 = unj ± ∆x ux |nj + ∆x2 uxx |nj ± O(∆x3 ) + ∆x4 u4x |nj + .... 2 2 Substitusikan (5.1.3) dan (5.1.4) ke dalam (5.1.1), ( ) 1 2 2 n n n n 2 ut |j + ∆t utt |j + .... = d uxx |j + ∆x u4x |j + .... , 2 2! 4!

(5.1.3)

selanjutnya dapat disederhanakan menjadi ( )n 1 2 (ut − duxx )nj + ∆td2 u4x − ∆x2 du4x + .... = 0. 2 4! j Jelas bahwa skema konsisten. Selanjutnya suku pertama ( 2 d 1 2 2 ∆td u4x − ∆x du4x = d∆t − 2 4! 2

truncation term yaitu ) 1 2 ∆x u4x , 6

1 yang berupa suku difusi akan bernilai nol jika S = . 6

5.2

Metode Implisit BTCS (Backward Time Center Space)

Metode BTCS memiliki akurasi O(∆t, ∆x2 ). Persamaan beda untuk persamaan difusi dengan menggunakan metode BTCS adalah ( n+1 ) n+1 n+1 n u un+1 − 2u − u + u j j+1 j j j−1 =d . (5.2.1) ∆t (∆x)2

j, n+1 j+1, n+1

j-1, n+1

j, n

Gambar 5.2.1: Stencil BTCS untuk persamaan difusi. Kestabilan Substitusikan unj = ρn eia∆xj ke dalam (5.2.1) menghasilkan: ρ − 1 = Sρ(eia∆x + e−ia∆x − 2), ρ − 1 = 2Sρ(cos a∆x − 1), dengan S =

d∆t . Dengan demikian, (∆x)2 ρ=

1 1 = 1 + 2S(1 − cos a∆x) 1 + 4S sin2

a∆x 2

.

BAB 5. METODE BEDA HINGGA BAGI PERSAMAAN DIFUSI

33

Karena untuk setiap S dan a penyebut selalu lebih besar dari satu, maka jelas bahwa d∆t . |ρ| < 1 Jadi skema ini stabil ubtuk setiap S = (∆x)2 Konsistensi Substitusikan hampiran Taylor un+1 dan un+1 j j±1 ke dalam (5.2.1), sehingga menghasilkan: ut |nj

(∆t) + utt |nj + .... = d 2!

( uxx |nj

) (∆x2 ) n + u4x |j + .... , 4!

ut |nj +

SRP '12

yang juga dapat dituliskan sebagai berikut: ( ) d(∆x2 ) ( ) (∆t) 2 d u4x |nj +.... = d uxx |nj + d∆t u4x |nj + .... + u4x |nj + ∆t u4xt |nj + .... , 2! 4!

dan dapat disederhanakan lagi menjadi ( ) 1 2 1 2 (ut − duxx ) + − d ∆t − d∆x u4x + .... = 0. 2 12

Jelas bahwa persamaan beda (5.2.1) konsisten. Lebih jauh, suku pertama truncation term d∆t 1 1 akan bernilai nol jika − 12 d2 ∆t − 12 d∆x2 = 0, atau ketika S = =− . 2 (∆x) 6 Selanjutnya persamaan beda (5.2.1) dapat dituliskan sebagai berikut : n+1 n (1 + 2S)un+1 − Sun+1 j j−1 − Suj+1 = uj .

(5.2.2)

Misalkan syarat batasnya bertipe Dirichlet yaitu u(0, t) = f1 dan u(L, t) = f2 . Titik-titik yang harus dihitung adalah un+1 , dengan j = 1, 2, ...., Nx − 1, j n = 0, 1, ...., Nt − 1.

Jika persamaan (5.2.2) dievaluasi untuk j = 1 hingga j = N x − 1 berturut-turut akan menghasilkan (1 + 2S)un+1 − Sun+1 = Su(0, t) + un1 , 1 2 n+1 n+1 n+1 −Su1 + (1 + 2S)u2 − Su3 = un2 , .. .. .. . . . n+1 n −Sun+1 + (1 + 2S)u Nx −2 Nx −1 = uNx −1 + Su(L, t) Persamaan-persamaan di atas dapat disusun membentuk SPL: A(Nx −1)×(Nx −1) un+1 (Nx −1)×1 = b(Nx −1)×1 , dengan

     A=    

(1 + 2S) −S 0 .. . 0

−S

···

···

(1 + 2S) −S ··· −S (1 + 2S)

0 .. .

. ···

···

. −S −S (1 + 2S)

     ,    

BAB 5. METODE BEDA HINGGA BAGI PERSAMAAN DIFUSI 

dan

b(Nx −1)×1

   =   

un1 + Sf1 .. . unj .. .

34

       

unNx −1 + Sf2

5.3

Metode Crank - Nicholson

n+ 1 ut |j 2

=

un+1 − unj j ∆t

+ O(∆t2 ),

(5.3.1)

Sedangkan uxx di titik grid j, (n + 21 ) dihampiri dengan n+ 21

uxx |j

SRP '12

Perhatikan persamaan difusi (5.0.1). Terapkan center time untuk menghampiri ut di titik grid j, (n + 21 )

1 + uxx |nj ), = (uxx |n+1 j 2

(5.3.2)

dan selanjutnya uxx |n+1 dan uxx |nj ) dihampiri dengan center difference. Sehingga diperj oleh persamaan beda untuk metoda Crank-Nicholson berikut ( n+1 ) n+1 n+1 un+1 − unj unj+1 − 2unj + unj−1 d uj+1 − 2uj + uj−1 j = . (5.3.3) + ∆t 2 (∆x)2 (∆x)2 Dengan demikian metode Crank - Nicholson memiliki akurasi O(∆t2 , ∆x2 ).

j, n+1 j+1, n+1

j-1, n+1 j, n+1/2 j-1, n

j+1, n j, n

Gambar 5.3.1: Stencil Crank - Nicholson untuk persamaan difusi.

Kestabilan Substitusikan unj = ρn eia∆xj ke dalam (5.3.3) dan dengan sedikit manipulasi aljabar dapat diperoleh: ρ = ρ =

1 + S(cos a∆x − 1) , 1 − S(cos a∆x − 1) ( ) 1 − 2S sin2 a∆x 2 ( ), 1 + 2S sin2 a∆x 2

(5.3.4)

(5.3.5)

d∆t . Dari (5.3.5) jelas bahwa |ρ| ≤ 1 untuk setiap S > 0. Jadi, skema (∆x)2 stabil tanpa syarat. dengan S =


35

yang dapat disederhanakan menjadi (ut − duxx ) −

d ∆x2 u4x + ...u6x + .... = 0. 12

SRP '12

Berikut ini adalah persamaan differensial termodifikasi dari persamaan difusi dengan penerapan metode Crank-Nicholson. ( d 2 1 n n ut |j + ∆t utt |j + .... = uxx |n+1 + ∆x2 u4x |n+1 + .... j j 2 2 4! ) 2 + uxx |nj + ∆x2 u4x |nj + .... , 4! ( ) n 1 d∆t 1 n n 2 u4x + .... , ut |j + ∆t utt |j + .... = d uxx + ∆x u4x + 2 4! 2 j

Dari persamaan diferensial termodifikasi di atas jelas bahwa metode Crank-Nicholson konsisten, dan suku pertama truncation errornya berupa suku difusi. Algoritma metoda implisit bagi persamaan difusi ut = kuxx , 0 < x < L, t > 0, dengan syarat batas % (Example on) constructing the coefficient matrix A, size n × n, tridiagonal. (The size of A depends on the type of left and right boundary, Dirichlet or Neumann. D = sparse(1:n,1:n,(1+2*s)*ones(1,n),n,n); E = sparse(2:n,1:n-1,-s*ones(1,n-1),n,n); A=full(E+D+E’) % Inserting initial condition u(x, 0) % Computing u(x, tn+1 ) using u(x, tn ) for n=1:Nt % Construct a vector RHS, values of u(x, tn ) needed. (Left and right boundary condition should be incorporated.) % Solving a tridiagonal system Au(x, tn ) = RHS end

5.4

Metode - θ

Metode - θ merupakan modifikasi dari metode Crank-Nicholson. Persamaan beda untuk persamaan difusi dengan metode ini adalah ( n+1 ) uj+1 − 2un+1 + un+1 un+1 − unj unj+1 − 2unj + unj−1 j j−1 j + (1 − θ) =d θ . (5.4.1) ∆t (∆x)2 (∆x)2 Perhatikan persamaan (5.4.1). Jika θ = 0, maka persamaan tersebut menjadi persamaan beda untuk metode FTCS dengan akurasi O(∆t, ∆x2 ). Jika θ = 12 , maka persamaan tersebut menjadi persamaan beda untuk metode Crank-Nicholson dengan akurasi O(∆t2 , ∆x2 ). Dapat dikatakan bahwa metode-θ terdiri dari (1 − θ) metoda eksplisit, dan θ metoda implisit. Kestabilan Substitusikan unj = ρn eia∆xj ke dalam (5.4.1) dan dengan sedikit manipulasi aljabar dapat


j, n+1 è

j-1, n+1

j-1, n

j, n

36

j+1, n+1

j+1, n

diperoleh amplification factor: ρ=

1 − 2S + 2θS(1 − cos a∆x) , 1 + 2θS(1 − cos a∆x)

dengan S =

d∆t . Perhatikan bahwa |ρ| ≤ 1 , untuk setiap S. Jadi, skema stabil tanpa (∆x)2

syarat.

5.5

SRP '12

Gambar 5.4.1: Stencil metode - θ untuk persamaan difusi.

Soal-soal

1. [PR] Diketahui persamaan difusi ut = uxx

0<x<4

u(0, t) = 4,

ux (4, t) = 0

u(x, 0) = (x + 1)(4 − x)

(a) Tentukan solusi eksaknya dengan metoda separasi variable.

(b) Implementasikan metode FTCS. Plot hasilnya unumerik (x, ti ) bersama-sama dengan ueksak (x, ti ) untuk waktu-waktu ti tertentu.

(c) Implementasikan metode BTCS (implisit). Plot hasilnya unumerik (x, ti ) bersamasama dengan ueksak (x, ti ) untuk waktu-waktu ti tertentu. 2. (a) Carilah solusi u(x, t) dari persamaan difusi ut = Duxx , 0 < x < 1, t > 0, dengan syarat awal u(x, 0) = sin πx dan syarat batas u(0, t) = u(1, t) = 0, t > 0. (Petunjuk : Coba solusi berbentuk f (t) sin πx.) ∫1 (b) Definisikan q(t) ≡ 0 u(x, t)dx, bagaimanakah q(t) berubah terhadap waktu?

3. [Penney 9.5 no. 17] (Steady-state temperature) Misalkan suatu batang penghantar yang terisolasi secara lateral mempunyai suhu awal u(x, 0) = f (x), dan suhu kedua ujungnya u(0, t) = A dan u(L, t) = B, dengan A, B konstanta tak nol. (a) Teramati bahwa setelah waktu yang cukup lama t → ∞, u(x, t) mendekati steady-state temperature uss (x) yang tak lain adalah solusi equilibrium dari persamaan panas serta syarat batas ∂u ∂2u = k 2, ∂t ∂x Tentukan uss (x).

u(0, t) = A, u(L, t) = B.


37

(b) Selanjutnya didefinisikan transient temperature utr (x, t) ≡ u(x, t) − uss (x). Tunjukkan bahwa utr (x, t) memenuhi persamaan panas dengan syarat batas homogen u(0, t) = 0, u(L, t) = 0 dan syarat awal utr (x, 0) = g(x) ≡ f (x) − uss (x). (c) Pada akhirnya tentukan solusi masalah semula.

SRP '12

Diketahui sebuah dinding logam dengan tebal L dan koefisien penghantar panas k. Misalkan u(x, t) menyatakan suhu suatu titik di posisi x pada dinding (Suhu dinding logam kita anggap 4. homogen pada kedua arah lainnya). Misalkan tepi dinding kiri bersentuhan dengan suhu udara luar, sedangkan tepi kanan dinding bersentuhan dengan suhu kamar, lihat Gambar.

Suhu udara luar = B

Suhu kamar = A

0

L

(a) Jika pada saat awal suhu dinding adalah T0 , maka pdp beserta syarat batas yang sesuai adalah ut = kuxx , u(0, t) = B,

pada

0 < x < L,

t>0

u(L, t) = A.

(∗)

Syarat batas (*) disebut syarat batas tipe Dirichlet tak homogen. Jika A dan B bernilai nol, maka disebut syarat batas Dirichlet homogen.

(b) Jika diketahui A = B, kita tahu bahwa masalah distribusi suhu dinding akan simetris terhadap titik tengah dinding, yaitu x = L/2. Jadi di sini kita bisa memformulasikan masalah distribusi suhu dinding cukup pada domain separo dinding saja, menjadi ut = kuxx ,

pada

u(0, t) = A, ux (L/2, t) = 0.

0 < x < L/2,

t>0

(∗∗)

Perhatikan perubahan syarat batas di x = L/2! Syarat batas (**) disebut syarat batas tipe Neumann homogen.

x

Bab 6

SRP '12

Persamaan Laplace & Poisson Perhatikan masalah distribusi panas pada lempeng persegi [0, 1]×[0, 1]. Misalkan u(x, y, t) menyatakan suhu lempeng di posisi (x, y) saat t. Setelah waktu yang cukup lama, suhu lempeng telah mencapai kondisi equilibrium, sehingga sudah tidak bergantung lagi pada waktu t. Dalam kondisi equilibrium, suhu lempeng hanya bergantung pada kondisi batas. Suhu lempeng saat kondisi equilibrium dinyatakan sebagai fungsi dari x dan y saja, misal u(x, y). Misalkan tepi kiri kanan lempeng terisolasi, sehingga syarat batasnya ux (0, y) = ux (1, y) = 0. Sepanjang tepi bawah lempeng suhu dipertahankan nol, jadi syarat batasnya u(x, 0) = 0. Sepanjang tepi atas lempeng, suhu dipertahankan 1000 , jadi syarat batasnya u(x, 1) = 1000 . Persamaan diferensial parsial beserta syarat batas yang sesuai untuk masalah ini adalah

△u = 0

pada

[0, 1] × [0, 1]

ux (0, y) = ux (1, y) = 0 u(x, 0) = 0, u(x, 1) = 1000 Melalui substitusi langsung, tunjukkan bahwa u(x, y) = 100y merupakan solusi dari masalah di atas. Plot dari solusi u(x, y) diberikan pada gambar di samping. Periksa kesesuaian antara gambar permukaan dengan keempat syarat batasnya.

Dengan demikian dapat kita simpulkan bahwa persamaan Laplace △u = 0 adalah persamaan bagi solusi steady persamaan difusi ut = △u. Bentuk umum dari persamaan Laplace adalah △u = 0

pada

Ω,

(6.0.1)

dengan △ = ∂xx + ∂yy adalah operator Laplace. Sedangkan △u = f

pada

Ω,

(6.0.2)

dengan f fungsi sebarang, dikenal sebagai persamaan Poisson. Solusi dari persamaan Laplace disebut sebagai fungsi harmonik. Buktikan bahwa fungsi-fungsi berikut adalah fungsi harmonik pada domain R2 : 1. Semua polinom berderajat satu: f (x, y) = ax + by + c, a, b, c ∈ R 38

BAB 6. PERSAMAAN LAPLACE & POISSON

39

2. Polinom berderajat dua: f (x, y) = a(x2 − y 2 ) + b, a, b ∈ R 3. Polinom berderajat tiga: f (x, y) = a(x3 − 3xy 2 ) + d(y 3 − 3x2 y), a, d ∈ R

SRP '12

4. Fungsi f (x, y) = sinh nx sin ny, n ∈ R. (Catatan: Fungsi harmonik lainnya cosh nx sin ny, sinh nx cos ny, cosh nx cos ny. 2 2 Perhatikan bahwa fungsi-fungsi e−kn t cos nx, e−kn t sin nx, n ∈ R merupakan solusi persamaan difusi ut = kuxx . Suku-suku tersebut muncul pada solusi dengan metoda separasi variabel. Hal serupa juga berlaku bagi fungsi-fungsi cos cnt cos nx, sin cnt cos nx, cos cnt sin nx, sin cnt sin nx, n ∈ R yang merupakan solusi persamaan gelombang utt = c2 uxx . Persamaan Laplace dan Poisson ditemui pada masalah-masalah seperti • Masalah elektrostatik • Aliran fluida steady • Fungsi kompleks yang analitik • Gerak Brown

Soal [Strauss 1.5 no. 4] Perhatikan persamaan Poisson dengan syarat batas tipe Neumann homogen { △u = f pada D (6.0.3) ∂u ∂n = 0 pada ∂D a. Apa yang dapat kita tambahkan pada solusi sehingga menghasilkan solusi lain? (Ketunggalan solusi tidak dipenuhi) b. Interpretasikan hasil a. untuk masalah aliran panas dan difusi tinta.

6.1

Prinsip maksimum dan ketunggalan solusi persamaan Laplace

Prinsip Maksimum Misalkan D adalah himpunan buka, tersambung, dan terbatas. Fungsi harmonik u di domain D yang kontinu di D(= D ∪ ∂D) memenuhi: Versi lemah: min u(x, y) pada ∂D ≤ u(x, y) ≤ max u(x, y) pada ∂D, Bukti: Strauss 6.1 Versi Kuat:

min u(x, y) pada ∂D < u(x, y) < max u(x, y) pada ∂D,kecuali jika u(x, y) = konstan. Bukti: Menggunakan mean value property Strauss 6.3 Perhatikan persamaan Poisson dengan syarat batas tipe Dirichlet berikut { △u = f pada D u = h pada ∂D

(6.1.1)


40

Akan ditunjukkan bahwa solusi (6.1.2) adalah tunggal. Misalkan u1 dan u2 adalah masingmasing solusi (6.1.2). Tuliskan w ≡ u1 − u2 , maka w memenuhi persamaan Laplace pada domain D dengan syarat batas tipe Dirichlet homogen: △w = 0 pada D,

w = 0 pada ∂D.

Menurut prinsip maksimum 0 ≤ w(x) ≤ 0

SRP '12

untuk setiap titik x ∈ D. Ini berarti w(x) = 0 atau u1 (x) = u2 (x). Soal: [Strauss 6.1 no. 11] Tunjukkan bahwa persamaan Poisson dengan syarat batas tipe Neumann berikut { △u = f pada D (6.1.2) ∂u ∂n = g pada ∂D ∫∫∫ ∫∫ punya solusi jika dan hanya jika D f dV = ∂D g dS. Petunjuk terapkan teorema divergensi Gauss pada △u = f .

Aplikasi pada aliran fluida

Tinjau suatu aliran fluida steady pada belokan seperti pada gambar. Perhatikan lebih detail domain persegi panjang tepat pada sudut belokan. Dalam keadaan setimbang, total massa fluida yang masuk menembus sisi AB akan sama dengan total massa fluida yang menembus sisi BC, sehingga pada domain ABCD berlaku persamaan konservasi massa ∇ · u = 0.

Jika aliran fluida bersifat irrotasional (∇ × u = 0, berarti medan vektor u bersifat konservatif. Bagi medan vektor konservatif, dijamin terdapat fungsi skalar ϕ(x, y) sehingga u = ∇ϕ. Substitusikan hubungan ini ke dalam persamaan konservasi massa akan menghasilkan ∇ · ∇ϕ = ∆ϕ = 0

pada

persegi ABCD.

Kondisi yang harus menyertai keberlakuan persamaan konservasi massa adalah fluks fluida menembus AB = fluks fluida menembus BC Mengingat fluks adalah u · n× panjang sisi, maka syarat keberlakuan konservasi massa di atas menjadi

D

∇ϕ · n AB ≡ ∇ϕ · n BC. Syarat ini tepat merupakan syarat eksistensi solusi (6.1.2) bagi persamaan Laplace ∆ϕ = 0 pada domain persegi panjang ABCD. Coba yakinilah hal ini.

A

C

B


41

Sifat invarian di 2-D Buktikan bahwa operator Laplace ∂xx +∂yy bersifat invarian terhadap translasi dan rotasi. Dengan demikian Laplacian △ merupakan model bagi situasi fisik yang bersifat isotropic, yang artinya fungsi harmonik tidak memilih arah tertentu. (Bandingkan dengan solusi persamaan transport aux + but = 0 yang bertranslasi pada arah (a, b)T .) Sifat invarian terhadap rotasi mengilhami penulisan △ dalam koordinat polar. Uraian mengenai persamaan Laplace dalam koordinat polar akan disajikan pada Subbab6.3.

SRP '12

Soal Latihan 6.1.1. Perhatikan persamaan Laplace ∆u = 0 (uxx + uyy = 0). Cari semua polinom homogen derajat dua berbentuk p(x, y) = ax2 + bxy + cy 2 yang memenuhi persamaan Laplace! Sketsakan solusi polinom derajat dua yang relatif sederhana.

Cari solusi dari △u = 0 pada persegi satuan [0, 1] × [0, 1] yang memenuhi syarat batas u(x, 1) = sin πx dan u = 0 pada ketiga sisi lainnya. Mengingat syarat batasnya, cobalah solusi berbentuk u(x, y) = a(y) sin πx. Sketsa solusinya diberikan pada gambar di samping. Bandingkan dengan solusi yang diperoleh melalui metoda separasi variabel Kerjakan soal-soal Strauss 6.1 no 1-12.

6.2

Persamaan Laplace pada domain persegi panjang

Diskusikan solusi persamaan Laplace △u = 0 pada domain D = [0, a]×[0, b] dengan syarat batas seperti pada gambar y

b

u = f ( y)

u= g(x)

Du = 0

ux = k ( y)

u y + u = h( x )

a

x

Kerjakan soal-soal Strauss 6.2 no 1-7. Soal-soal: 1. Diketahui persamaan Laplace beserta syarat batasnya uxx + uyy = 0, ux (0, y) = 0,

pada pers. panj. [0, 2] × [0, 1]

ux (2, y) = −1,

uy (x, 0) = 0,

uy (x, 1) = 1/2,

(a) Periksa bahwa syarat batas Neumann yang diberikan memenuhi 0. Ini berarti ’syarat perlu’ eksistensi solusi dipenuhi.

∫

∂u batas pp ∂n

ds =

(b) Tentukan solusi u(x, y) berbentuk polinom berderajat-2 yang memenuhi persamaan Laplace beserta syarat batasnya.


42

2. Perhatikan persamaan Laplace uxx + uyy = 0 dengan syarat batas ux (1, y) = −1 uy (x, 0) = 0, uy (x, 1) = 1. ∫ ∂u (a) Periksa keberlakuan batas pp ∂n ds = 0, yang mana berarti ’syarat perlu’ eksistensi solusi dipenuhi. ux (0, y) = 0,

(b) Carilah solusi eksak masalah di atas yang berbentuk polinom derajat dua

6.3

Rumus Poisson

SRP '12

Kerjakan soal-soal Strauss 6.1 no 1 - 12

Perhatikan masalah temperatur steady pada lempeng lingkaran berjari-jari a yang diatur oleh persamaan Laplace dalam koordinat polar berikut. 1 1 △u = ∂r2 u + ∂r u + 2 ∂θ2 u = 0, r r u(a, θ) = h(θ).

pada (r, θ) ∈ (0, a) × (0, 2π)

(6.3.1)

(6.3.2)

Perhatikan bahwa di sini h adalah fungsi berperiode 2π. Selain itu mengingat (r, θ) dan (r, θ + 2π) adalah dua titik yang sama, maka di sini berlaku syarat batas natural u(r, θ) = u(r, θ + 2π).

(6.3.3)

Seperti biasa, metoda pemisahan variabel diawali dengan memisalkan u(r, θ) = R(r)Θ(θ).

(6.3.4)

Substitusikan (6.3.4) ke dalam (6.3.1) menghasilkan 1 1 R′′ Θ + R′ Θ + 2 RΘ′′ = 0, r r jika dibagi dengan RΘ/r2 akan diperoleh r2 R′′ + rR′ Θ′′ + = 0. R Θ Selanjutnya kita peroleh Θ′′ r2 R′′ + rR′ =− =λ R Θ

untuk suatu konstanta λ. Dengan demikian kita peroleh dua buah persamaan diferensial biasa r2 R′′ + rR′ − λR = 0

(6.3.5)

Θ′′ + λΘ = 0

(6.3.6)

Tunjukkan bahwa solusi (6.3.6) yang memenuhi (6.3.3) hanya menghasilkan pasang nilai eigen dan fungsi eigen berikut λ0 = 0, 2

λn = n ,

Θ0 (θ) = 1

(6.3.7)

Θn (θ) = An cos nθ + Bn sin nθ,

(6.3.8)


43

untuk n = 1, 2, · · · . Untuk λ0 = 0 persamaan (6.3.5) menghasilkan solusi R0 (r) = C0 + D0 ln r

(6.3.9)

Jika kita menginginkan solusi yang kontinu di r = 0, maka fungsi ln r tidak terpakai, sehingga R0 (r) = 1. Untuk λn = n2 persamaan (6.3.5) menghasilkan solusi Rn (r) = Cn rn + Dn r−n

(6.3.10)

Solusi yang kontinu di r = 0 hanyalah rn , sehingga Rn (r) = rn . Jadi telah kita dapat solusi

SRP '12

un (r, θ) = Rn (r)Θ(θ) = rn (An cos nθ + Bn sin nθ),

untuk n = 0, 1, 2, · · · . Pandang kombinasi linier dari semua bentuk di atas A0 n + Σ∞ n=1 r (An cos nθ + Bn sin nθ). 2 Gunakan syarat batas (6.3.2) menghasilkan u(r, θ) =

A0 n + Σ∞ n=1 a (An cos nθ + Bn sin nθ). 2 Jadi An , Bn adalah koefisien deret Fourier dari h(θ), suatu fungsi berperiode 2π ∫ 2π 1 h(ϕ) cos nϕ dϕ (6.3.11) An = πan 0 ∫ 2π 1 Bn = h(ϕ) sin nϕ dϕ (6.3.12) πan 0 h(θ) =

untuk n = 0, 1, 2, · · · . Sehingga solusi (6.3.1,6.3.2) adalah ∫ 2π ∫ 2π 1 rn h(ϕ) dϕ + Σ∞ h(ϕ)(cos nϕ cos nθ + sin nϕ sin nθ) u(r, θ) = n=1 2π πan 0 0 { } ∫ 2π ( r )n 1 ∞ = h(ϕ) 1 + 2Σn=1 cos n(θ − ϕ) dϕ a 2π 0 Perhatikan bahwa ( r )n 1 1 + 2Σ∞ cos n(θ − ϕ) = real(∗) n=1 a 2π

a2 − r2 a2 − 2ar cos(θ − ϕ) + r2 Sehingga pada akhirnya kita peroleh rumus Poisson berikut. ∫ a2 − r2 2π h(ϕ) u(r, θ) = dϕ, untuk r ≤ a. 2 2π a − 2ar cos(θ − ϕ) + r2 0 (∗) =

(6.3.13)

Perhatikan bahwa rumus Poisson menyajikan solusi persamaan Laplace (6.3.1) dengan syarat batas tipe Dirichlet (6.3.2) secara eksplisit. Sayangnya rumus di atas tidak ada hubungannya dengan persamaan Poisson. Berikut ini rumus Poisson dituliskan secara geometris: a2 − |x|2 u(x) = 2πa

∫ |x′ |=a

dengan x = (r, θ), x′ = (a, ϕ)

u(x′ ) ds′ , |x − x′ |2

(6.3.14)


44

Mean value property Misalkan h(ϕ) merupakan fungsi kontinu pada lingkaran C = ∂D, maka rumus Poisson memberikan fungsi harmonik (satu-satunya) pada D yang memenuhi lim u(x) = h(x0 ), untuk setiap x0 ∈ ∂D

x→x0

SRP '12

Ini berarti u(x) merupakan fungsi kontinu pada D = D ∪ ∂D. Lebih jauh, u(x) diferensiabel tak hingga kali di D. Teorema 6.3.1. (Mean value property) Misalkan u(x) fungsi harmonik di D yang kontinu di D, maka nilai dari u di pusat D sama dengan nilai rata-rata dari u pada ∂D. Proof: Poisson formula dihitung pada x = 0 menghasilkan ∫ ∫ a2 u(x′ ) ′ 1 u(0) = ds = u(x′ ) ds′ , 2πa |x′ |=a a2 2πa |x′ |=a

yang tak lain adalah nilai rata-rata dari u di keliling lingkaran |x′ | = a.

Dengan menggunakan teorema nilai rata-rata inilah prinsip maksimum versi kuat dibuktikan. Kerjakan soal-soal Strauss 6.3 Soal Latihan 6.3.2.

1. Tunjukkan bahwa dalam koordinat polar

1 1 △ = ∂r2 + ∂r + 2 ∂θ2 . r r

2. Tunjukkan bahwa u(r, θ) = a ln |r| + b, dengan a, b konstanta sebarang, adalah solusi dari △u = 0. Fungsi ln |r| dan fungsi konstan yang merupakan solusi dikenal sebagai solusi fundamental. Kedua solusi ini akan muncul pada formulasi metoda separasi variable pada persamaan Laplace dalam koordinat polar. Perhatikan bahwa solusi ln |r| tak terdefinisi di r = 0, lihat gambar. 3. Tentukan solusi (6.3.1) dan (6.3.2) jika diketahui { T0 untuk 0 < θ < π h(θ) = T1 untuk π < θ < 2π

4. [Kreyszig 11.9 no. 13] (Temperatur pada lempeng setengah lingkaran.) Tentukan temperatur u pada lempeng {(r, θ)|r < R, 0 < θ < π}, dimana temperatur pada sepanjang tepi bawah adalah 00 dan pada batas setengah lingkaran atas adalah u0 , dengan u0 konstanta. 5. (Frekwensi fundamental genderang) Getaran membran pada genderang memenuhi persamaan 1 1 T utt = c2 (urr + ur + 2 uθθ ), c2 = , pada domain cakram lingkaran jari-jari R. r r ρ


45

Jika dimisalkan getaran membran simetris secara radial, maka u tak bergantung pada θ sehingga simpangan membran tiap saat dinyatakan sebagai u(r, t), dan persamaan diferensialnya 1 T utt = c2 (urr + ur ), c2 = , pada domain cakram lingkaran jari-jari R. r ρ

6.4

SRP '12

Jika genderang diikat erat, maka syarat batasnya u(R, t) = 0. Seringkali syarat awalnya bersifat radial simetri: u(r, 0) = h(r), ut (r, 0) = g(r). Penentuan solusi deret u(r, t) memerlukan formulasi menggunakan fungsi Bessel [Kreyszig 11.10].

Formula Green bagi solusi persamaan Laplace

Optional, sumber: Strauss Section 7 1. Misalkan diberikan syarat batas sebagai berikut: ∆u = 0 u = φ pada ∂Ω1 ∂u = ψ pada ∂Ω2 , ∂n maka solusinya adalah ∫ 2πu(x, y) =

φ(ξ) ∂Ω1

∂G (ξ, x)dξ + ∂n

∫

ψ(ξ)G(ξ, x)dξ. ∂Ω2

2. (Strauss hal. 181) Misalkan diberikan syarat batas sebagai berikut: ∆u = f pada D u = h pada ∂D, maka solusinya adalah ∫ ∫ u(x0 ) =

∂G(x, x0 ) dS + h(x) ∂n bdyD

∫ ∫ ∫

f (x)G(x, x0 )dx. D

Metode beda hingga bagi persamaan Laplace 2D 7.1

Metode langsung

SRP '12

Bab 7

Perhatikan persamaan Laplace uxx + uyy = 0 pada domain D. Terapkan metode beda pusat bagi uxx dan uyy , sehingga diperoleh ui+1,j − 2ui,j + ui−1,j ui,j+1 − 2ui,j + ui,j−1 + = 0. ∆x2 ∆y 2

-1

Pemilihan ∆x dan ∆y ternyata tidak terlalu berpengaruh sehingga dapat dipilih ∆x = ∆y, sehingga persamaan bedanya:

-1

-1

4ui,j − (ui+1,j + ui−1,j + ui,j+1 + ui,j−1 ) = 0.

4

(7.1.1)

-1

Baby example: Selesaikan persamaan Laplace pada domain persegi panjang dengan syarat awal tipe Dirichlet berikut: 0 0 0 0

0 ∗ ∗ 0

0 0 ∗ 24 ∗ 0 0 0

Terdapat 4 persamaan dengan 4 variabel sehingga solusinya dapat dicari sebagai solusi SPL biasa sehingga diperoleh 0 0 0 0 0 2 7 24 0 1 2 0 0 0 0 0 Cek bahwa (7.1.1) dipenuhi bagi keempat titik interiornya. Contoh 1: Perhatikan persamaan Laplace pada domain persegi D dengan syarat batas tipe Dirichlet seperti pada gambar. Partisikan persegi panjang tersebut dengan lebar selang partisi ∆x = ∆y seperti pada gambar. Perhatikan cara penomoran titik-titik partisinya. Solusi u dapat diperoleh sebagai solusi SPL Au = b, 46

BAB 7. METODE BEDA HINGGA BAGI PERSAMAAN LAPLACE 2D

47

dimana



u12×1

    =   

u1 u2 u3 u4 .. .





       

    =   

u12

b12×1

a+d a a a+b .. .

        

??

SRP '12

A12×12 =

Bagaimanakah A dan b jika terdapat Nx × Ny titik partisi? Contoh 2: Persamaan Laplace berlaku pada domain D dengan syarat batas tipe Dirichlet pada tepi bawah dan atas, syarat batas tipe Neumann pada tepi kiri dan kanan. Perhatikan gambar berikut. UN

UNXNY Ux = 0

UW

D

Ny buah

U2 U1

Up

UE

US

UNX

Perhatikan bahwa up , suatu titik partisi pada tepi kiri, ux |p = 0, sehingga uE − uW = 0 ⇔ uE = uW 2∆x . Persamaan beda menjadi: 4up − (uN + uS + 2uE ) = 0. • Coba rumuskan persamaan beda bagi titik pada tepi kanan! • Bagaimanakah bentuk matriks A? • Bagaimana jika syarat batas Neumann

∂u = Dn u = 0 berlaku pada seluruh ∂D? ∂n

BAB 7. METODE BEDA HINGGA BAGI PERSAMAAN LAPLACE 2D

7.2

48

Metode Iterasi

Keuntungan dari metode iterasi adalah membutuhkan memori komputer lebih sedikit. Selain itu, dengan metoda iterasi kesalahan formulasi lebih mudah dihindari, terutama dalam membangun matriks A yang mana sangat bergantung pada bentuk domain D dan syarat batas yang diberikan. Metode iterasi Jacobi (1) Metode ini diawali dengan sebuah nilai tebakan bagi ui,j , ∀i, j. Selanjutnya kita hitung (2)

(n+1)

ui,j

=

) 1 ( (n) (n) (n) (n) ui+1,j + ui−1,j + ui,j+1 + ui,j−1 . 4 (3)

SRP '12

ui,j menggunakan rumus rekursif berikut (7.2.1)

(n)

Dilanjutkan dengan menghitung ui,j , ∀i, j dan seterusnya untuk ui,j , ∀i, j untuk n = 4, 5, · · · . Dua kriteria yang biasa digunakan untuk menghentikan iterasi adalah (n+1)

max |ui,j i,j

(n)

− ui,j | < ε

atau n > Nmaksimum .

Perhatikan bahwa persamaan (7.2.1) tak lain adalah persamaan beda bagi persamaan difusi 2 − D : ut = uxx + uyy (center difference bagi uxx dan uyy serta forward difference bagi ut , dengan ∆x = ∆y, ∆t = 41 (∆x)2 ). Buktikan!

Metode Jacobi konvergen secara lambat. Metode iterasi konvergen jika matriks A memenuhi dominan diagonal: ∑ |aij | , i = 1, ...., n. |aii | ≥ j=1,j̸=i

Metode Gauss - Seidel Metode Gauss - Seidel merupakan perbaikan dari metode Jacobi. Misalkan dimulai dari titik apit kiri bawah bergerak ke kanan, lalu ke atas. Persamaan (7.2.1) tetap digunakan, (n) (n+1) hanya saja ui−1,j sebagai titik kirinya diganti dengan ui−1,j (karena toh sudah diketahui) (n)

(n+1)

dan ui,j−1 sebagai titik bawahnya diganti dengan ui,j−1 , sehingga menjadi (n+1)

ui,j

=

) 1 ( (n) (n+1) (n) (n+1) ui+1,j + ui−1,j + ui,j+1 + ui,j−1 . 4

Metode Gauss - Seidel dua kali lebih cepat daripada metode Jacobi. Kerjakan soal-soal Strauss 8.4

(7.2.2)

1 Pendahuluan pdp 2. 4 Persamaan Difusi Prinsip Maksimum Fungsi Green Metoda separasi variable, recall

Recommend Documents