Persamaan Difusi [PDF]

Citation preview

Persamaan Difusi Penurunan, Solusi Analitik, Solusi Numerik (Beda Hingga, RBF)



M. Jamhuri UIN Malang



April 7, 2013



M. Jamhuri



Persamaan Difusi



Penurunan Persamaan Difusi Misalkan u (x, t) menyatakan konsentrasi dari zat pada posisi x dan pada waktu t. Pada selang [x0 , x1 ] , massa zat ˆ x1 u (x, t) dx M= x0



dan perubahan massa



ˆ x1 dM ut (x, t) dx (1) = dt x0 Massa pada selang tersebut akan berubah bila ada zat yang masuk atau keluar selang tersebut. Hukum Fick mengatakan rata-rata penyebaran sebanding dengan gradien konsentrasi dM = zat masuk − zat keluar dt = kux (x1 , t) − kux (x0 , t) (2) dimana k adalah konstanta pembanding. dengan menyamakan dM pada persamaan (1) dan (2) diperoleh dt ˆ x1 ut (x, t) dx = kux (x1 , t) − kux (x0 , t) x0



atau



ˆ



x1



ut (x, t) dx = k x0 M. Jamhuri



ˆ



x1



uxx (x, t) dx x0



Persamaan Difusi



(3)



Jika integral kedua ruas dari (3) dihilangkan diperoleh ut = kuxx yang biasa disebut sebagai persamaan difusi atau persamaan panas.



M. Jamhuri



Persamaan Difusi



(4)



Solusi Analitik Sebelum menentukan solusi persamaan difusi (4) pada daerah −∞ < x < ∞ dan t > 0, kita tinjau lebih dahulu solusi persamaan difusi dalam bentuk khusus Q (x, t) = g (p) dengan



x . 4kt Permasalahan disini adalah bagaimana bentuk dari g , untuk itu akan kita lakukan langkah-langkah sebagai berikut: substitusikan Q pada (4), dengan p= √



∂Q ∂t



= =



∂Q ∂x



= =



∂2Q ∂x 2



= =



M. Jamhuri



dg ∂p dp ∂t 1 − pg ′ (p) 2t dg ∂p dp ∂x 1 √ g ′ (p) 4kt



(6)



1



∂ ′ g (p) ∂x 4kt 1 ′′ g (p) 4kt √



(5)



Persamaan Difusi



(7)



sehingga diperoleh Qt 1 − pg ′ (p) 2t



=



pg ′ (p)



=



=



kQxx   1 ′′ k g (p) 4kt 1 ′′ − g (p) 2



g ′′ (p) + 2pg ′ (p) = 0 Solusi dari (8) dapat diperoleh sebagai berikut d2 d g (p) + 2p g (p) dp 2 dp   dg d + 2p dp dp misalkan



dan



Solusi dari ODE (10) adalah



=



0



=



0



dg =v dp   d + 2p v = 0 dp dv dp



=



−2pv



v



=



C1 e −p



M. Jamhuri



(8)



2



Persamaan Difusi



(9)



(10)



selanjutnya substitusikan v pada (9), sehingga diperoleh



ˆ



dg dp



=



dg



=



g



=



2



C1 e −p ˆ 2 C1 e −p dp  ˆ 2 e −p dp + C2 C1



dan Q (x, t) = C1



ˆ



√x 4kt



2



e −p dp + C2



0



