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
dan dengan uraian yang sama diperoleh Q (x − y , 0) = 0 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
