KNM XVI
3-6 Juli 2012
UNPAD, Jatinangor
SKEMA BEDA HINGGA TAK-STANDAR UNTUK MODEL EPIDEMI DENGAN LAJU PENULARAN TERSATURASI YANG DIMODIFIKASI AGUS SURYANTO Jurusan Matematika Universitas Brawijaya Jl. Veteran Malang 65145; Email:
[email protected]
Abstrak Skema beda hingga tak standar untuk model epidemi SIR dengan laju penularan tersaturasi yang dimodifikasi akan dikonstruksi. Konstruksi skema dilakukan dengan mengaproksimasi turunan pertama dengan generalisasi beda maju dan mendekati sukusuku lainnya dengan pendekatan tak-lokal. Fungsi penyebut untuk generalisasi beda maju dikonstruksi berdasarkan hukum konservasi populasi model epidemi. Berbeda dengan metode Euler dan metode Runge-Kutta, skema yang dikembangkan tetap menjaga sifatsifat dinamik dari persamaan diferensial SIR seperti hukum konservasi populasi, kepositifan penyelesaian, titik tetap dan kestabilannya tanpa bergantung pada langkah waktu integrasi. Hasil-hasil simulasi numerik juga menunjukkan bahwa skema yang dikembangkan lebih akurat dibandingkan dengan skema beda hingga tak-standar tanpa mempertimbangkan hukum konservasi. Kata Kunci: model epidemi, skema beda hingga tak standar, konsistensi dinamik, hukum konservasi, fungsi penyebut.
1. Pendahuluan Dalam beberapa dekade terakhir telah banyak dikembangkan model matematika yang menjelaskan dinamika penyakit menular.Pada model-model matematika epidemi umumnya populasi dibagi menjadi beberapa kelas yang berbeda. Salah satu model matematika epidemi yang termasuk dalam kategori ini adalah model epidemi SIR yang membagi populasi menjadi kelompok individu rentan terhadap infeksi suatu penyakit yang dinotasikan dengan S (Susceptible), kelompok individu terinfeksi yang dinotasikan dengan I (Infective) dan kelompok individu sembuh dari infeksi yang dinotasikan dengan R (Removed). Salah satu mekanisme penting dalam model matematika epidemi SIR adalah laju penularan yaitu laju perubahan individu rentan (S) menjadi individu terinfeksi (I).Pada Model klasik yang dikembangkan oleh Kermack and McKendrik [1], laju penularan dimodelkan dalam bentuk fungsi bi-linear. Dalam hal ini, penularan penyakit berbanding lurus baik dengan jumlah individu rentan maupun individu terinfeksi, sehingga semakin besar jumlah individu rentan maupun individu terinfeksi akan menyebabkan tingkat penularan yang semakin besar pula.Secara umum jika jumlah individu terinfeksi relatif kecil maka laju penularan dapat dimodelkan secara bi-linear. Tetapi jika jumlah individu terinfeksi sangat banyak maka perilaku sosial, psikologi masyarakat maupun mekanisme lain akan mengakibatkan perubahan laju penularan. Berkaitan dengan hal tersebut beberapa peneliti membuat model epidemi dengan laju
ISBN : 978-602-19590-2-2
879
Suryanto A.
Skema Beda Hingga …
SI . Di sini 1 1S efek faktor saturasi menjelaskan kontrol epidemi akibat tindakan preventif dari subpopulasi rentan, lihat [2-3].Selain itu, terdapat juga model yang menggunakan laju SI penularan tersaturasi yang berbentuk . Pada model kedua ini kontak efektif 1 2I antara subpopulasi terinfeksi dengan subpopulasi rentan dapat tersaturasi pada tingkat infeksi yang tinggi akibat padatnya subpopulasi terinfeksi ataupun akibat proteksi subpopulasi rentan, lihat [4-5].Kombinasi kedua model telah dilakukan oleh Pathak [6]. Jika diasumsikan bahwa individu yang sembuh memiliki imunitas sehingga tidak dapat tertular lagi, maka model tersebut berbentuk: dS SI A S dt 1 1S 2 I
penularan tersaturasi.Salah satunya adalah laju penularan yang berbentuk
dI SI I dt 1 1S 2 I
(1)
dR I R. dt Laju penularan pada model SIR (1) menggambarkan efek saturasi dimana kontak efektif antar individu rentan dan terinfeksi dapat mencapai suatu nilai maksimum akibat perubahan perilaku sosial atau psikologi populasi.
Perlu diperhatikan bahwa model matematika penyebaran penyakit (1) berbentuk sistem persamaan diferensial nonlinear dimana penyelesaian umumnya sulit diperoleh.Oleh karena itu, metode pendekatan numerik diperlukan untuk mendapatkan penyelesaian yang akurat. Metode numerik akan mentransformasikan sistem persamaan diferensial ke dalam sistem persamaan beda. Tentu saja sistem persamaan beda tersebut diharapkan konsisten secara dinamik dengan sistem persamaan diferensialnya, yaitu jika kedua sistem tersebut mempunyai sifat-sifat dinamik yang sama, misalnya titik tetap dan kestabilannya, kepositifanpenyelesaian, konservasi jumlah populasi dan lain sebagainya [7-9]. Dalam artikel ini akan dikonstruksi skema numerik yang konsisten secara dinamik untuk model epidemi SIR (1). Skema dikonstruksi menggunakan beda hingga tak-standar. Hasil-hasil numerik menggunakan skema yang dikonstruksi akan dibandingkan dengan hasil-hasil metode Euler dan metode Runge-Kutta orde empat.
2. Sifat-Sifat Dinamik Model SIR dengan Laju Penularan Tersaturasi yang Dimodifikasi Pada bagian ini akan dibahas secara ringkas tentang sifat-sifat dinamik model epidemi SIR dengan laju penularan terinfeksi yang dimodifikasi (1). Sifat-sifat tersebut akan digunakan sebagai acuan dalam membahas sifat-sifat dinamik model diskrit hasil diskritisasi model (1) menggunakan beda hingga tak standar. Pertama, perlu diperhatikan bahwa total populasi merupakan jumlah seluruh individu S, I dan R, yaitu N = S + I + R. Persamaan (1) menjelaskan dinamika populasi sehingga perilaku dinamik akan dipelajari pada domain S 0, I 0, R 0 | S I R N .
KNM XVI - 3-6 Juli 2012 – UNPAD, Jatinangor
880
KNM XVI
3-6 Juli 2012
UNPAD, Jatinangor
Dengan demikian semua subpopulasi S, I dan R tidak boleh bernilai negatif.Kondisi ini disebut kondisi kepositifan populasi.Jika persamaan-persamaan pada sistem (1) dijumlahkan maka diperoleh persamaan diferensial untuk total populasi dN A N dt . (2) Persamaan (2) disebut hukum konservasi untuk populasi yang berkaitan dengan sistem (1).Solusi persamaan (2) adalah 1 N t A A N 0 e t (3) dimana N 0 N 0 S0 I 0 R0 dengan masing-masing S 0 S0 , I 0 I 0 dan
R0 R0 menyatakan jumlah subpopulasi awal. Dari persamaan (3) dapat dilihat bahwa total populasi tidak bernilai konstan. Lebih tepatnya total populasi merupakan fungsi monoton naik jika N0 A / dan monoton turun jika N0 A / . Pada kedua kasus tersebut berlaku lim N t A / , sehingga bidang S I R A / adalah manifold t
invarian yang stabil pada oktan pertama. Selanjutnya, Pathak [6] mendefinisikan angka reproduksi dasar untuk model SIR (1), yaitu A 1 A 0 . (4) 1 Selain itu, Pathak [6] juga menunjukkan bahwa jika 0 maka model SIR (1) mempunyai titik kesetimbangan tunggal yang bersifat stabil asimtotik, yaitu
S0*, I0*, R0* A / , 0, 0 .Titik
kesetimbangan tersebutdikenal sebagai titik kesetimbangan bebas penyakit. Apabila 0 1 maka selain titik kesetimbangan bebas
* * * penyakit, model SIR (1) juga mempunyai titik kesetimbangan endemi Se , I e , Re dengan
A 0 1 A 2 A I * ; I e* Re* e 0 2 A 0 2 A dan .Pada kasus ini, titik bebas penyakit tidak stabil sedangkan titik endemi stabil asimtotik. S e*
3. Skema Beda Hingga Tak Standar Untuk mengkonstruksi skema numerik, variabel waktu t didiskritisasi pada titik-titik tn nh n 0,1,2,... dengan h 0 adalah lebar grid atau sering disebut lebar langkah waktu integrasi. Solusi pada titik tn masing-masing adalah S tn , I tn dan Rtn . Solusi numerik pada titik yang sama t n yang berkaitan dengan solusi tersebut dinyatakan sebagai Sn , I n dan Rn . Selanjutnya skema numerik dikonstruksi dengan menerapkan metode beda hingga tak standar yang dikembangkan oleh Mickens [7-9]. Dengan menerapkan generalisasi beda maju dan pendekatan tak lokal terhadap persamaan (1), diperoleh model numerik tak standar untuk model epidemi (1) yaitu
ISBN : 978-602-19590-2-2
881
Suryanto A.
Skema Beda Hingga …
S n 1 S n I n S n 1 A S n 1 , h 1 1 S n 2 I n I n 1 I n I n S n 1 I n 1 , h 1 1 S n 2 I n
(5)
Rn 1 Rn I n 1 Rn 1. h Perhatikan bahwa turunan didekati dengan beda maju tetapi dengan menggunakan fungsi penyebut h , sehingga dikenal sebagai perluasan beda maju.Jika dipilih fungsi penyebut h h , maka skema (5) tetap merupakan skema beda hingga tak standar karena menggunakan pendekatan tak lokal. Pada artikel ini akan dipilih fungsi penyebut sehingga total populasi eksak tetap terjaga, yaitu dengan cara berikut.Dengan menjumlah semua persamaan pada sistem diskrit (5), dapat diperoleh persamaan N n 1 N n (6) A N n 1 . h Perhatikan bahwa persamaan (6) merupakan bentuk diskritisasi tak-standar dari persamaan konservasi populasi (2) sehingga fungsi penyebut h dapat diperoleh dengan membandingkan solusi persamaan beda (6) dengan solusi persamaan (2), yaitu 1 solusi pada persamaan (3) pada saat t tn 1 , yaitu N n 1 A A N n exp h . Dengan cara ini diperoleh bahwa fungsi penyebut harus berbentuk berbentuk exph 1 . (7) h Dengan menggunakan fungsi penyebut (7), persamaan beda (6) mempunyai solusi yang tepat sama dengan solusi (3) yaitu 1 N n A A N 0 e nh . (8) Persamaan (8) menunjukkan bahwa hukum konservasi populasi selalu terpenuhi untuk sembarang nilai h. Selanjutnya perhatikan bahwa skema (5) merupakan skema implisit karena menggunakan pendekatan tak-lokal. Tetapi skema tersebut dapat disusun kembali dengan mudah untuk mendapatkan skema eksplisit yaitu S n A h S n 1 1 h S n , I n
I n h S n , I n S n 1 1 h R h I n 1 Rn 1 n , 1 h I n 1
(9)
dengan S , I I / 1 1S 2 I . Karena semua parameter dan fungsi penyebut bernilai tk-negatif, skema (9) akan selalu menghasilkan penyelesaian tak-negatif. Analog dengan kasus kontinu, kondisi kepositifan penyelesaian tersebut bersama-sama dengan hukum konservasi populasi menjamin keterbatasan penyelesaian numerik yang dihasilkan oleh skema (5) atau (9).
KNM XVI - 3-6 Juli 2012 – UNPAD, Jatinangor
882
KNM XVI
3-6 Juli 2012
UNPAD, Jatinangor
4. Simulasi Numerik Pada bagian ini ditunjukkan beberapa hasil numerik untuk mengkonfirmasikan hasil analisis sebelumnya.Simulasi dilakukan dengan mengunakan parameter seperti pada Tabel 1. Tabel 1. Nilai-nilai parameter untuk simulasi Parameter A Titik bebas Titik endemi penyakit Nilai 0.13 0.01 0.94
0.05
0.1
0.01
0.05
4.1 Titik Kesetimbangan Bebas Penyakit ( 0 1 )
40 h = 0.1
30
R 20 10 Titik Bebas Penyakit 0 15 10
I
5 0
0
5
10
15
20
S
Gambar 1. Diagram fasa untuk skenario titik bebas penyakit menggunakan skema beda hingga tak standar (11) dengan h = 0.1. Skema beda hingga tak standar (5) atau bentuk eksplisitnya (9) mempunyai titik kesetimbangan bebas penyakit yang stabil asimtotik. Untuk sifat kestabilan, pada Gambar 1 ditunjukkan hasil simulasi numerik menggunakan parameter seperti pada Tabel 1 dengan 1 0.13 , h = 0.1 dan beberapa nilai awal. Dengan menggunakan parameter tersebut diperoleh 0 0.9742 1 sehingga titik bebas penyakit (18.8, 0.0, 0.0) merupakan titik kesetimbangan yang stabil.Gambar 1 menunjukkan bahwa solusi semua numerik dengan berbagai nilai awal konvergen ke titik bebas penyakit. Selanjutnya pada Gambar 2 ditunjukkan perbandingan solusi numerik menggunakan skema beda hingga tak standar (9), beda hingga tak standar seperti pada skema (9) tetapi dengan fungsi penyebut h h , metode Euler dan metode Runge-Kutta orde 4. Pada gambar tersebut digunakan langkah waktu integrasi yang cukup besar yaitu h = 5 dan nilai awal (20, 2, 15). Terlihat bahwa solusi numerik yang dihasilkan oleh seluruh metode tersebut mempunyai perilaku dinamik yang sama yaitu semua solusi bernilai tak negatif dan konvergen ke titik bebas penyakit. Tetapi, error total populasi dari masing-masing metode menunjukkan perbedaan yang cukup signifikan, lihat Gambar 3. Pada Gambar 3
ISBN : 978-602-19590-2-2
883
Suryanto A.
Skema Beda Hingga …
tersebut ditunjukkan error mutlak dari total populasi dalam skala logaritma. Terlihat bahwa error yang dihasilkan oleh skema beda hingga tak standar (9) sangat kecil yaitu dalam orde 10 14 . Metode beda hingga tak standar dengan h h dan metode Euler mempunyai error yang hampir sama dan sangat besar karena keduanya mempunyai akurasi orde pertama. Metode Runge-Kutta orde 4 mempunyai error yang relatif kecil dibandingkan dengan skema beda hingga tak standar dengan h h maupun skema Euler tetapi relatif sangat besar dibandingkan dengan metode beda hingga tak standar (9). (a) Skema Tak Standar; = (exp(h) - 1) / 20
20
15
15 S I R
10
(b) Skema Tak Standar; = h
S I R
10 5
5 0 0
50
100 t
150
0 0
200
(c) Skema Euler
20
50
15
S I R
10
150
S I R
10
5
200
(d) Skema Runge Kutta order 4
20
15
100 t
5
0 0
50
100 t
150
0 0
200
50
100 t
150
200
Gambar 2. Perbandingan solusi numerik untuk skenario titik bebas penyakit dengan h = 5. Solusi numerik dari masing-masing metode konvergen ke titik bebas penyakit. 0
10
-5
10 error mutlak
Tak Standar; = (exp(h)-1)/ Tak Standar; = h Euler -10
10
Runge-Kutta orde 4 error bernilai nol
-15
10
0
25
50
75
100 t
125
150
175
200
Gambar 3. Perbandingan error mutlak total populasi untuk simulasi seperti pada Gambar 4.
4.2 Titik Kesetimbangan Endemi ( 0 1 ) Selanjutnya dilakukan simulasi dengan metode beda hingga tak standar (9) untuk skenario titik endemi. Di sini juga parameter seperti pada Tabel 1 dengan 1 0.13 . Pada Gambar 4 ditunjukkan hasil simulasi menggunakan h = 0.1 dengan beberapa nilai
KNM XVI - 3-6 Juli 2012 – UNPAD, Jatinangor
884
KNM XVI
3-6 Juli 2012
UNPAD, Jatinangor
awal yang berbeda. Pada kasus ini terdapat dua titik kesetimbangan yaitu titik kesetimbangan bebas penyakit(18.8, 0.0, 0.0) dan titik endemi (5.89, 1.17, 11.74).Meskipun demikian, dapat dilihat bahwa semua solusi konvergen ke titik endemi. Hal ini sesuai dengan analisis bahwa karena 0 1 , titik endemi merupakan titik kesetimbangan yang stabil.
20 15
R 10 Titik Endemi 5 Titik Bebas Penyakit
0 15 10 5
I
0
10
5
0
20
15
S
Gambar 4. Diagram fasa untuk skenario titik endemi menggunakan skema beda hingga tak standar (11) dengan h = 0.1.
(a) Skema Tak Standar; = (exp(h) - 1) / 25
25
20
20
15
15
10
10
5
5
0 0
50
100 t
150
200
(c) Skema Euler
400
0 0
400
200
200
0
0
-200
-200
-400 0
5
10
t
15
20
25
-400 0
(b) Skema Tak Standar; = h
50
100 t
150
200
(d) Skema Runge Kutta order 4
2
4 t
6
8
Gambar 5. Perbandingan solusi numerik untuk skenario titik endemi dengan h = 5. Solusi numerik dari metode beda hingga baik skema (11) maupun skema dengan h h konvergen ke titik endemi; sedangkan metode Euler dan Runge-Kutta menghasilkan solusi bernilai negatif dan tak terbatas (plot hanya pada interval [-400, 400].
Pada Gambar 5 dibandingkan solusi numerik dari skema beda hingga tak standar (9), skema beda hingga tak standar dengan h h , skema Euler dan skema RungeKutta orde 4 menggunakan h = 5. Gambar 5 menunjukkan bahwa kedua beda hingga tak
ISBN : 978-602-19590-2-2
885
Suryanto A.
Skema Beda Hingga …
standar mempunyai perilaku dinamik yang sama yaitu kedua solusi konvergen ke titik endemi. Tetapi metode Euler dan Runge-Kutta orde 4 tidak stabil.Dalam hal ini, solusi metode Euler maupun Runge-Kutta dapat bernilai negatif dan tak terbatas. Meskipun kedua metode beda hingga tak standar mempunyai perilaku dinamik yang sama tetapi dapat dilihat pada Gambar 6 bahwa total populasi dari metode beda hingga dengan h h mempunyai error mutlak yang relatif jauh lebih besar dibandingkan dengan skema (9). Hal ini menunjukkan bahwa skema (6) menjaga hokum konservasi populasi.
5. Kesimpulan Dalam tulisan ini telah dikembangkan metode numerik untuk menyelesaikan model epidemi SIR dengan laju penularan tersaturasi yang dimodifikasi menggunakan metode beda hingga tak-standar. Skema numerik diturunkan dengan memperkenalkan fungsi penyebut dalam mendekati turunan pertama dan menerapkan pendekatan tak-lokal untuk suku-suku lainnya.Fungsi penyebut diturunkan dari hukum konservasi populasi sehingga total populasi eksak selalu terpenuhi. Hasil-hasil simulasi numerik menunjukkan bahwa skema yang dikembangkan konsisten secara dinamik dengan model epidemi SIR kontinu untuk berapapun langkah waktu integrasi yang digunakan. 5
10
0
10
error mutlak
Tak Standar; = (exp(h)-1)/ Tak Standar; = h -5
10
Euler Runge-Kutta orde 4 error bernilai nol
-10
10
-15
10
0
20
40
60
80
100 t
120
140
160
180
200
Gambar 6. Perbandingan error mutlak total populasi untuk simulasi seperti pada Gambar 5.
Ucapan Terima Kasih Penelitian ini dibiayai oleh Direktorat Jenderal Pendidikan Tinggi, Kementerian Pendidikan dan Kebudayaan melalui DIPA Universitas Brawijaya No. : 0636/02304.2.16/15/2012, dan berdasarkan SK Rektor Universitas Brawijaya No. : 058/SK/2012.
KNM XVI - 3-6 Juli 2012 – UNPAD, Jatinangor
886
KNM XVI
3-6 Juli 2012
UNPAD, Jatinangor
Daftar Pustaka [1] Kermack, W.O. dan A.G. McKendrick, A Contribution to the Mathematical Theory of Epidemic, Proc. R. Soc. London A 115: 700-721, 1927. [2] C. Wei dan L. Chen, A delayed epidemic model with pulse vaccination, Discrete Dynamics in Nature and Society, 2008, Article ID 746951. [3] Zhang, J.-Z., Jin, Z., Liu, Q.-X. dan Zhang, Z.-Y., Analysis of a delayed SIR model with nonlinear incidence rate, Discrete Dynamics in Nature and Society2008, Article ID 66153. [4] Xu, R. dan Ma, Z., Stability of a delayed SIRS epidemic model with a nonlinear incidence rate, Chaos, Solitons and Fractals, 41(5): 2319-2325, 2009. [5] A. Suryanto, A Dynamically Consistent Nonstandard Numerical Scheme for Epidemic Model with Saturated Incidence Rate, Int. J. of Mathematics and Computation 13 (D11), 112-123, 2011. [6] Pathak, S., Maiti, A. dan Samanta, G.P., Rich Dynamics of an SIR Epidemic Model, Nonlinear Analysis: Modelling and Control15(1), 71–81, 2010. [7] Mickens, R.E., Nonstandard Finite Difference Schemes, World Scientific, Singapore 1994. [8] Mickens, Application of Nonstandard Finite Difference Schemes, World Scientific, Publishing Co Pte. Ltd., 2000. [9] Mickens,R.E., Dynamic consistency: a fundamental principle for constructing nonstandard finite difference schemes for differential equations, Journal of Difference Equations and Applications11 (7): 645-653, 2005.
ISBN : 978-602-19590-2-2
887
Suryanto A.
KNM XVI - 3-6 Juli 2012 – UNPAD, Jatinangor
Skema Beda Hingga …
888