Konstanta C1 dan C2 diperoleh dengan menggunakan syarat awal khusus, yang diberikan dalam bentuk ( 1, untuk x > 0 Q (x, 0) = 0, untuk x < 0



M. Jamhuri



Persamaan Difusi



Hitung limit t → 0+ Kasus x > 0



lim Q (x, t) = C1



t→0+



ˆ



∞



e



−p 2



dp + C2 = C1



0



√ π + C2 = 1 2



Dalam menghitung integral tak wajar, kita gunakan distribusi normal berbentuk ˆ ∞ 2 1 e −p dp = 1 √ π −∞ Kasus x < 0 lim Q (x, t) = C1



t→0−



ˆ



∞



0



2



e −p dp+C2 = −C1



ˆ



0



−∞



2



e −p dp+C2 = −C1



√



π +C2 = 0 2



Dari dua limit diatas diperoleh 1 C1 = √ π



dan



C2 =



sehingga Q (x, t) =



1 1 +√ 2 π



ˆ



√x 4kt



1 2 2



e −p dp



0



untuk t > 0. Dari Q yang sudah diperkenalkan di atas, kita akan menentukan solusi u terkait dengan Q. Tetapi lebih dahulu kita perhatikan sifat-sifat berikut.



M. Jamhuri



Persamaan Difusi



Jika u memenuhi ut − kuxx = 0 maka v = ux juga memenuhi persamaan tersebut. Kita dapat menunjukkan dengan memeriksa apakah v memenuhi persamaan, turunan dari v   ∂u ∂ vt = ∂t ∂x = vx



= =



vxx



=



∂2u ∂t∂x   ∂ ∂u ∂x ∂x



∂2u ∂x 2  2  ∂ u ∂ ∂x ∂x 2



∂3u ∂x 3 diatas pada persamaan difusi, yaitu =



Selanjutnya terapkan vt , dan vxx vt − kvxx



= = = =



∂2u ∂3u −k 3 ∂t∂x ∂x   ∂ ∂u ∂2u −k 2 ∂x ∂t ∂x ∂ ·0 ∂x 0 memenuhi persamaan difusi. M. Jamhuri



Persamaan Difusi



Dengan Q seperti didefinisikan diatas, S (x, t) =



∂Q ∂x



juga solusi persamaan panas. Hal ini dapat ditunjukkan, karena Q memenuhi persamaan panas, dan sifat sebelum ini, Begitu juga S (x, y ) memenuhi persamaan panas, dan juga ˆ ∞ S (x − y , t) g (y ) dy W (x, t) = −∞



untuk sebarang g (y ) asalkan integral konvergen. Dengan sifat-sifat diatas dan pendefinisian S terkait dengan Q, maka u dapat didefinisikan sebagai ˆ ∞ u (x, t) = S (x − y , t) φ (y ) dy −∞



untuk t > 0, yang memenuhi persamaan panas. Masalah sekarang adalah apakah u tersebut memenuhi kondisi awal u (x, 0) = φ (x) . Untuk itu, kita tuliskan u dalam dalam Q ˆ ∞ ∂Q u (x, t) = (x − y , t) φ (y ) dy −∞ ∂x sedangkan ∂Q ∂Q ∂ (x − y ) ∂y ∂ (x − y ) ∂Q = =− ∂x ∂y ∂ (x − y ) ∂x ∂y ∂x M. Jamhuri



Persamaan Difusi



Selanjutnya gunakan integral parsial, sehingga diperoleh   ˆ ∞ Q (x − y , t) φ′ (y ) dy u (x, t) = − Qφ|∞ −∞ − −∞



Suku pertama pada ruas kanan bernilai nol dengan menggunakan asumsi φ → 0 untuk |y | → ∞, sehingga diperoleh ˆ ∞ Q (x − y , 0) φ′ (y ) dy u (x, 0) = −∞



Sekarang kita gunakan Q (x, 0) = 1



untuk



x > 0 ⇔ Q (x − y , 0) = 1



untuk



y x.



Bila hal ini diterapkan pada integral, didapat ˆ x φ′ (y ) dy = φ (x) u (x, 0) = −∞



memenuhi syarat yang ada, dan secara eksplisit solusinya u (x, t) = √



1 4πkt



M. Jamhuri



ˆ



∞



e



−



(x−y )2 4kt



φ (y ) dy



−∞



Persamaan Difusi



(11)



Contoh Tentukan solusi ut − kuxx = 0 untuk −∞ < x < ∞, dengan syarat awal u (x, 0) = e −x



Dari persamaan 11 diperoleh u (x, t)



= =



(x − y )2 + 4kty 4kt



= = = =



√ √



1



ˆ



4πkt 1



ˆ



4πkt



∞



e



−



(x−y )2 4kt



−







e −y dy



−∞ ∞



e



(x−y )2 +4kty 4kt







dy



−∞



i 1 h (x − y )2 + 4kty 4kt
 1  2 x − xy + y 2 + 4kty 4kt i 1 h (x − y − 2kt)2 + 4ktx − 4k 2 t 2 4kt   x − y − 2kt 2 + (x − kt) √ 4kt



sehingga (12) menjadi e −(x −kt) u (x, t) = √ 4πkt M. Jamhuri



ˆ



∞



2



e −s ds = e −(x −kt)



−∞ Persamaan Difusi



(12)



Metode Pemisahan Variabel Diberikan persamaan difusi ut = 3uxx



pada



0 < x < π,



t>0



(13)



dengan kondisi batas u (0, t)



=



u (π, t) = 0



(14)



u (x, 0)



=



4 sin (2x)



(15)



Misalkan u (x, t) = X (x) T (t) dan substitusikan pemisalan tersebut pada (13), sehingga diperoleh XT ′ = 3X ′′ T X ′′ T′ = (16) 3T X Ruas kiri dari (16) hanya bergantung pada variabel t saja, sedangkan ruas kanan hanya bergantung pada variabel x saja, kondisi tersebut hanya mungkin dipenuhi jika keduanya merupakan konstan yaitu X ′′ T′ = = −λ 3T X



(17)



Misalkan λ = β 2 , maka persamaan (17) dapat dituliskan menjadi dua buah ODE yaitu X ′′ + β 2 X = 0 (18) dan



T ′ + 3λT = 0 M. Jamhuri



Persamaan Difusi



(19)



Solusi dari (18) adalah X (x) = C1 e i βx + C2 e −i βx atau dalam bentuk sinusoidal X (x) = A cos (βx) + B sin (βx)



(20)



Kondisi u (0, t) = 0 memberikan A = 0, sehingga X (x) = B sin (βx) selanjutnya kondisi u (π, t) = 0 memberikan sin (βπ)



=



0



βπ



=



arcsin 0



βπ



=



nπ,



β



=



n



{n = 0, 1, 2, . . . }



sehingga diperoleh Xn (x) = sin (nx)



(21)



Solusi dari persamaan (19) adalah T (t)



=



Ce −3λt



karena λ = β 2 = n2 , maka Tn (t) = Ce −3n M. Jamhuri



2



t



Persamaan Difusi



(22)



Dari persamaan (21) dan (22), maka diperoleh solusi un (x, t) = Cn e −3n



2



t



sin (nx)



Karena kombinasi linier dari solusi persamaan difusi adalah solusi, maka u (x, t) =



∞ X



Cn e −3n



2



t



sin (nx)



n=1



Selanjutnya gunakan kondisi awal (15) u (x, 0) = 4 sin (2x) sehingga diperoleh 4 sin (2x) =



∞ X



Cn sin (nx)



n=1



dimana Cn



= =



8 π (



ˆ



π



sin (2x) sin (nx) dx 0



0, 4



jika n 6= 2 n lainnya



Substitusikan kembali Cn pada (23) sehingga diperoleh u (x, t) = 4e −12t sin (2x) M. Jamhuri



Persamaan Difusi



(23)



Metode Numerik dengan RBF Persamaan difusi (13) yaitu ut = 3uxx kita aproksimasi dengan jaringan RBF sebagai N X



αj



N X ∂2 ∂ φ (x, t) = 3 αj 2 φ (x, t) ∂t ∂x j=1



N X



αj







j=1



j=1



 ∂ ∂2 φ (x, t) − 3 2 φ (x, t) = 0 ∂t ∂x



dimana φ (x, t) =



(24)



q (x − c)2 + (t − d)2 + ǫ2



∂ φ (x, t) ∂t



=



∂2 φ (x, t) ∂x 2



=



q



t −d 2



(x − c) + (t − d)2 + ǫ2



(t − d)2 + ǫ2 h i3 2 (x − c)2 + (t − d)2 + ǫ2



{α}N j=1 adalah koefisien interpolan atau bobot jaringan yang akan ditentukan, sedangkan c dan d adalah center dari jaringan, dan ǫ adalah parameter bebas yang harus dipilih. M. Jamhuri



Persamaan Difusi



Berikutnya aproksimasi kondisi batas (14) memberikan N X



αj φ (0, t) = 0



(25)



N X



αj φ (π, t) = 0



(26)



αj φ (x, 0) = 4 sin (2x)



(27)



j=1



dan



j=1



Dari kondisi batas (15) diperoleh N X j=1



Untuk mendapatkan solusi numerik dari persamaan difusi (13) dengan kondisi batas (14) dan (15), pertama kita harus menentukan koefisien α dari sistem persamaan (24), (25), (26), dan (27). Selanjutnya gunakan α yang didapat untuk menentukan solusi u dengan cara mengaproksimasi u sebagai u (x, t) ≈



M. Jamhuri



N X



αj φ (x, t) .



j=1



Persamaan Difusi



Hasil Simulasi



Hasil simulasi metode RBF diatas diperoleh dengan menggunakan 16 buah titik untuk 0 < x < π dan 21 buah titik untuk 0 < t < 1. Parameter ǫ dipilih sebagai ǫ=



var (x) + var (y ) 2



M. Jamhuri



Persamaan Difusi



Plot error mutlak antara metode RBF vs hasil eksak



M. Jamhuri



Persamaan Difusi



Metode Beda Hingga: FTCS Pada tulisan ini akan dibahas beberapa metode beda hingga untuk persamaan difusi ut = kuxx (28) dengan k suatu konsatnta. Metode FTCS (Forward Time Central Space) biasa disebut sebagai metode eksplisit untuk 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 ∆x 2 , sehingga diperoleh persamaan beda sebagai berikut: ujn+1 − ujn ∆t



=k



n n − 2ujn + uj−1 uj+1



∆x 2



(29)



Persamaan (29) dapat disederhanakan sebagai ujn+1 = atau



dengan S =



k∆t . ∆x 2



 k∆t  n n uj+1 − 2ujn + uj−1 + ujn ∆x 2



  n n + uj−1 ujn+1 = (1 − 2S) ujn + S uj+1



M. Jamhuri



Persamaan Difusi



(30)



Stencil untuk metode FTCS pada persamaan difusi dapat dilihat pada gambar berikut:



Kestabilan: Substitusikan ujn = ρn e iaj pada persamaan (30), sehingga diperoleh   ρn+1 e iaj = (1 − 2S) ρn e iaj + S ρn e ia(j+1) + ρn e ia(j−1) ρe iaj ,



Bagi kedua ruas dari persamaan (31) dengan   ρ = (1 − 2S) + S e ia + e −ia =



= =



sehingga diperoleh



(1 − 2S) + S ([cos a + i sin a] + [cos a − i sin a])



(1 − 2S) + 2S cos a 1 + 2S (cos a − 1)



Agar skema stabil, maka |ρ| ≤ 1, yaitu |ρ| −1 −2 −1 0



= ≤ ≤ ≤ ≤



|1 + 2S (cos a − 1)| 1 + 2S (cos a − 1) 2S (cos a − 1) S (cos a − 1) (1 − cos a) S



M. Jamhuri



Persamaan Difusi



≤ ≤ ≤ ≤ ≤



1 1 0 0 1



(31)



min (1 − cos a) = 0, dan max (1 − cos a) = 2, sehingga 2S



≤



S



≤



1 1 2



Jadi skema akan stabil jika S=k



1 ∆t ≤ ∆x 2 2



Konsistensi: Diberikan dua hampiran berikut: ujn+1



=



n uj±1



=



1 1 1 2 ∆t utt |nj + ∆t 3 uttt |nj + utttt |nj + · · · (32) 2 3! 4! 1 1 1 ujn ± ∆x ux |nj + ∆x 2 uxx |nj ± ∆x 3 uxxx |nj + uxxxx |nj + · (33) ·· 2 3! 4!



ujn + ∆t ut |nj +



n n uj+1 + uj−1 = 2ujn + ∆x 2 uxx |nj +



1 uxxxx |nj + · · · 12



(34)



Substitusikan (32) dan (34) pada persamaan (28), sehingga diperoleh ujn + ∆t ut |nj +



1 2 ∆t utt |nj + · · · 2



M. Jamhuri



=



(1 − 2S) ujn +   1 uxxxx |nj + · · · S 2ujn + ∆x 2 uxx |nj + 12 Persamaan Difusi



Contoh Penerapan Metode FTCS Diberikan persamaan difusi ut = 3uxx



pada



0 < x < π,



t>0



(35)



dengan kondisi batas u (0, t)



=



u (π, t) = 0



(36)



u (x, 0)



=



4 sin (2x)



(37)



Persamaan difusi (35) dengan kondisi batas (36), dan (37) diatas akan kita selesaikan secara numerik menggunakan skema FTCS dengan langkah-langkah sebagai berikut. Persamaan (35) kita diskritkan dengan menggunakan persamaan beda (30), yaitu   3∆t n n , S= + uj−1 ujn+1 = (1 − 2S) ujn + S uj+1 (38) ∆x 2 sedangkan kondisi batas (36) dan (37) sebagai u1n = 0 uj1



dan =



n =0 uM x



4 sin 2xj








dimana {n = 1, . . . Nt , j = 1, . . . , Mx } dengan Nt =



j



T −0 ∆t



k



dan Mx =



Contoh, misalkan untuk j = 2 dan n = 1, maka (38) menjadi
u22 = (1 − 2S) u21 + S u31 + u11 M. Jamhuri



Persamaan Difusi



 π−0  ∆x



.



Simulasi metode beda hingga FTCS



M. Jamhuri



Persamaan Difusi



Error mutlak: metode beda hingga vs hasil eksak



M. Jamhuri



Persamaan Difusi



Metode Implisit BTCS
Metode BTCS memiliki akurasi O ∆t, ∆x 2 , persamaan beda untuk persamaan difusi dengan menggunakan metode BTCS adalah ujn+1 − ujn ∆t



=k



n+1 n+1 − 2ujn+1 + uj−1 uj+1



∆x 2



 k∆t  n+1 n+1 uj+1 − 2ujn+1 + uj−1 2 ∆x n+1 n+1 = ujn + (2S + 1) ujn+1 − Suj+1 −Suj−1



(39)



ujn+1 − ujn =



dengan S =



k∆t . ∆x 2



Kestabilan: Substitusikan ujn = ρn e iaj ke dalam (40) sehingga diperoleh −Sρe −ia + (2S + 1) ρ − Sρe ia   −S e −ia + e ia + (2S + 1)



1 ρ 1 (1 − cos a) 2S + 1 − ρ (1 − cos a) 2Sρ + ρ



−2S cos a + 2S + 1 −



ρ M. Jamhuri



= =



1 1 ρ



=



0



=



0



=



1



=



1 (1 − cos a) 2S + 1



Persamaan Difusi



(40)



Karena untuk setiap S dan a penyebut selalu lebih besar atau sama dengan 1, maka jelas bahwa |ρ| ≤ 1 jadi skema stabil untuk setiap S =



k∆t . ∆x 2



Perhtikan persamaan beda (40) diatas, jika diberikan syarat batas bertipe dirichlet yaitu u (0, t) = f1 dan u (L, t) = f2 . Titik-titik yang harus dihitung adalah ujn+1



M. Jamhuri



Persamaan Difusi



Contoh penerapan metode BTCS Diberikan persamaan difusi ut = 3uxx



pada



0 < x < π,



t>0



(41)



dengan kondisi batas u (0, t)



=



u (π, t) = 0



(42)



u (x, 0)



=



4 sin (2x)



(43)



Persamaan beda skema BTCS untuk persamaan (41) adalah ujn − ujn−1



=



ujn−1



=



h i n n ujn − S uj+1 − 2ujn + uj−1



=



∆t



ujn



−



n n + (1 + 2S) ujn − Suj+1 −Suj−1



atau



=



3



n n − 2ujn + uj−1 uj+1



∆x 2 i 3∆t h n n uj+1 − 2ujn + uj−1 2 ∆x 3∆t ujn−1 , S= ∆x 2 ujn−1



n n = −ujn−1 − (1 + 2S) ujn + Suj+1 Suj−1



M. Jamhuri



Persamaan Difusi



(44)



Kondisi batas (42) kita diskritkan sebagai u1n = 0,



dan



n =0 uM x



(45)








(46)



dan (43) kita diskritkan sebagai uj1 = 4 sin 2xj



dimana {j = 1, . . . , Mx , n = 1, . . . , Nt } dengan Mx =



 π−0  ∆x



dan Nt =



j



T −0 ∆t



k



.



Dalam bentuk matrik dapat kita gambarkan persamaan beda (44), (45), dan (46) sebagai j \n 1 2 3 . .. Mx − 1 Mx



1 0 4 sin (2x2 ) 4 sin (2x3 ) . .. 4 sin (2xMx −1 ) 0



2 0 u22 u32 . ..



3 0 u23 u33 . ..



2 uM x −1 0



3 uM x −1 0



M. Jamhuri



··· ··· ··· ··· .. . ··· 0



Persamaan Difusi



Nt − 1 0 u2Nt −1 u3Nt −1 . .. Nt −1 uMx −1 0



Nt 0 u2Nt u3Nt . .. Nt uM x −1 0



Sebagai contoh, untuk j = 2, 3, . . . , Mx − 1 dan n = 2 akan kita tentukan ujn , yaitu u22 =? maka dengan menggunakan (44) diperoleh n n = −ujn−1 − (1 + 2S) ujn + Suj+1 Suj−1



j =2 j =3 j =4 .. . j = Mx − 1



⇒ ⇒ ⇒ .. . ⇒



Su12 − (1 + 2S) u22 + Su32 Su22 − (1 + 2S) u32 + Su42 Su32 − (1 + 2S) u42 + Su52 .. . 2 2 2 + SuM − (1 − 2S) uM SuM x x −1 x −2



= = = .. . =



−u21 −u31 −u41 .. . 1 −uM x −1



so we have matrix         



− (1 + 2S) S 0



S − (1 + 2S) S



. . . 0



. . . 0



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



··· ··· ··· .



.



. ···



M. Jamhuri



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



         



2 u2 2 u3 2 u4 . . .



2 uM x −1



Persamaan Difusi











         =         



1 − Su 2 −u2 1 1 −u3 1 −u4 . . . 1 2 −uM − SuM x −1 x



         