TUGAS AKHIR - SM141501
EKSISTENSI BIFURKASI MUNDUR DAN KENDALI OPTIMAL PADA MODEL PENYAKIT VEKTOR-BORNE YANG DISEBABKAN NYAMUK CHARISMA JUNI KUMALASARI NRP 1211 100 032 Dosen Pembimbing: Subchan, M.Sc, Ph.D Drs. M. Setijo Winarko,M.Si JURUSAN MATEMATIKA Fakultas Matematika dan Ilmu Pengetahuan Alam Institut Teknologi Sepuluh Nopember Surabaya 2015
FINAL PROJECT - SM141501
EXISTENCE OF BACKWARD BIFURCATION AND OPTIMAL CONTROL VECTOR-BORNE DISEASE MODELS ON CAUSED MOSQUITO CHARISMA JUNI KUMALASARI NRP 1211 100 032 Supervisors: Subchan, M.Sc, Ph.D Drs. M. Setijo Winarko,M.Si DEPARTMENT OF MATHEMATICS Faculty of Mathematics and Natural Sciences Sepuluh Nopember Institute of Technology Surabaya 2015
EKSISTENSI BIFURKASI MUNDUR DAN KENDALI OPTIMAL PADA MODEL PENYAKIT VEKTOR-BORNE YANG DISEBABKAN NYAMUK Nama Mahasiswa NRP Jurusan Pembimbing
: : : :
Charisma Juni Kumalasari 1211 100 032 Matematika FMIPA-ITS 1. Subchan, M.Sc, Ph.D 2. Drs. M. Setijo Winarko,M.Si
Abstrak Penyakit vector-borne merupakan penyakit yang terjadi pada manusia yang penularannya melalui vektor (perantara) seperti serangga. Contoh penyakit vector-borne seperti demam berdarah, virus West Nile, virus encephalitis dan malaria [1]. Model yang digunakan merupakan kombinasi dari dua model non linear dari populasi individu dan vektor [1]. Populasi individu dikelompokkan menjadi susceptible, infected, dan recovered, sedangkan vektor dikelompokkan menjadi susceptible dan infected. Dalam Tugas Akhir ini membahas tentang analisa pada model penyakit vector-borne dengan menyelidiki adanya bifurkasi mundur, menentukan basic reproduction number, kestabilan dari setiap titik kesetimbangan, kendali optimal, dan solusi numerik dari model interaksi dinamis yang disimulasikan dengan menggunakan MATLAB, sehingga didapatkan hasil untuk meminimalkan jumlah host (inang) yang terinfeksi serta jumlah populasi vektor. Kata-kunci:
Model epidemik, Bifurkasi Mundur, Kendali Optimal, Prinsip Minimum Pontriagyn (PMP)
vii
Halaman ini sengaja dikosongkan.
EXISTENCE OF BACKWARD BIFURCATION AND OPTIMAL CONTROL VECTOR-BORNE DISEASE MODELS ON CAUSED MOSQUITO Name NRP Department Supervisors
: : : :
Charisma Juni Kumalasari 1211 100 032 Mathematics FMIPA-ITS 1. Subchan, M.Sc, Ph.D 2. Drs. M. Setijo Winarko,M.Si
Abstract Vector-borne diseases are diseases that occur in humans are transmitted through a vector (intermediaries) such as insects. Examples of vector-borne diseases such as dengue fever, West Nile virus, encephalitis virus and malaria [1]. The model is used a combination of two non-linear models of the individual and vector population [1]. The population of individuals is classified into susceptible, infected, and recovered, while the vectors classified into susceptible and infected. In this final project is about the analysis on the model of vector-borne diseases to investigate the existence of backward bifurcation, determine the basic reproduction number, the stability of each equilibrium point, optimal control, and the numerical solution of the dynamic interaction models are simulated using MATLAB, so we get the results to minimize the number of hosts (host) and the number of infected vector population. Key-words:
Epidemic model, Optimal control, Principle (PMP)
ix
Backward bifurcation, Pontryagins Minimum
Halaman ini sengaja dikosongkan.
KATA PENGANTAR
Assalamu’alaikum Wr. Wb. Alhamdulillaahirobbil’aalamiin, segala puji dan syukur penulis panjatkan ke hadirat Allah SWT yang telah memberikan limpahan rahmat, petunjuk serta hidayah-Nya, sehingga penulis dapat menyelesaikan Tugas Akhir yang berjudul ”EKSISTENSI BIFURKASI MUNDUR DAN KENDALI OPTIMAL PADA MODEL PENYAKIT VEKTOR-BORNE YANG DISEBABKAN NYAMUK” sebagai salah satu syarat kelulusan Program Sarjana Jurusan Matematika FMIPA Institut Teknologi Sepuluh Nopember (ITS) Surabaya. Tugas Akhir ini dapat terselesaikan dengan baik berkat bantuan dan dukungan dari berbagai pihak. Oleh karena itu, penulis menyampaikan ucapan terima kasih dan penghargaan kepada: 1. Ibu Dr. Erna Apriliani, M.Si selaku Ketua Jurusan Matematika FMIPA ITS. 2. Bapak Drs. M. Setijo Winarko, M.Si dan Bapak Subchan, M.Sc, Ph.D selaku dosen pembimbing atas segala bimbingan dan motivasinya kepada penulis dalam mengerjakan Tugas Akhir ini sehingga dapat terselesaikan dengan baik. 3. Bapak Dr. Hariyanto, M.Si, Drs. Iis Herisman, M.Si dan Drs. Daryono Budi Utomo, M.Si selaku dosen xi
penguji atas semua saran yang telah diberikan demi perbaikan Tugas Akhir ini. 4. Bapak Drs. Chairul Imron, MI.Komp selaku Ketua Prodi S-1 Jurusan Matematika FMIPA ITS. 5. Bapak Drs. Soetrisno, MI.Komp selaku dosen wali yang telah memberikan arahan akademik selama penulis menempuh pendidikan di Jurusan Matematika FMIPA ITS. 6. Bapak dan Ibu dosen serta seluruh staff Tata Usaha dan Laboratorium Jurusan Matematika FMIPA ITS. 7. Teman teman angkatan 2011 Jurusan Matematika atas dukungan yang telah diberikan kepada penulis. Penulis juga menyadari bahwa dalam Tugas Akhir ini masih terdapat kekurangan. Oleh sebab itu, kritik dan saran yang bersifat membangun sangat penulis harapkan demi kesempurnaan Tugas Akhir ini. Akhirnya, penulis berharap semoga Tugas Akhir ini dapat bermanfaat bagi banyak pihak. Surabaya, Januari 2015
Penulis
xii
Special Thank’s To Keberhasilan penulisan Tugas Akhir ini tidak lepas dari orang-orang terdekat penulis. Oleh sebab itu, penulis mengucapkan terima kasih kepada : 1. Allah SWT yang telah memberi rahmat, petunjuk, kekuatan, dan kesabaran dalam setiap langkah kehidupan penulis serta kepada Nabi Muhammad SAW yang telah membimbing umat-Nya dari zaman jahiliyah menuju zaman yang penuh ilmu. 2. Mbok e, Mama, dan Ayah yang tercinta, terima kasih atas segala doanya, juga kasih sayang dan pendidikan yang selalu dicurahkan kepada penulis selama ini. 3. Kakak-kakak yang sangat ku sayangi, mas Gatot, mas Gepri, dan mas Rizky. Terima kasih banyak atas segala doa, dukungan, motivasi, dan nasihat-nasihatnya kepada penulis. 4. Mas Tony, Mas Ipin, Mas Fahim, Mas Danang, Mbak Fahmi, Mbak Irma, dan Mbak Firdha atas semua bantuan yang diberikan dalam pengerjaan Tugas Akhir ini. 5. Teh Gita, Mb Fia, Veda, Mila dan teman satu atap lainnya yang selalu memberikan semangat serta saransaran kepada penulis. Terima kasih juga anak-anak lab lainnya dan mas Tux yang telah membantu dan menyemangati penulis. Makasih banyak ya semuanya.. 6. Teman-teman seperjuangan Tugas Akhir, Nilam, Ka Marmel, Mas Lulu, Mas Romi, Mb Nadia, Muna, Zamroji, Andika, dan lain-lain yang saling mendukung dan memotivasi satu sama lain. xiii
7. Teman-teman angkatan 2011, terima kasih atas doa dan dukungan kalian selama ini. Kalian merupakan keluarga baru ku disini. 8. Seluruh keluarga besar HIMATIKA ITS dan UKM Taekwondo ITS atas dukungan dan semangat yang diberikan kepada penulis. 9. Semua pihak yang tak bisa penulis sebutkan satupersatu, terima kasih telah membantu sampai terselesaikannya Tugas Akhir ini.
xiv
DAFTAR ISI
HALAMAN JUDUL
i
LEMBAR PENGESAHAN
vi
ABSTRAK
vii
ABSTRACT
ix
KATA PENGANTAR
xi
DAFTAR ISI
xv
DAFTAR GAMBAR
xix
DAFTAR TABEL
xxi
DAFTAR SIMBOL
xxiii
DAFTAR SIMBOL
xxv
BAB I 1.1 1.2 1.3 1.4 1.5 1.6
PENDAHULUAN Latar Belakang . . . . . . . . . . . . . . . . . . . . . . . . . . Rumusan Masalah . . . . . . . . . . . . . . . . . . . . . . . Batasan Masalah . . . . . . . . . . . . . . . . . . . . . . . . Tujuan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Manfaat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Sistematika Penulisan . . . . . . . . . . . . . . . . . . . .
1 1 2 2 3 3 3
BAB II 2.1 2.2 2.3 2.4 2.5
TINJAUAN PUSTAKA Penyakit Vector-Borne . . . . . . . . . . . . . . . . . . . . Sistem Kompartemen . . . . . . . . . . . . . . . . . . . . . Bilangan Reproduksi Dasar (R0 ) . . . . . . . . . . . Kestabilan Titik Tetap . . . . . . . . . . . . . . . . . . . Stabil Asimtotik Lokal . . . . . . . . . . . . . . . . . . . .
5 5 6 6 9 9
xv
2.6 2.7 2.8 2.9
Bifurkasi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Teori Kendali Optimal . . . . . . . . . . . . . . . . . . . . Prinsip Minimum Pontryagin . . . . . . . . . . . . . . Metode Beda Hingga . . . . . . . . . . . . . . . . . . . . .
10 11 12 15
BAB III 3.1 3.2 3.3
METODE PENELITIAN Studi Literatur . . . . . . . . . . . . . . . . . . . . . . . . . . Mengkaji Model Interaksi Dinamis . . . . . . . . . Mencari Titik Kesetimbangan, Menentukan Bilangan Reproduksi Dasar dari Model, dan Bifurkasi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Menentukan Formulasi Masalah Kendali Optimal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Menentukan Penyelesaian Kendali Optimal . . Membuat Simulasi . . . . . . . . . . . . . . . . . . . . . . . Kesimpulan dan Saran . . . . . . . . . . . . . . . . . . . .
19 19 19
3.4 3.5 3.6 3.7 BAB IV 4.1 4.2 4.3 4.4 4.5 4.6
4.7 4.8 4.9 4.10
ANALISIS DAN PEMBAHASAN Deskripsi Model dan Asumsi . . . . . . . . . . . . . . Daerah Penyelesaian Model . . . . . . . . . . . . . . . Titik Kesetimbangan Bebas Penyakit . . . . . . . Menentukan Bilangan Reproduksi Dasar . . . . Titik Kesetimbangan Endemik . . . . . . . . . . . . . Kestabilan Lokal Model Interaksi Dinamis . . 4.6.1 Kestabilan Lokal Titik Setimbang Bebas Penyakit . . . . . . . . . . . . . . . . . . . . 4.6.2 Kestabilan Lokal Titik Setimbang Endemik . . . . . . . . . . . . . . . . . . . . . . . . . . Analisa Bifurkasi . . . . . . . . . . . . . . . . . . . . . . . . . Model Kendali Optimal . . . . . . . . . . . . . . . . . . . Penyelesaian Kendali Optimal . . . . . . . . . . . . . Solusi Numerik . . . . . . . . . . . . . . . . . . . . . . . . . . xvi
19 20 20 20 20 21 21 26 31 33 38 44 53 57 66 71 76 83
BAB V PENUTUP 91 5.1 Kesimpulan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 5.2 Saran . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 DAFTAR PUSTAKA
95
LAMPIRAN A
Kurva Bifurkasi Mundur
97
LAMPIRAN B
Kurva Bifurkasi Maju
99
LAMPIRAN C Simulasi dari populasi individu yang terinfeksi dan Simulasi dari jumlah populasi vektor dengan kendali 101 LAMPIRAN D
Biodata Penulis
xvii
107
Halaman ini sengaja dikosongkan.
DAFTAR TABEL
Tabel 4.1 Input Parameter . . . . . . . . . . . . . . . . . . . . . . 68 Tabel 4.2 Input Parameter . . . . . . . . . . . . . . . . . . . . . . 70 Tabel 4.3 Input Parameter . . . . . . . . . . . . . . . . . . . . . . 88
xxi
Halaman ini sengaja dikosongkan.
DAFTAR GAMBAR
Gambar 2.1 Bifurkasi Maju . . . . . . . . . . . . . . . . . . . . . 10 Gambar 2.2 Bifurkasi Mundur . . . . . . . . . . . . . . . . . . . 10 Gambar 4.1 Diagram Kompartemen Dari Model Interaksi Dinamis . . . . . . . . . . . . . . . . . . . Gambar 4.2 Kurva Bifurkasi Mundur . . . . . . . . . . . . . Gambar 4.3 Kurva Bifurkasi Maju . . . . . . . . . . . . . . . Gambar 4.4 Simulasi dari populasi individu yang terinfeksi dengan kendali . . . . . . . . . . . . . Gambar 4.5 Simulasi dari jumlah populasi vektor dengan kendali . . . . . . . . . . . . . . . . . . . . .
xix
23 69 71 89 90
Halaman ini sengaja dikosongkan.
Daftar Simbol Sh (t) Ih (t) Rh (t) Sv (t) Iv (t) Nh (t) Nv (t) b1 b2 β1 β2 β3 µh µv γh δh δv R0 Ef E1
Individu yang rentan terhadap penyakit pada waktu t. Individu yang telah terinfeksi penyakit pada waktu t. Individu yang telah sembuh dari penyakit pada waktu t. Vektor yang rentan terhadap penyakit pada waktu t. Vektor yang telah terinfeksi penyakit pada waktu t. Jumlah populasi individu pada waktu t. Jumlah populasi vektor pada waktu t. angka kelahiran populasi individu yang rentan. angka kelahiran populasi vektor yang rentan. angka kejadian penularan melalui kontak langsung dengan individu yang terinfeksi. angka kejadian penularan melalui kontak dengan vektor yang terinfeksi. angka individu yang terinfeksi menjadikan vektor rentan terinfeksi. angka kematian natural dari populasi individu. angka kematian natural dari populasi vektor. angka kesembuhan pada populasi individu yang terinfeksi. angka kematian individu yang disebabkan penyakit. angka kematian dari vektor yang terinfeksi yang disebabkan penyakit. Bilangan reproduksi dasar. Titik kesetimbangan bebas penyakit. Titik kesetimbangan endemik. xxiii
Daftar Simbol
J λ H J(u1 , u2 , u3 ) I = [t0 , tf ] λ1 , λ2 , λ3 λ4 , λ5 u1 , u2 , u3 δ
Matrik Jacobian. Nilai eigen. Fungsi Hamiltonian. Fungsi Objektif. Periode waktu yang direncanakan. Variabel ko-keadaan (costate). Variabel pengendali. Variasi.
xxv
BAB I PENDAHULUAN
Pada bab ini dijelaskan hal-hal yang melatarbelakangi munculnya permasalahan yang dibahas dalam Tugas Akhir ini. Kemudian permasalahan tersebut disusun kedalam suatu rumusan masalah. Selanjutnya dijabarkan juga batasan masalah untuk mendapatkan tujuan yang diinginkan serta manfaat yang dapat diperoleh. Adapun sistematika penulisan Tugas Akhir diuraikan pada bagian akhir bab ini. 1.1
Latar Belakang Penularan penyakit dapat terjadi pada manusia melalui vektor (perantara) seperti serangga dikenal sebagai ”arthropod borne disease” atau kadang-kadang disebut juga dengan ”vector borne disease”. Penyakit yang tergolong ”arthropod borne disease” merupakan penyakit yang bersifat berbahaya dan dapat bersifat endemis serta dapat menimbulkan wabah dengan ancaman kematian [3]. Penyakit vektor-borne seperti demam berdarah, virus West Nile, virus encephalitis dan malaria merupakan infeksi menular ke manusia dan hewan lain yang disebabkan oleh arthropoda pemakan-darah. Arthropoda (serangga atau arakhnida) yang kebanyakan menyediakan vektor termasuk serangga pengisap darah seperti, nyamuk, kutu, dan lalat penggigit. Mayoritas penyakit bawaan vektor bertahan di alam dengan memanfaatkan hewan sebagai host (inang) vertebratanya, dan zoonosis. Untuk sejumlah kecil zoonosis, seperti malaria dan demam berdarah, manusia adalah host (inang) utama. Vektor pembawa patogen dari host (inang) 1
2 yang terinfeksi dan mengirimkannya untuk host (inang) perantara atau langsung ke host (inang) manusia [1]. Penyakit vector-borne penting bagi WHO (World Health Organization) khususnya wilayah asia bagian selatan-timur. Penyakit vector-borne adalah ancaman serius bagi kesehatan, yang mempengaruhi populasi yang paling rentan secara ekonomi. Dengue adalah salah satu penyakit menular yang paling cepat menyebar dari abad kedua puluh satu. Penyebarannya berubah ketika bergerak dari kota ke daerah pedesaan dan wilayah geografis baru akibat perubahan iklim. Lebih dari 1,3 miliar orang di kawasan ini beresiko malaria, seperti lebih dari 75% penduduk tinggal di daerah rawan malaria [5]. Dalam Tugas Akhir ini akan dilakukan analisa pada model penyakit vector-borne dengan menyelidiki adanya bifurkasi mundur. Selain itu, digunakan kendali optimal untuk meminimalkan jumlah host (inang) yang terinfeksi serta jumlah populasi vektor. Selanjutnya, akan ditentukan solusi yang optimal dari model tersebut. 1.2
Rumusan Masalah Permasalahan dalam Tugas Akhir ini adalah :
1. Bagaimana menentukan basic reproduction number, kestabilan dari setiap titik kesetimbangan, dan bifurkasi mundur ? 2. Bagaimana menentukan solusi yang optimal dari model tersebut ? 3. Bagaimana interpretasi hasil analisa dari model tersebut beserta simulasinya ? 1.3
Batasan Masalah Batasan masalah yang digunakan dalam Tugas Akhir ini antara lain:
3 1. Menerapkan Prinsip Minimum penyelesaian yang optimal.
Pontryagin
dalam
2. Model interaksi dinamik antara populasi manusia dan populasi vektor nyamuk adalah kombinasi dari dua model non linear dari populasi individu dan vektor [1]. 1.4
Tujuan Tujuan yang dicapai dari penulisan Tugas Akhir ini antara lain: 1. Menentukan basic reproduction number, kestabilan dari setiap titik kesetimbangan, dan bifurkasi mundur. 2. Menentukan solusi yang optimal dari model tersebut. 3. Mengintrepetasikan hasil analisa dari model tersebut beserta simulasinya. 1.5
Manfaat Manfaat dari Tugas Akhir ini antara lain:
1. Membantu mempelajari dampak dari ditentukannya basic reproduction number, kestabilan dari setiap titik kesetimbangan, dan bifurkasi mundur pada model penyakit vektor-borne. 2. Membantu menentukan cara yang efektif untuk mengendalikan pencegahan penyakit vektor-borne. 1.6
Sistematika Penulisan Penulisan Tugas Akhir ini disusun dalam lima bab, yaitu:
1. BAB I PENDAHULUAN Pada bab ini berisi tentang gambaran umum dari penulisan Tugas Akhir yang meliputi latar belakang, rumusan masalah, batasan masalah, tujuan, manfaat, dan sistematika penulisan.
4 2. BAB II TINJAUAN PUSTAKA Pada bab ini berisi tentang materi-materi yang mendukung Tugas Akhir ini, antara lain penyakit vektor-borne, sistem kompartemen, bilangan reproduksi dasar, kestabilan titik tetap, bifurkasi, teori kendali optimal, prinsip minimum pontryagin, dan metode beda hingga. 3. BAB III METODE PENELITIAN Pada bab ini dibahas tentang langkah langkah dan metode yang digunakan untuk menyelesaikan tugas akhir ini. 4. BAB IV ANALISIS DAN PEMBAHASAN Pada bab ini akan menguraikan bagaimana memperoleh daerah penyelesaian model, kestabilan lokal di setiap titik kesetimbangan, analisa bifurkasi berdasarkan bilangan reproduksi dasar, penerapan prinsip minimum pontryagin dan hamiltonian untuk mencari kendali optimal, mencari solusi numerik dengan menggunakan metode beda hingga dan simulasi dari model tersebut 5. BAB V PENUTUP Bab ini berisi kesimpulan akhir yang diperoleh dari Tugas Akhir serta saran untuk pengembangan penelitian selanjutnya.
BAB II TINJAUAN PUSTAKA
Pada bab ini akan dijelaskan mengenai tinjauan pustaka yang menjadi dasar materi dalam penyusunan tugas akhir serta menunjang metode metode yang digunakan dalam pembahasan tugas akhir ini. 2.1
Penyakit Vector-Borne Vektor adalah organisme yang mengirimkan patogen dan parasit dari satu orang yang terinfeksi (atau hewan) yang lain, menyebabkan penyakit serius pada populasi manusia [5]. Penyakit vector-borne adalah penyakit yang disebabkan oleh patogen dan parasit pada populasi manusia yang paling sering ditemukan di daerah tropis dan tempat di mana akses terhadap air minum dan sanitasi yang sehat terhambat. Penyakit vector-borne paling mematikan adalah malaria yang menyebabkan sekitar 660.000 kematian pada tahun 2010, sebagian besar adalah anak-anak Afrika. Namun, penyakit vector-borne paling cepat berkembang di dunia adalah demam berdarah, dengan peningkatan 30 kali lipat dalam kejadian penyakit selama 50 tahun terakhir [4]. Terjadinya penyebaran penyakit vektor-borne adalah melalui kontak langsung yaitu Artropoda secara langsung memindahkan penyakit atau investasi dari satu orang ke orang lain melalui kontak langsung, melalui transmisi secara mekanis yaitu Artropoda sebagai vektor mekanis yang membawa agen penyakit dari manusia yang berasal dari tinja, darah, ulkus superficial, atau eksudat, serta melalui transmisi secara biologis yaitu agen penyakit mengalami perubahan siklus 5
6 dengan atau tanpa multiplikasi di dalam tubuh artropoda. Penyakit yang ditularkan oleh vektor (nyamuk) adalah Malaria, filarial, yellow fever, ensefalitis, dan DHF. Vektor (nyamuk) dapat mudah berkembang pesat jika kondisi lingkungan mendukung. Beberapa kondisi yang mengakibatkan vektor (nyamuk) dapat berkembang pesat sebagai berikut [11]: a. Perubahan lingkungan fisik seperti pertambangan, industri dan pembangunan perumahan yang mengakibatkan berkembangbiaknya vektor(nyamuk) penyakit. b. Sistem drainase permukiman dan perkotaan yang tidak memenuhi syarat sehingga menjadi tempat perindukan vektor. c. Sistem pengelolaan sampah yang belum memenuhi syarat sehingga sampah menjadi sarang vektor(nyamuk). 2.2
Sistem Kompartemen Sistem kompartemen merupakan sebuah susunan kerja atau proses yang menunjukkan aliran individu dari satu kompartemen ke kompartemen lainnya seperti saat individu tersebut sehat, tertular penyakit atau sembuh dari penyakit [2]. 2.3
Bilangan Reproduksi Dasar (R0 ) Bilangan reproduksi dasar (Basic Reproduction Number) atau biasa disebut R0 adalah suatu parameter yang digunakan untuk mengetahui tingkat penyebaran suatu penyakit. Bilangan reproduksi dasar adalah bilangan yang menunjukkan jumlah individu rentan yang dapat menderita penyakit disebabkan oleh satu individu infeksi. Namun
7 adapula yang mengartikan bilangan yang menyatakan banyaknya rata-rata individu infected sekunder akibat tertular individu infected primer yang berlangsung dalam populasi individu rentan penyakit. Untuk menentukan bilangan reproduksi dasar, digunakan metode Driessche dan Watmough [12]. Dengan mengasumsikan bahwa populasi dapat dikelompokkan ke dalam n kompartemen. Diberikan x = (x1 , ..., xn )t , dengan xi ≥ 0 adalah bilangan dari individu pada masing-masing kompartemen sehingga kompartemen m < n pertama sesuai dengan individu terinfeksi. Diberikan Xs adalah himpunan dari semua titik kesetimbangan bebas penyakit. Didefinisikan Xs = {x ≥ 0|xi = 0, i = 1, ..., m} Selanjutnya untuk menghitung R0 , penting untuk membedakan infeksi baru dari semua perubahan dalam populasi. Untuk itu didefinisikan Fi (x) adalah laju dari kemunculan infeksi baru pada kompartemen i,Vi− (x) adalah laju dari perpindahan individu keluar dari kompartemen i, dan Vi+ (x) adalah laju dari perpindahan individu masuk ke kompartemen i. Model penyebaran penyakit terdiri dari kondisi awal non-negatif dengan persamaan sistem berikut: x˙i = fi = Fi (x) − Vi (x), i = 1, ..., n, dengan Vi = Vi− − Vi+ dan memenuhi asumsi-asumsi berikut: a. Jika x ≥ 0, maka Fi , Vi− , Vi+ ≥ 0 untuk i=1,...,n. Karena masing-masing merepresentasikan perpindahan langsung dari individu, maka semua non-negatif. b. Jika xi = 0, maka Vi− = 0. Secara khusus, jika x ∈ Xs , maka Vi− = 0 untuk i=1,...,n. Hal ini berarti jika sebuah kompartemen kosong, maka tidak ada perpindahan dari
8 individu yang keluar dari kompartemen dikarenakan kematian atau infeksi. c. Fi = 0 jika i > m. Hal ini mengakibatkan timbulnya penyakit adalah nol. d. Jika x ∈ Xs , maka Fi = 0 dan Vi+ = 0 untuk i=1,...,m. Jika populasi bebas dari penyakit, maka populasi akan tetap bebas dari penyakit (tidak ada infeksi). e. Jika F(x) menuju ke nol, maka semua nilai eigen dari Df (x0 ) mempunyai bagian real negatif. Hal ini berdasarkan turunan dari f di dekat titik kesetimbangan bebas penyakit (DFE), didefinisikan DFE dari f adalah penyelesaian kestabilan lokal dari titik kesetimbangan bebas penyakit, dengan f terbatas ke Xs . Jika populasi ada di sekitar DFE, maka populasi akan kembali ke DFE menurut linearisasi sistem : ∂fi x˙ = Df (x0 )(x − x0 ), dengan Dfi (x0 ) = ∂xi x=0
adalah matriks Jacobian yang dihitung di sekitar x = x0 . Oleh karena itu, asumsi ini mengakibatkan DFE stabil. Didefinisikan K = F V −1 sebagai next generation matriks dan ∂Vi (x0 ) ∂Fi (x0 ) −1 ,V = , R0 = ρ(F V ) dengan F = ∂xi ∂xi 1 ≤ i, j ≤ m, dan ρ(A) adalah nilai eigen yang dominan dari matriks A [6]. Jika model hanya mempunyai dua titik kesetimbangan yaitu titik kesetimbangan bebas penyakit dan titik kesetimbangan endemik, maka tidak terjadi endemik jika R0 ≤ 1 dan terjadi endemik jika R0 ≥ 1 [9].
9 2.4
Kestabilan Titik Tetap Pandang persamaan diferensial sebagai berikut : dx = f (x, y) dt dy = g(x, y) (2.1) dt Sebuah titik (x0 , y0 ) merupakan titik kesetimbangan dari persamaan (2.1) jika memenuhi f (x0 , y0 ) = 0 dan g(x0 , y0 ) = 0. Karena turunan suatu konstanta sama dengan nol, maka sepasang fungsi konstan [6]. x(t) ≡ x0 dan y(t) ≡ y0 adalah penyelesaian kesetimbangan dari persamaan (2.1) untuk semua t. 2.5
Stabil Asimtotik Lokal Kestabilan asimtotis lokal pada titik keseimbangan ditentukan oleh tanda pada bagian real dari akar-akar karakteristik sistem. Teorema 2.1 Titik setimbang (x0 , y0 ) stabil asimtotis jika dan hanya jika nilai karakteristik dari ∂f ∂f (x0 , y0 ) (x0 , y0 ) ∂x ∂y M atriks J = ∂g ∂g (x0 , y0 ) (x0 , y0 ) ∂x ∂y mempunyai tanda negatif pada bagian realnya dan tidak stabil jika sedikitnya satu dari nilai karakteristik mempunyai tanda positif pada bagian realnya. Analisis kestabilan dilakukan untuk mengetahui laju penyebaran suatu penyakit. Analisis ini dilakukan pada titik setimbang bebas penyakit (Disease Free Equilibrium) dan titik setimbang endemik (Endemic Equilibrium).
10 2.6
Bifurkasi
Pada sistem dinamik non linear sering dijumpai kestabilan di sekitar titik kesetimbangan suatu sistem persamaan yang menunjukkan fenomena bifurkasi. Bifurkasi secara umum adalah perubahan kualitatif yang meliputi perubahan stabilitas dan perubahan banyaknya titik kesetimbangan karena perubahan nilai-nilai parameter. Dalam epidemiologi, fenomena bifurkasi berhubungan dengan parameter ambang batas, yang sering disebut bilangan reproduksi dasar dan biasanya disimbolkan dengan R0 [7]. Ada dua jenis bifurkasi dalam model penyebaran penyakit menular yaitu bifurkasi maju dan bifurkasi mundur. Eksistensi bifurkasi maju dan mundur pada model penyebaran penyakit ditunjukkan oleh diagram bifurkasi pada Gambar 2.1 dan Gambar 2.2 dengan merupakan parameter bifurkasi dan merupakan populasi individu yang terinfeksi penyakit. Fenomena bifurkasi maju terjadi pada saat R0 > 1 dimana hanya ada satu titik kesetimbangan endemik. Sedangkan fenomena bifurkasi mundur terjadi pada saat R0 < 1 mempunyai dua titik kesetimbangan endemik [7].
Gambar 2.1 Bifurkasi Maju
Gambar 2.2 Bifurkasi Mundur
11 2.7
Teori Kendali Optimal Pada prinsipnya, tujuan dari pengendali optimal adalah menentukan signal atau kendali yang akan diproses dalam sistem dinamik dan memenuhi beberapa konstrain, dengan tujuan memaksimumkan atau meminimumkan fungsi tujuan (J) yang sesuai [10]. Adapun formulasi masalah kendali optimal terdiri dari: a. Mendiskripsikan (model).
secara
matematis
suatu
sistem
b. Menentukan fungsi objektif (performance index). c. Menentukan kendala dan kondisi batas yang harus dipenuhi. Secara umum, masalah kendali optimal diformulasikan sebagai berikut: misalkan suatu sistem dinamik diberikan oleh persamaan: x˙ = f (x (t), u(t), t) (2.2) dengan keadaan awal x (t0 ) = x 0 dan keadaan akhir x (tf ) = x f serta u(t) yang menyatakan pengendali keadaan pada waktu t. Dalam hal ini, masalah kendali optimal adalah mencari pengendali optimal u ∗ (t) yang memenuhi persamaan keadaan (state) dengan syarat nilai J berikut ini: Z
tf
J = S (x (tf ), tf ) +
V (x (t), u(t), t)
(2.3)
t0
adalah minimum atau maksimum. Bentuk umum persamaan J di atas disebut fungsi tujuan bentuk Bolza dengan S adalah bentuk Mayer dan V adalah bentuk Lagrange. Dengan kondisi sistem yaitu waktu akhir tetap atau bebas dan keadaan (state) akhir seluruhnya atau sebagian bebas atau tetap.
12 2.8
Prinsip Minimum Pontryagin Prinsip Minimum Pontryagin merupakan salah satu cara dalam menyelesaikan masalah kendali optimal dengan kendala yang terbatas. Metode tersebut digunakan untuk memperoleh kendali terbaik pada sistem dinamik dari state awal hingga akhir, yaitu dengan memininumkan fungsi objektif yang ingin dicapai. Hal ini dikembangkan oleh L. S. Pontryagin pada tahun 1950. Oleh karena itu, prinsip ini disebut sebagai Prinsip Minimum Pontryagin. Prinsip ini menyatakan secara informal bahwa persamaan Hamiltonian akan diminimumkan sepanjang U yang merupakan himpunan dari semua kendali [13]. Hasilnya juga dapat dinamakan Prinsip Maksimum Pontryagin yaitu dengan mengalikan (-1) pada performance index. Penyelesaian masalah kendali optimal dengan menggunakan metode tidak langsung dilakukan dengan menyelesaikan kondisi perlu kendali optimal. Berdasarkan Prinsip Minimum Pontryagin, kondisi perlu dari masalah kendali optimal yang harus diselesaikan adalah persamaan stasioner, persamaan state, dan persamaan costate serta kondisi transversality [9]. Dengan memperhatikan persamaan keadaan dan fungsi tujuan yang telah diberikan pada persamaan (2.2) dan (2.3), langkah dalam menyelesaikan masalah kendali optimal adalah sebagai berikut [8] : a. Langkah 1 Bentuk fungsi Hamiltonian yang disimbolkan dengan H, yaitu: H = H(x(t), u(t), λ(t), t) = V (x(t), u(t), t) + λ0 (t) f (x(t), u(t), t) dengan tanda ”0” menyatakan suatu transpose.
13 b. Langkah 2 Memaksimumkan H terhadap u(t) dengan cara: ∂H =0 ∂u(t) sehingga diperoleh kondisi stasioner u∗ (t). c. Langkah 3 Dengan menggunakan u∗ (t) yang telah dihasilkan pada langkah 2, akan didapatkan fungsi Hamiltonian baru yang optimal, H ∗ , yaitu: H ∗ = H(x∗ (t), u∗ (t), λ∗ (t), t) ≤ H(x∗ (t), u(t), λ∗ (t), t) d. Langkah 4 Selesaikan 2n persamaan, dengan n adalah jumlah variabel keadaan: x˙ ∗ (t) =
∂H ∗ ∂λ
dan persamaan costate yaitu: ∗
∂H ∗ λ˙ (t) = − ∂x
Dengan kondisi batas diberikan oleh keadaan awal dan keadaan akhir yang disebut kondisi transversality. Kondisi batas secara umum yaitu: 0 ∂S ∂S ∗ ∗ δtf + −λ H + δxf = 0 ∂t tf ∂x ∗ tf Dengan S adalah bentuk Mayer dari fungsi objektif, H adalah persamaan Hamiltonian, δ menunjukkan variasi dan tanda ∗ menunjukkan keadaan saat variabel pengendalinya stasioner.
14 e. Langkah 5 Substitusi hasil-hasil yang diperoleh pada langkah 4 ke dalam persamaan u∗ (t) pada langkah 2 untuk mendapatkan kendali optimal yang dicari. Dalam menentukan kondisi transversality yang sesuai, terdapat macam-macam kondisi batas, yaitu [8]: a. Fixed-final time and Fixed-final state system Artinya waktu akhir dan state saat waktu akhir telah ditentukan/diketahui. x(t0 ) = x0 , x(tf ) = xf b. Free-final time and Fixed-final state system Artinya waktu akhir belum ditentukan/tidak diketahui dan state saat waktu akhir telah ditentukan/diketahui. ∂S ∗ x(t0 ) = x0 , x(tf ) = xf , H + =0 ∂t tf c. Fixed-final time and Free-final state system Artinya waktu akhir telah ditentukan/diketahui sedangkan state saat waktu akhir belum diketahui/tidak ditentukan. ∂S ∗ x(t0 ) = x0 , λ (tf ) = ∂x ∗tf d. Free-final time and dependent free-final state system Artinya waktu akhir belum ditentukan/tidak diketahui dan state saat akhir belum ditentukan/tidak diketahui dan nilainya bergantung pada sesuatu. x(t0 ) = x0 , x(tf ) = θ f , 0 ∂S ∂S ∗ ˙ =0 H+ + − λ (t) θ(t) ∂t ∗ ∂x ∗ tf
15 e. Free-final time and independent free-final state system Artinya waktu akhir belum ditentukan/tidak diketahui dan state saat akhir belum ditentukan/tidak diketahui dan nilainya tidak bergantung pada sesuatu. δx(t0 ) = x0 ∂S ∂S ∗ H+ = 0, − λ (t) =0 ∂t ∗tf ∂x ∗tf 2.9
Metode Beda Hingga Jika u = u(x) maka turunan pertama dari u terhadap x didefinisikan du u(x + h) − u(x) = lim dx h→0 h dengan mensubtitusikan nilai x = x − h, maka du u(x) − u(x − h) = lim dx h→0 h dengan mensubtitusikan nilai x = x + h2 , maka u(x + h2 ) − u(x − h2 ) du = lim dx h→0 h Jika u = u(x) diekspansikan menurut deret taylor, maka diperoleh 1. Persamaan persamaan beda hingga maju adalah sebagai berikut: hdu h2 d2 u (x) + (x) + ... 1!dx 2! dx2 hdu u(x + h) − u(x) = (x) + o(h2 ) dx u(x + h) − u(x) du ≈ h dx u(x + h) = u(x) +
(2.4)
16 2. Persamaan persamaan beda hingga mundur adalah sebagai berikut: hdu h2 d2 u (x) + (x) + ... 1!dx 2! dx2 hdu u(x) − u(x − h) = (x) + o(h2 ) dx u(x) − u(x − h) du ≈ h dx u(x − h) = u(x) −
(2.5)
3. Persamaan persamaan beda hingga tengah jika persamaan (2.4) dikurangkan dengan persamaan (2.5), maka du + ... dx du u(x + h) − u(x − h) = 2h + o(h2 ) dx u(x + h) − u(x − h) du ≈ 2h dx u(x + h) − u(x − h) = 2h
Dengan tiga syarat batas, yaitu: 1. Syarat batas Diriclet, contoh: u(0) = 100 2. Syarat batas Neumann, contoh :
du (1) = 0 dx
3. Syarat batas Robbins, contoh: u(0) +
du (0) = 3 dx
Selain itu persamaan beda dapat dituliskan dalam bentuk lain yaitu sebagai berikut: diberikan u(x) = u(x = ih) = ui , maka
17
du ui+1 − ui = adalah beda maju dx h du ui − ui−1 = adalah beda mundur dx h du ui+1 − ui−1 = adalah beda tengah dx 2h
Halaman ini sengaja dikosongkan.
BAB III METODE PENELITIAN
Bab ini menguraikan metode yang akan digunakan dalam penelitian secara rinci. Metodologi penelitian yang digunakan berguna sebagai acuan sehingga penelitian ini dapat disusun secara sistematis. 3.1
Studi Literatur
Tahap ini merupakan tahap untuk melakukan identifikasi permasalahan, yaitu mencari referensi yang menunjang penelitian. Referensi bisa berupa tugas akhir, jurnal, buku, maupun artikel terkait. 3.2
Mengkaji Model Interaksi Dinamis
Untuk memahami model interaksi dinamik disusun asumsi-asumsi tertentu sehingga dapat dibuat model kompartemen dengan susceptible, infected, dan recovered untuk populasi manusia dan susceptible dan infected untuk populasi nyamuk. 3.3
Mencari Titik Kesetimbangan, Menentukan Bilangan Reproduksi Dasar dari Model, dan Bifurkasi
Dari model interaksi dinamis akan dicari titik kesetimbangan bebas penyakit (I = 0) dan titik kesetimbang -an endemik (I 6= 0) melalui substitusi persamaan model. Kemudian menentukan nilai bilangan reproduksi dasar (R0 ) dan selanjutnya menentukan kurva bifurkasi melalui nilai R0 . 19
20 3.4
Menentukan Formulasi Masalah Kendali Optimal Pada tahapan ini, ditentukan formulasi masalah kendali optimal yang meliputi menentukan model kendali optimal, fungsi objektif, dan kondisi syarat batas yang harus dipenuhi. 3.5
Menentukan Penyelesaian Kendali Optimal Pada tahap ini, dilakukan penyelesaian kendali optimal yang telah diformulasikan pada tahapan sebelumnya. Metode yang digunakan dalam penyelesaian masalah tersebut adalah Prinsip Minimum Pontryagin. Langkah-langkah yang dilakukan dalam tahap ini antara lain: 1. Membentuk fungsi Hamiltonian, 2. Menentukan persamaan state dan costate, 3. Menentukan kondisi batas yang harus dipenuhi, 4. Menentukan pengendali optimal. 3.6
Membuat Simulasi Dari Bifurkasi yang diperoleh akan dibuat simulasinya dan pada kendali optimal diselesaikan dengan menggunakan metode beda hingga semi implisit. Langkah-langkahnya adalah menyelesaikan model sistem dengan menggunakan metode beda hingga maju, dan fungsi adjoin yang diperoleh diselesaikan dengan menggunakan metode beda hingga mundur karena adanya kondisi transversality. Selanjutnya, dapat diperoleh penyelesaian model sistem secara numerik. 3.7
Kesimpulan dan Saran Setelah dilakukan analisis dan pembahasan maka dapat ditarik suatu kesimpulan dan saran sebagai masukan untuk pengembangan penelitian lebih lanjut.
BAB IV ANALISIS DAN PEMBAHASAN
Pada bab ini, akan dibahas tentang daerah penyelesaian model, titik kesetimbangan bebas penyakit,titik kesetimbangan endemik, kemudian akan dicari kestabilan lokal dari setiap titik kesetimbangan tersebut, dan bilangan reproduksi dasar, kemudian menentukan bifurkasinya (bifurkasi mundur) berdasarkan nilai bilangan reproduksi dasar, dan menentukan solusi yang optimal dari model dengan menerapkan Prinsip Minimum Pontryagin. Selanjutnya akan ditentukan penyelesaian solusi numerik dari model dan mensimulasikannya dengan menggunakan MATLAB. 4.1
Deskripsi Model dan Asumsi Model interaksi dinamis yang akan dibahas pada Tugas Akhir ini memiliki asumsi sebagai berikut: 1. Individu dikelompokkan menjadi tiga kelompok, yaitu susceptible (Sh (t)) yaitu individu yang rentan terhadap penyakit pada saat t , infected (Ih (t)) adalah individu yang terinfeksi penyakit pada saat t, dan recovered (Rh (t)) adalah individu yang sembuh dari penyakit pada saat t sedangkan untuk vektor (nyamuk) dikelompokkan menjadi dua kelompok, yaitu susceptible (Sv (t)) adalah vektor yang rentan terhadap penyakit pada saat t , dan infected (Iv (t)) adalah vektor yang terinfeksi penyakit pada saat t. Vektor yang sembuh akan langsung mati, sehingga tidak terjadi penularan. Jumlah populasi individu dinyatakan sebagai Nh (t) dengan Nh (t) = 21
22 Sh (t) + Ih (t) + Rh (t), sedangkan untuk jumlah populasi vektor dinyatakan sebagai Nv (t) = Sv (t) + Iv (t). 2. Model mendiskripsikan interaksi antara individu (rentan) yang berkaitan dengan kematian yang disebabkan karena kontak dengan penyakit maupun vektor yang terifeksi penyakit. Individu (rentan) terhadap penyakit dapat terinfeksi melalui dua cara, yaitu melakukan kontak langsung dengan penyakit β1 dan melakukan kontak dengan vektor yang terinfeksi β2 . 3. Berikut merupakan definisi parameter-parameter yang terdapat dalam model dinamik, yaitu a. b1 menyatakan angka kelahiran populasi individu yang rentan dan b2 menyatakan angka kelahiran populasi vektor yang rentan. b. β3 menyatakan angka individu yang terinfeksi menjadikan vektor rentan terinfeksi. c. µh menyatakan angka kematian natural dari populasi individu dan µv menyatakan angka kematian natural dari populasi vektor. d. γh menyatakan angka kesembuhan pada populasi individu yang terinfeksi, δh menyatakan angka kematian individu yang disebabkan penyakit, dan δv menyatakan angka kematian dari vektor yang terinfeksi yang disebabkan penyakit. Dari asumsi di atas, dapat digambarkan diagram kompartemen dari model interaksi dinamis sebagai berikut:
23
Gambar 4.1: Diagram Kompartemen Dari Model Interaksi Dinamis Dari Gambar 4.1, diperoleh model interaksi dinamis sebagai berikut: 1. Besarnya laju populasi individu yang rentan terhadap penyakit (susceptible) dipengaruhi oleh angka kelahiran populasi individu yang rentan, sedangkan populasi akan menurun dengan adanya beberapa kejadian penularan penyakit, seperti populasi individu yang terinfeksi penyakit karena melakukan kontak langsung dengan individu yang terinfeksi, populasi individu yg terinfeksi penyakit karena melakukan kontak dengan vektor yang terinfeksi, dan kematian alami dari individu. β1 Sh Ih β2 Sh Iv dSh = b1 − − − µh S h dt Nh Nh
24 2. Besarnya laju populasi individu yang terinfeksi penyakit (infected) akan bertambah saat terdapat populasi individu yang terinfeksi penyakit akibat kontak langsung dengan individu yang terinfeksi maupun dengan vektor yang terinfeksi dan populasi akan menurun dengan adanya kejadian individu yang terinfeksi namun telah sembuh dan juga menurun karena kematian yang disebabkan karena terinfeksi penyakit maupun secara alami. dIh β1 Sh Ih β2 Sh Iv = + − γh Ih − δh Ih − µh Ih dt Nh Nh 3. Besarnya laju populasi individu yang sembuh dari penyakit (recovered) bergantung pada jumlah individu yang terinfeksi namun telah sembuh dan populasi akan berkurang saat terdapat kejadian kematian individu yang sembuh. dRh = γh Ih − µh Rh dt 4. Besarnya laju populasi vektor yang rentan terhadap penyakit (susceptible) bergantung pada angka kelahiran populasi vektor yang rentan dan populasi akan berkurang saat angka populasi individu yang terinfeksi yang menjadikan vektor rentan terinfeksi dan kematian populasi vektor. dSv β3 Sv Ih = b2 − − µv Sv dt Nh 5. Besarnya laju populasi vektor yang terinfeksi penyakit (infected) bergantung pada angka populasi individu yang terinfeksi yang menjadikan vektor rentan terinfeksi
25 dan populasi akan berkurang akibat angka kematian dari vektor yang terinfeksi yang disebabkan penyakit dan kematian alami yang terinfeksi. dIv β3 Sv Ih = − δ v Iv − µ v Iv dt Nh
Dari penjelasan diatas, maka sistem persamaan dapat ditulis sebagai berikut [1]: dSh dt dIh dt dRh dt dSv dt dIv dt
β1 Sh Ih β2 Sh Iv − − µh Sh (4.1) Nh Nh β1 Sh Ih β2 Sh Iv + − γh Ih − δh Ih − µh Ih (4.2) Nh Nh
= b1 − =
= γh Ih − µh R h β3 Sv Ih − µv Sv Nh β3 Sv Ih − δv Iv − µv Iv Nh
(4.3)
= b2 −
(4.4)
=
(4.5)
dengan kondisi awal: Sh (0) > 0, Ih (0) > 0, Rh (0) > 0, Sv (0) > 0, Iv (0) > 0 Diketahui bahwa jumlah populasi individu dinyatakan sebagai Nh (t) dengan Nh (t) = Sh (t) + Ih (t) + Rh (t) sehingga dapat dituliskan Rh (t) = Nh (t) − Sh (t) − Ih (t). Supaya memudahkan perhitungan analisis selanjutnya, persamaan h (4.3) digantikan oleh persamaan dN dt = b1 − µh Nh − δh Ih dengan Nh = Sh + Ih + Rh , kemudian turunkan persamaan
26 tersebut terhadap t, sehingga: Nh0 (t) = Sh0 (t) + Ih0 (t) + Rh0 (t) β1 Sh Ih β2 Sh Iv = b1 − − − µh Sh Nh Nh β1 Sh Ih β2 Sh Iv + + − γh Ih − δh Ih − µh Ih Nh Nh + (γh Ih − µh Rh ) = b1 − µh Sh − δh Ih − µh Ih − µh Rh = b1 − µh Sh − δh Ih − µh Ih − µh (Nh − Sh − Ih ) = b1 − µh Sh − δh Ih − µh Ih − µh Nh + µh Sh + µh Ih = b1 − µh Nh − δh Ih
(4.6)
Maka sistem persamaan baru dari model menjadi : dSh dt dIh dt dNh dt dSv dt dIv dt
β1 Sh Ih β2 Sh Iv − − µh S h (4.7) Nh Nh β1 Sh Ih β2 Sh Iv + − γh Ih − δh Ih − µh Ih (4.8) Nh Nh
= b1 − =
= b1 − µh Nh − δh Ih β3 Sv Ih − µv Sv Nh β3 Sv Ih − δv Iv − µv Iv Nh
(4.9)
= b2 −
(4.10)
=
(4.11)
Dengan kondisi awal: Sh (0) > 0, Ih (0) > 0, Nh (0) > 0, Sv (0) > 0, Iv (0) > 0 (4.12) 4.2
Daerah Penyelesaian Model Dalam sistem persamaan baru dari model, diketahui h bahwa dN dt = b1 − µh Nh − δh Ih selanjutnya dengan cara yang
27 v sama akan dicari nilai dN dt . Karena diketahui bahwa Nv = Sv + Iv Kemudian turunkan persamaan tersebut terhadap t, diperoleh :
Nv0 (t) = Sv0 (t) + Iv0 (t) β3 Sv Ih β 3 S v Ih = b2 − − µv Sv + − δv Iv − µv Iv Nh Nh = b2 − µv Sv − δv Iv − µv Iv selanjutnya substitusikan Sv = Nv − Iv Nv0 (t) = b2 − µv (Nv − Iv ) − δv Iv − µv Iv = b2 − µv Nv + µv Iv − δv Iv − µv Iv = b2 − µv Nv − δv Iv
(4.13)
Selanjutnya akan dicari nilai Nh (t) dan Nv (t) melalui persamaan (4.6) dan (4.13) dNh = b1 − µh Nh − δh Ih dt ≤ b1 − µh Nh Sehingga persamaan diferensialnya dapat dituliskan dNh = b1 − µh Nh dt Maka dapat diselesaikan dengan menggunakan persamaan diferensial peubah terpisah. dNh = b1 − µh Nh dt dNh = dt b1 − µh Nh
(4.14)
Selanjutnya persamaan (4.14) akan diintegralkan di kedua ruasnya
28 Z
dNh = b1 − µh Nh
Z dt
Misal: u = b1 − µh Nh
(4.15)
maka dNh = − µduh , sehingga Z Z du − = dt uµh 1 1 ⇔ − ln|u| = t − ln|c| µh µh 1
⇔ ln|u| = −µh t + µh ln|c| µh 1
⇔ ln|u| − µh ln|c| µh = −µh t |u| ⇔ ln 1 = −µh t µ |c| µh h |u| = −µh t ⇔ ln |c| u ⇔ = e−µh t c ⇔ u = ce−µh t Selanjutnya substitusikan nilai u ke persamaan (4.15), sehingga b1 − µh Nh = ce−µh t µh Nh = b1 − ce−µh t b1 − ce−µh t µh b1 − ce−µh t lim Nh (t) = lim t→∞ t→∞ µh b1 = µh Nh (t) =
29 Karena itu, maka 0 < Nh (t) ≤ Sedangkan
b1 µh
dNv = b2 − µv Nv − δv Iv dt ≤ b2 − µv Nv Sehingga persamaan diferensialnya dapat dituliskan dNv = b2 − µv Nv dt Maka dapat diselesaikan dengan menggunakan persamaan diferensial peubah terpisah. dNv = b2 − µv Nv dt dNv = dt b2 − µv Nv
(4.16)
Selanjutnya persamaan (4.16) akan diintegralkan di kedua ruasnya Z
dNv = b2 − µv Nv
Z dt
Misal: u = b2 − µv Nv maka du dNv = − µv
(4.17)
30 sehingga Z
Z du − = dt uµv 1 1 ⇔ − ln|u| = t − ln|c| µv µv 1
⇔ ln|u| = −µv t + µv ln|c| µv 1
⇔ ln|u| − µv ln|c| µv = −µv t |u| ⇔ ln 1 = −µv t µ |c| µv v |u| ⇔ ln = −µv t |c| u ⇔ = e−µv t c ⇔ u = ce−µv t Selanjutnya substitusikan nilai u ke persamaan (4.17) sehingga b2 − µv Nv = ce−µv t µv Nv = b2 − ce−µv t b2 − ce−µv t µv b2 − ce−µv t lim Nv (t) = lim t→∞ t→∞ µv b2 = µv Nv (t) =
Sehingga 0 < Nv (t) ≤ µb2v Jadi daerah yang mungkin untuk sistem persamaan (4.7) sampai (4.11) adalah Ω = b b 5 , (N ≤ 1 , N ≤ 2 ) (Sh , Ih , Nh , Sv , Iv ) ∈ R+ dikarenakan v h µh µv kondisi awal pada persamaan (4.12) bernilai positif dan pada
31 5 , maka Ω merupakan invarian positif. Selain itu, dN h dan R+ dt dN v ≤ 0 memenuhi Ω yang merupakan invarian positif dan dt t → ∞, sehingga dapat ditulis 0 < (Nh , Nv ) ≤ ( µb1h , µb2v ).
4.3
Titik Kesetimbangan Bebas Penyakit
Titik kesetimbang bebas penyakit adalah suatu keadaan tidak terjadi penyebaran penyakit menular dalam suatu populasi sehingga I = 0. Untuk memperoleh titik kesetimbangan bebas penyakit dengan menyatakan ruas kiri pada persamaan (4.7) sampai (4.11) bernilai nol kemudian mensubstitusikannya untuk memperoleh titik Ef = Sh0 , Ih0 , Nh0 , Sv0 , Iv0 . Karena pada kondisi tidak ada penyebaran penyakit menular, maka Ih0 = Iv0 = 0. Selanjutnya akan dicari nilai Sh0 , Nh0 dan Sv0 melalui persamaan (4.7),(4.9), dan (4.10) yang ruas kanannya bernilai nol kemudian mensubstitusikannya dengan Ih0 = Iv0 = 0. Menentukan nilai Sh0 dSh =0 dt β1 Sh Ih β2 Sh Iv − − µh Sh = 0 Nh Nh ⇔ b1 − 0 − 0 − µh Sh = 0
⇔ b1 −
⇔ b1 − µh Sh = 0 ⇔ µh Sh = b1 b1 ⇔ Sh = , maka µh b1 Sh0 = µh
(4.18)
32 Menentukan nilai Nh0 dNh =0 dt ⇔ b1 − µh Nh − δh Ih = 0 ⇔ b1 − µh Nh − 0 = 0 ⇔ b1 − µh Nh = 0 ⇔ µh Nh = b1 b1 , maka ⇔ Nh = µh b1 Nh0 = µh
(4.19)
Menentukan nilai Sv0 dSv =0 dt β3 Sv Ih − µv S v = 0 Nh ⇔ b2 − 0 − µv Sv = 0 ⇔ b2 −
⇔ b2 − µv Sv = 0 ⇔ µv Sv = b2 b2 ⇔ Sv = , maka µv b2 Sv0 = µv
(4.20)
Berdasarkan persamaan (4.18),(4.19),dan (4.20), diketahui bahwa Ih0 = Iv0 = 0, maka diperoleh titik kesetimbanganbebas b1 b2 b1 , 0, , , 0 . penyakit Ef = Sh0 , Ih0 , Nh0 , Sv0 , Iv0 = µh µh µv
33 4.4
Menentukan Bilangan Reproduksi Dasar
Dalam model epidemiologi, bilangan reproduksi dasar yang dilambangkan dengan R0 adalah konsep kunci dan didefinisikan sebagai jumlah rata-rata infeksi sekunder yang timbul dari individu yang terinfeksi primer yang masuk ke kelas susceptible selama periode infeksi susceptible[2]. Dengan menggunakan metode Driessche dan Watmough [12] akan ditentukan bilangan reproduksi dasar untuk itu didefinisikan sebagai berikut :
dIh dt dIv dt dNh dt dSh dt dSv dt
β1 Sh Ih β2 Sh Iv + − (γh + δh + µh )Ih Nh Nh β3 Sv Ih = − (δv + µv )Iv Nh =
= b1 − µh Nh − δh Ih
(4.21)
β1 Sh Ih β2 Sh Iv − − µh S h Nh Nh β3 Sv Ih = b2 − − µv Sv Nh = b1 −
Fi adalah laju kemunculan infeksi baru pada kompartemen i, Vi− adalah laju dari perpindahan individu keluar dari kompartemen i, Vi+ adalah laju dari perpindahan individu masuk ke dalam kompartemen i Vi = Vi− − Vi+ Untuk sistem persamaan (4.21) F dan V adalah:
34 F =
β1 Sh Ih Nh β3 Sv Ih Nh
0 0 0
Sh Iv − β2N + (γh + δh + µh )Ih h
(δv + µv )Iv −b1 + µh Nh + δh Ih dan V = Sh Ih Sh Iv −b1 + β1N + β2N + µh Sh h h S v Ih −b2 + β3N + µv S v h
Kompartemen yang terinfeksi adalah Ih ke Iv , sehingga perhatikan Ih dan Iv dengan mensubtitusikan nilai titik b1 b1 b2 kesetimbangan bebas penyakit Ef = , 0, , , 0 . µh µh µv
β1 Ih
β3 b2 µh µv b1 Ih F = 0 0 0
−β2 Iv + (γh + δh + µh )Ih
(δv + µv )Iv −b1 + µh Nh + δh Ih dan V = β S I β S I −b1 + 1Nhh h + 2Nhh v + µh Sh Sv Ih −b2 + β3N + µv Sv h
Sehingga diperoleh, F =
β 1 Ih β3 b2 µh µv b1 Ih
−β2 Iv + (γh + δh + µh )Ih
dan V =
(δv + µv )Iv
35 Selanjutnya akan dicari matriks F dan V yang didapat dari F dan V dengan cara berikut sebagai berikut: ∂V ∂F ∂F ∂V Ih
Ih
F =
∂FIv Ih
Ih
Ih
Iv ∂FIv Iv
Ih
Ih
dan V =
∂VIv Ih
Iv ∂VIv Iv
maka akan didapatkan matriks sebagai berikut [12] β1 0 γh + δ h + µh −β2 dan V = F = β3 b2 µh 0 0 δv + µv µv b1 dari F tersebut dapat dituliskan sebagai berikut 0 0 β1 0 + F = β3 b2 µh 0 0 0 µv b1 Selanjutnya dicari V −1 , sehingga diperoleh δ v + µv β2 1 V −1 = (δv + µv )(γh + δh + µh ) 0 γ h + δh + µ h 1 (γh + δh + µh ) = 0
V −1
β2 (δv + µv )(γh + δh + µh ) 1 (δv + µv )
Selanjutnya mencari nilai R0 yaitu R0 = ρ(F V −1 ) sehingga a. untuk
β1 0
F =
0
0
36 maka 1 0 (γh + δh + µh ) 0 0
F V −1 =
β1 0
β1 (γh + δh + µh ) = 0
β2 (δv + µv )(γh + δh + µh ) 1 (δv + µv )
β1 β2 (δv + µv )(γh + δh + µh ) 0
Jadi nilai eigen dari matriks next generation diperoleh dengan menyelesaikan det (F V −1 − λI) = 0, sehingga didapat β1 β2 β1 − λ (γh + δh + µh ) (δv + µv )(γh + δh + µh ) det =0 0 −λ
β1 − λ (−λ) = 0 (γh + δh + µh ) β1 − λ + λ2 = 0 (γh + δh + µh ) β1 λ − +λ =0 (γh + δh + µh )
Dari perhitungan diatas bisa didapatkan nilai eigen untuk λ = 0 dan didapat λ=
β1 (γh + δh + µh )
(4.22)
37 b. Untuk F =
0
0
β3 b2 µh µv b1
0
maka F V −1 =
" =
0
0
β3 b 2 µ h µv b1
0
1
β2
(γh + δh + µh )
(δv + µv )(γh + δh + µh )
0
1
0
(δv + µv ) 0
β3 b2 µh
β2 β3 b2 µh
(µv b1 )(γh + δh + µh )
(µv b1 )(δv + µv )(γh + δh + µh )
#
Jadi nilai eigen dari matriks next generation diperoleh dengan menyelesaikan det (F V −1 − λI) = 0, sehingga didapat 0−λ 0 det =0 β β b µ β b µ (µv b1 )(γ3h 2+ δhh + µh ) (µv b1 )(δv +2µv3)(γ2 hh+ δh + µh ) − λ
β2 β3 b2 µh − λ (−λ) = 0 (µv b1 )(δv + µv )(γh + δh + µh ) β2 β3 b2 µh − λ + λ2 = 0 (µv b1 )(δv + µv )(γh + δh + µh ) β2 β3 b2 µh λ − +λ =0 (µv b1 )(δv + µv )(γh + δh + µh )
Dari perhitungan diatas bisa didapatkan nilai eigen untuk λ = 0 dan didapat β2 β3 b2 µh λ= (4.23) (µv b1 )(δv + µv )(γh + δh + µh ) Berdasarkan perolehan nilai eigen dari persamaan (4.22) dan (4.23), maka R0 diperoleh sebagai berikut β2 β3 b2 µh β1 R0 = + (µv b1 )(δv + µv )(γh + δh + µh ) (γh + δh + µh )
38 4.5
Titik Kesetimbangan Endemik Titik Kesetimbangan Endemik digunakan untuk menunjukkan adanya penyebaran penyakit pada suatu populasi sehingga Ih∗ 6= Iv∗ 6= 0. Untuk memperoleh titik kesetimbangan endemik dengan menyatakan ruas kiri bernilai nol pada persamaan (4.7) dIh dNh dSv h sampai (4.11), sehingga dS dt = 0, dt = 0, dt = 0, dt = 0, v dan dI = 0. Kemudian mensubstitusikannya untuk dt memperoleh titik E1 = (Sh∗ , Ih∗ , Nh∗ , Sv∗ , Iv∗ ). Pertama-tama akan dicari nilai Nh∗ sehingga dNh =0 dt ⇔ b1 − µh Nh − δh Ih = 0 ⇔ µh Nh = b1 − δh Ih b1 − δh Ih∗ Nh∗ = µh
(4.24)
Selanjutnya akan dicari nilai Sh∗ dengan cara berikut dIh =0 dt β1 Sh Ih β2 Sh Iv ⇔ + − γh Ih − δh Ih − µh Ih = 0 Nh Nh β1 Ih + β2 Iv ⇔ Sh = (γh + δh + µh )Ih (4.25) Nh Selanjutnya dSh =0 dt β1 Sh Ih β2 Sh Iv − − µh Sh = 0 ⇔ b1 − Nh Nh β1 Ih + β2 Iv ⇔ b1 − µh Sh = Sh Nh
(4.26)
39 Dengan mensubtitusikan persamaan (4.25) ke persamaan (4.26) , sehingga b1 − µh Sh = (γh + δh + µh )Ih b1 − (γh + δh + µh )Ih = µh Sh b1 − (γh + δh + µh )Ih∗ Sh∗ = (4.27) µh Selanjutnya akan dicari nilai Sv∗ dengan cara berikut dSv β3 Sv Ih = b2 − − µ v Sv = 0 dt Nh β3 Sv Ih = b2 − µv Sv ⇔ Nh β3 Sv Ih ⇔ + µv Sv = b2 Nh β3 Ih + µv Nh ⇔ Sv = b2 Nh b2 Nh∗ Sv∗ = ∗ β3 Ih + µv Nh∗
(4.28)
Dengan mensubtitusikan persamaan (4.24) ke persamaan (4.28), sehingga
Sv∗ =
b2
b1 −δh Ih∗ µh
b −δ I ∗ β3 Ih∗ + µv 1 µhh h b2 b1 − b2 δh Ih∗ µh = µh β3 Ih∗ µh + µv (b1 − δh Ih∗ ) b2 (b1 − δh Ih∗ ) = β3 Ih∗ µh + µv (b1 − δh Ih∗ )
40 Selanjutnya akan dicari nilai Iv∗ , sehingga dIv =0 dt β 3 S v Ih ⇔ − δv Iv − µv Iv = 0 Nh β 3 S v Ih = Iv (δv + µv ) ⇔ Nh 1 −δh Ih ) β3 β3 Ih µb2h(b+µ Ih v (b1 −δh Ih ) ⇔ = Iv (δv + µv ) b1 −δ I h h
µh
⇔
β3 b2 Ih (b1 − δh Ih ) µh β3 Ih + µv (b1 − δh Ih )
µh b1 − δh Ih
= Iv (δv + µv ) Iv∗ =
β3 b2 Ih∗ µh (δv + µv ) µh β3 Ih∗ + µv (b1 − δh Ih∗ )
(4.29)
Misalkan K1 = (γh + δh + µh ) dan K2 = (δv + µv ), sehingga persamaan (4.27) menjadi Sh∗ =
b1 − K1 Ih∗ µh
(4.30)
dan persamaan (4.29) menjadi Iv∗ =
K2
β3 b2 Ih∗ µh µh β3 Ih∗ + µv (b1 − δh Ih∗ )
(4.31)
Karena Nilai Ih∗ 6= 0, maka selanjutnya akan dicari nilai g(Ih ) dengan cara mensubtitusikan nilai Sh∗ , Nh∗ , Sv∗ , dan Iv∗ pada persamaan (4.8) dan sisi kiri persamaan (4.8) disama
41 dengankan nol, sehingga dIh β1 Sh Ih β2 Sh Iv = + − γ h I h − δh I h − µ h I h dt Nh Nh β1 Sh∗ Ih∗ β2 Sh∗ Iv∗ + − (γh + δh + µh )Ih∗ ⇔0= Nh∗ Nh∗ b −K I ∗ β1 1 µh1 h Ih∗ ⇔0= b −δ I ∗ 1
h h
µh
β2 +
⇔0=
⇔0=
b1 −K1 Ih∗ µh
β3 b2 Ih∗ µh K2 (µh β3 Ih∗ +µv (b1 −δh Ih∗ ))
β1 (b1 − K1 Ih∗ )Ih∗ + (b1 − δh Ih∗ ) K2
⇔0=
b1 −δh Ih∗ µh
− K1 Ih∗
β2 (b1 − K1 Ih∗ ) (b1 − δh Ih∗ ) ! (β3 b2 Ih∗ µh ) − K1 Ih∗ µh β3 Ih∗ + µv (b1 − δh Ih∗ )
β1 (b1 − K1 Ih∗ )Ih∗ (K2 (µh β3 Ih∗ + µv (b1 − δh Ih∗ ))) (b1 − δh Ih∗ ) K2 µh β3 Ih∗ + µv (b1 − δh Ih∗ ) β2 (b1 − K1 Ih∗ )β3 b2 Ih∗ µh + (b1 − δh Ih∗ ) K2 µh β3 Ih∗ + µv (b1 − δh Ih∗ ) K1 Ih∗ (b1 − δh Ih∗ ) (K2 (µh β3 Ih∗ + µv (b1 − δh Ih∗ ))) − (b1 − δh Ih∗ ) K2 µh β3 Ih∗ + µv (b1 − δh Ih∗ ) (β1 b1 Ih∗ − β1 K1 Ih∗2 ) (K2 (µh β3 Ih∗ + µv (b1 − δh Ih∗ ))) (K2 b1 − K2 δh Ih∗ ) µh β3 Ih∗ + µv (b1 − δh Ih∗ ) (β2 β3 b1 b2 µh Ih∗ − β2 β3 b1 b2 K1 µh Ih∗2 ) (b1 − δh Ih∗ ) K2 µh β3 Ih∗ + µv (b1 − δh Ih∗ ) (K1 b1 Ih∗ − K1 δh Ih∗ ) (K2 (µh β3 Ih∗ + µv (b1 − δh Ih∗ ))) − (b1 − δh Ih∗ ) K2 µh β3 Ih∗ + µv (b1 − δh Ih∗ ) +
42 ⇔0=
β1 b1 K2 µh β3 Ih∗2 + K2 µv b21 β1 Ih∗ (b1 − δh Ih∗ ) K2 µh β3 Ih∗ + µv (b1 − δh Ih∗ ) −K2 µv b1 β1 δh Ih∗2 − β1 K1 K2 µh β3 Ih∗3 (b1 − δh Ih∗ ) K2 µh β3 Ih∗ + µv (b1 − δh Ih∗ ) −β1 K1 K2 µv b1 Ih∗2 + β1 K1 K2 µv δh Ih∗3 (b1 − δh Ih∗ ) K2 µh β3 Ih∗ + µv (b1 − δh Ih∗ ) +β2 b1 β3 b2 µh Ih∗ − β2 K1 µh b2 β3 Ih∗2 (b1 − δh Ih∗ ) K2 µh β3 Ih∗ + µv (b1 − δh Ih∗ ) −K1 K2 b1 µh β3 Ih∗2 − K1 K2 b21 µv Ih∗ (b1 − δh Ih∗ ) K2 µh β3 Ih∗ + µv (b1 − δh Ih∗ ) +K1 K2 b1 µv δh Ih∗2 + K1 K2 µh δh β3 Ih∗3 (b1 − δh Ih∗ ) K2 µh β3 Ih∗ + µv (b1 − δh Ih∗ ) +K1 K2 µv δh b1 Ih∗2 − K1 K2 µv δh2 Ih∗3 (b1 − δh Ih∗ ) K2 µh β3 Ih∗ + µv (b1 − δh Ih∗ )
Selanjutnya pembilang dan penyebut dibagi dengan −Ih∗ , sehingga −β1 b1 K2 µh β3 Ih∗ − K2 µv b21 β1 + K2 µv b1 β1 δh Ih∗ 0 = 1 ∗ ∗ ∗ −I ∗ (b1 − δh Ih ) K2 (µh β3 Ih + µv (b1 − δh Ih )) h
+
β1 K1 K2 µh β3 Ih∗2 + β1 K1 K2 µv b1 Ih∗ − β1 K1 K2 µv δh Ih∗2 1 (b1 − δh Ih∗ ) K2 µh β3 Ih∗ + µv (b1 − δh Ih∗ ) ∗ −I
−
β b β b µ − β2 K1 µh b2 β3 Ih∗ − K1 K2 b1 µh β3 Ih∗ 2 1 3 2 h 1 ∗ ) K µ β I ∗ + µ (b − δ I ∗ ) (b − δ I ∗ 1 2 3 v 1 h h h h h h −I
h
h
+
K K b2 µ − K1 K2 b1 µv δh Ih∗ − K1 K2 µh δh β3 Ih∗2 1 2 1 v 1 ∗ ∗ ∗ −I ∗ (b1 − δh Ih ) K2 µh β3 Ih + µv (b1 − δh Ih ) h
−
K1 K2 µv δh b1 Ih∗ − K1 K2 µv δh2 Ih∗2
1 −Ih∗
(b1 − δh Ih∗ ) K2 µh β3 Ih∗ + µv (b1 − δh Ih∗ )
43 0 =
+
−
−
−
−
−
(−β1 K1 K2 µv δh + β1 K1 K2 µh β3 )Ih∗2 1 (b1 − δh Ih∗ ) K2 µh β3 Ih∗ + µv (b1 − δh Ih∗ ) ∗ −Ih (−K1 K2 µh δh β3 + K1 K2 µv δh2 )Ih∗2 1 −Ih∗
1 −Ih∗
1 −Ih∗
(b1 − δh Ih∗ ) K2 µh β3 Ih∗ + µv (b1 − δh Ih∗ ) (β1 b1 K2 µh β3 − β1 b1 K2 µv δh )Ih∗
(b1 − δh Ih∗ ) K2 µh β3 Ih∗ + µv (b1 − δh Ih∗ ) (β1 b1 K2 µh β3 − β1 b1 K2 µv δh )Ih∗
(b1 − δh Ih∗ ) K2 µh β3 Ih∗ + µv (b1 − δh Ih∗ ) (−β1 K1 K2 µv b1 − β2 K1 µh b2 β3 )Ih∗
1 ∗ ) K µ β I ∗ + µ (b − δ I ∗ ) (b − δ I 2 v 1 1 3 h h h h h h −Ih∗ (−β K K µ b + b1 K1 K2 µv δh + b1 K1 K2 µv δh )Ih∗ 3 1 2 h 1 1 (b1 − δh Ih∗ ) K2 µh β3 Ih∗ + µv (b1 − δh Ih∗ ) ∗ −Ih
β b2 K µ − β2 β3 µh b1 b2 + K1 K2 µv b21 1 1 2 v 1 (b1 − δh Ih∗ ) K2 (µh β3 Ih∗ + µv (b1 − δh Ih∗ )) ∗ −Ih
Sehingga diperoleh g(Ih ) = (K1 K2 µh β1 β3 − K1 K2 µh β3 δh + K1 K2 µv δh2 K1 K2 β1 µv δh )Ih∗2 + (K1 K2 µv b1 β1 −2K1 K2 µv b1 δh + K1 K2 µh b1 β3 + K2 b1 β1 µv δh −K2 b1 β1 µh β3 + K1 µh b2 β2 β3 )Ih∗ + µv b21 K1 K2 µh b2 β2 β3 β1 −µv b21 K1 K2 − µv b21 K1 K2 µv b1 K1 K2 K1
44 =
(K1 K2 (µh β3 (β1 − δh ) + µv δh (δh − β1 ))) Ih∗2 + (K1 K2 (µv b1 β1 − 2µv b1 δh + µh b1 β3 )) Ih∗ + (K2 b1 β1 (µv δh − µh β3 ) + K1 µh b2 β2 β3 ) Ih∗ µh b2 β2 β3 β1 2 +µv b1 K1 K2 1 − + µv b1 K1 K2 K1
=
(K1 K2 (µh β3 (β1 − δh ) + µv δh (δh − β1 ))) Ih∗2 + (K1 K2 (µv b1 β1 − 2µv b1 δh + µh b1 β3 )) + (K2 b1 β1 (µv δh − µh β3 ) + K1 µh b2 β2 β3 ) Ih∗ +µv b21 K1 K2 (1 − R0 )
Sehingga diperoleh persamaan g(Ih ) sebagai berikut g(Ih ) = AIh∗2 + BIh∗ + C = 0 dengan A = K1 K2 (µh β3 (β1 − δh ) + µv δh (δh − β1 )) , B = K1 K2 (µv b1 β1 − 2µv b1 δh + µh b1 β3 ) + K2 b1 β1 (µv δh − µh β3 ) + K1 µh b2 β2 β3 , C = µv b21 K1 K2 (1 − R0 ), µh b2 β2 β3 β1 R0 = + µv b1 K1 K2 K1 4.6
Kestabilan Lokal Model Interaksi Dinamis Setelah diperoleh titik kesetimbangan maka dilakukan analisis kestabilan. Analisis kestabilan dilakukan untuk mengetahui laju penyebaran suatu penyakit. Analisis ini dilakukan pada titik setimbang bebas penyakit (Disease Free Equilibrium) dan titik setimbang endemik (Endemic Equilibrium). Model interaksi dinamis merupakan model persamaan yang tak linier, sehingga perlu dilakukan linierisasi terlebih
45 dahulu sebelum melakukan analisis kestabilan. Untuk melakukan linierisasi digunakan ekspansi deret Taylor, pada persamaan (4.7) sampai (4.11) sehingga dapat dituliskan sebagai berikut. dSh = A(Sh , Ih , Nh , Sv , Iv ) dt β1 Sh Ih β2 Sh Iv = b1 − − − µh S h Nh Nh dIh = B(Sh , Ih , Nh , Sv , Iv ) dt β1 Sh Ih β2 Sh Iv + − γh Ih − δh Ih − µh Ih = Nh Nh dNh = C(Sh , Ih , Nh , Sv , Iv ) (4.32) dt = b1 − µh Nh − δh Ih dSv = D(Sh , Ih , Nh , Sv , Iv ) dt β3 Sv Ih = b2 − − µv Sv Nh dIv = E(Sh , Ih , Nh , Sv , Iv ) dt β3 Sv Ih = − δv Iv − µv Iv Nh Dengan titik tetap (Sh0 , Ih0 , Nh0 , Sv0 , Iv0 ), maka dSh = A(Sh0 , Ih0 , Nh0 , Sv0 , Iv0 ) = 0 dt dIh = B(Sh0 , Ih0 , Nh0 , Sv0 , Iv0 ) = 0 dt dNh = C(Sh0 , Ih0 , Nh0 , Sv0 , Iv0 ) = 0 dt dSv = D(Sh0 , Ih0 , Nh0 , Sv0 , Iv0 ) = 0 dt
(4.33)
46 dIv = E(Sh0 , Ih0 , Nh0 , Sv0 , Iv0 ) = 0 dt Misalkan:
Sh − Sh0 = u ⇒ S˙h = u˙ Ih − I 0 = v ⇒ I˙h = v˙
h Nh − Nh0 Sv − Sv0 Iv − Iv0
= x ⇒ N˙h = x˙ = y ⇒ S˙v = y˙
(4.34)
= z ⇒ I˙v = z˙
Deret Taylor dari sistem (4.32) disekitar titik tetap adalah
(Sh0 , Ih0 , Nh0 , Sv0 , Iv0 )
dSh ∂A = A(Sh0 , Ih0 , Nh0 , Sv0 , Iv0 ) + (Sh − Sh0 ) + dt ∂Sh ∂A ∂A + (Nh − Nh0 ) + (Ih − Ih0 ) ∂Ih ∂Nh ∂A ∂A (Sv − Sv0 ) + (Iv − Iv0 ) + ··· ∂Sv ∂Iv dIh ∂B = B(Sh0 , Ih0 , Nh0 , Sv0 , Iv0 ) + (Sh − Sh0 ) + dt ∂Sh ∂B ∂B (Ih − Ih0 ) + (Nh − Nh0 ) + ∂Ih ∂Nh ∂B ∂B (Sv − Sv0 ) + (Iv − Iv0 ) + ··· ∂Sv ∂Iv dNh ∂C = C(Sh0 , Ih0 , Nh0 , Sv0 , Iv0 ) + (Sh − Sh0 ) + dt ∂Sh ∂C ∂C (Ih − Ih0 ) + (Nh − Nh0 ) + ∂Ih ∂Nh ∂C ∂C (Sv − Sv0 ) + (Iv − Iv0 ) + ··· ∂Sv ∂Iv
47 dSv ∂D = D(Sh0 , Ih0 , Nh0 , Sv0 , Iv0 ) + (Sh − Sh0 ) + dt ∂Sh ∂D ∂D + (Nh − Nh0 ) + (Ih − Ih0 ) ∂Ih ∂Nh ∂D ∂D (Sv − Sv0 ) + (Iv − Iv0 ) + ··· ∂Sv ∂Iv ∂E dIv = E(Sh0 , Ih0 , Nh0 , Sv0 , Iv0 ) + (Sh − Sh0 ) + dt ∂Sh ∂E ∂E + (Nh − Nh0 ) + (Ih − Ih0 ) ∂Ih ∂Nh ∂E ∂E (Sv − Sv0 ) + (Iv − Iv0 ) + ··· ∂Sv ∂Iv Berdasarkan persamaan (4.33), maka linearisasi dari sistem (4.32) adalah ∂A ∂A dSh = (Sh − Sh0 ) + (Ih − Ih0 ) + dt ∂Sh ∂Ih ∂A ∂A (Nh − Nh0 ) + (Sv − Sv0 ) + ∂Nh ∂Sv ∂A (Iv − Iv0 ) ∂Iv dIh ∂B ∂B = (Sh − Sh0 ) + (Ih − Ih0 ) + dt ∂Sh ∂Ih ∂B ∂B (Nh − Nh0 ) + (Sv − Sv0 ) + ∂Nh ∂Sv ∂B (Iv − Iv0 ) ∂Iv dNh ∂C ∂C = (Sh − Sh0 ) + (Ih − Ih0 ) + dt ∂Sh ∂Ih ∂C ∂C (Nh − Nh0 ) + (Sv − Sv0 ) + ∂Nh ∂Sv
48 ∂C ∂Iv ∂D ∂D = (Sh − Sh0 ) + (Ih − Ih0 ) + ∂Sh ∂Ih ∂D ∂D + (Sv − Sv0 ) + (Nh − Nh0 ) ∂Nh ∂Sv ∂D (Iv − Iv0 ) ∂Iv ∂E ∂E = (Sh − Sh0 ) + (Ih − Ih0 ) + ∂Sh ∂Ih ∂E ∂E (Nh − Nh0 ) + (Sv − Sv0 ) + ∂Nh ∂Sv ∂E (Iv − Iv0 ) ∂Iv (Iv − Iv0 )
dSv dt
dIv dt
Dengan menggunakan permisalan (4.34), maka hasil linearisasi dari sistem (4.32) seperti yang tertulis tersebut menjadi: dSh dt dIh dt dNh dt dSv dt dIv dt
∂A ∂Sh ∂B =u ∂Sh ∂C =u ∂Sh ∂D =u ∂Sh ∂E =u ∂Sh =u
∂A ∂Ih ∂B +v ∂Ih ∂C +v ∂Ih ∂D +v ∂Ih ∂E +v ∂Ih +v
∂A ∂Nh ∂B +x ∂Nh ∂C +x ∂Nh ∂D +x ∂Nh ∂E +x ∂Nh +x
∂A ∂Sv ∂B +y ∂Sv ∂C +y ∂Sv ∂D +y ∂Sv ∂E +y ∂Sv +y
∂A ∂Iv ∂B +z ∂Iv ∂C +z ∂Iv ∂D +z ∂Iv ∂E +z ∂Iv +z
(4.35)
Persamaan (4.35) dapat ditulis dalam bentuk matriks sebagai berikut:
49
dSh
dIh
dt
dt dNh dt dSv dt dIv
=
dt
∂A
∂A
∂A
∂A
∂A
∂Sh
∂Ih
∂Nh
∂Sv
∂Iv
∂B
∂B
∂B
∂B
∂B
∂Sh
∂Ih
∂Nh
∂Sv
∂Iv
∂C
∂C
∂C
∂C
∂C
∂Sh
∂Ih
∂Nh
∂Sv
∂Iv
∂D
∂D
∂D
∂D
∂D
∂Sh
∂Ih
∂Nh
∂Sv
∂Iv
∂E
∂E
∂E
∂E
∂E
∂Sh
∂Ih
∂Nh
∂Sv
∂Iv
"
u v x y z
#
Matriks Jacobian dari matriks tersebut adalah ∂A ∂A ∂A ∂A ∂A J =
∂Sh
∂Ih
∂Nh
∂Sv
∂Iv
∂B
∂B
∂B
∂B
∂B
∂Sh
∂Ih
∂Nh
∂Sv
∂Iv
∂C
∂C
∂C
∂C
∂C
∂Sh
∂Ih
∂Nh
∂Sv
∂Iv
∂D
∂D
∂D
∂D
∂D
∂Sh
∂Ih
∂Nh
∂Sv
∂Iv
∂E
∂E
∂E
∂E
∂E
∂Sh
∂Ih
∂Nh
∂Sv
∂Iv
Selanjutnya akan dicari matriks Jacobian dari sistem (4.32) dengan mendiferensialkannya sebagai berikut ∂ β1 Sh Ih β2 Sh Iv ∂A = b1 − − − µh Sh ∂Sh ∂Sh Nh Nh β1 Ih β2 Iv − − µh (4.36) = − Nh Nh ∂A ∂ β1 Sh Ih β2 Sh Iv = b1 − − − µh S h ∂Ih ∂Ih Nh Nh β1 Sh = − (4.37) Nh ∂A ∂ β1 Sh Ih β2 Sh Iv = b1 − − − µh Sh ∂Nh ∂Nh Nh Nh β1 Sh Ih β2 Sh Iv = + (4.38) (Nh )2 (Nh )2
50 ∂A ∂Sv
= =
∂A ∂Iv
= =
∂B ∂Sh
= + =
∂B ∂Ih
= + =
∂B ∂Nh
= + =
∂B ∂Sv
= + =
∂ β1 Sh Ih β2 Sh Iv − − µh S h b1 − ∂Sv Nh Nh 0 (4.39) ∂ β1 Sh β2 Sh Iv b1 − − − µh Sh ∂Iv Nh Nh β 2 Sh − (4.40) Nh ∂ β1 Sh Ih β2 Sh Iv + ∂Sh Nh Nh ∂ (−γh Ih − δh Ih − µh Ih ) ∂Sh β1 Ih β2 Iv + (4.41) Nh Nh ∂ β1 Sh Ih β2 Sh Iv + ∂Ih Nh Nh ∂ (−γh Ih − δh Ih − µh Ih ) ∂Ih β1 Sh − γh − δh − µh (4.42) Nh ∂ β1 Sh Ih β2 Sh Iv + ∂Nh Nh Nh ∂ (−γh Ih − δh Ih − µh Ih ) ∂Nh β1 Sh Ih β2 Sh Iv − − (4.43) (Nh )2 (Nh )2 ∂ β1 Sh Ih β2 Sh Iv + ∂Sv Nh Nh ∂ (−γh Ih − δh Ih − µh Ih ) ∂Sv 0 (4.44)
51 ∂B ∂Iv
= + =
∂C ∂Sh
= =
∂C ∂Ih
= =
∂C ∂Nh
= =
∂C ∂Sv
= =
∂C ∂Iv
= =
∂D ∂Sh
= =
∂D ∂Ih
= =
∂D ∂Nh
= =
∂ β1 Sh Ih β2 Sh Iv + ∂Iv Nh Nh ∂ (−γh Ih − δh Ih − µh Ih ) ∂Iv β2 Sh Nh ∂ (b1 − µh Nh − δh Ih ) ∂Sh 0 ∂ (b1 − µh Nh − δh Ih ) ∂Ih −δh ∂ (b1 − µh Nh − δh Ih ) ∂Nh −µh ∂ (b1 − µh Nh − δh Ih ) ∂Sv 0 ∂ (b1 − µh Nh − δh Ih ) ∂Iv 0 ∂ β3 Sv Ih b2 − − µv Sv ∂Sh Nh 0 ∂ β3 Sv Ih b2 − − µv Sv ∂Ih Nh β3 Sv − Nh ∂ β3 Sv Ih b2 − − µv Sv ∂Nh Nh β3 Sv Ih (Nh )2
(4.45)
(4.46)
(4.47)
(4.48)
(4.49)
(4.50)
(4.51)
(4.52)
(4.53)
52 ∂D ∂Sv
= =
∂D ∂Iv
= =
∂E ∂Sh
= =
∂E ∂Ih
= =
∂E ∂Nh
= =
∂E ∂Sv
= =
∂E ∂Iv
= =
∂ β3 Sv Ih − µv Sv b2 − ∂Sv Nh β3 Ih − µv − Nh ∂ β3 Sv Ih b2 − − µv Sv ∂Sv Nh 0 ∂ β3 Sv Ih − δ v Iv − µ v Iv ∂Sh Nh 0 ∂ β3 Sv Ih − δ v Iv − µ v Iv ∂Ih Nh β3 Sv Nh ∂ β3 Sv Ih − δv Iv − µv Iv ∂Nh Nh β3 Sv Ih − (Nh )2 ∂ β3 Sv Ih − δv Iv − µv Iv ∂Sv Nh β 3 Ih Nh ∂ β3 Sv Ih − δv Iv − µv Iv ∂Iv Nh −δv − µv
(4.54)
(4.55)
(4.56)
(4.57)
(4.58)
(4.59)
(4.60)
53 Dari hasil turunan persamaan (4.36) sampai (4.60), dapat ditulis dalam bentuk matriks Jacobian sebagai berikut: β1 Ih β2 Iv β1 Sh β1 Sh Ih β2 Sh Iv −
J =
− − µh Nh Nh β2 Iv β1 Ih + Nh Nh 0
−
β1 Sh Nh
− γh − δh − µh
0 0
−δh β3 Sv − Nh β3 Sv Nh
0 0
−
Nh
0 β3 Ih
− µv Nh β3 Ih Nh
−
β2 Sh
Nh β2 Sh Nh 0 0
−
(Nh )2 β1 Sh Ih
+
−
(Nh )2 β2 Sh Iv
(Nh )2 (Nh )2 −µh β3 Sv Ih − (Nh )2 β3 Sv Ih − (Nh )2
−δv − µv
4.6.1
Kestabilan Lokal Titik Setimbang Bebas Penyakit Telah diketahui sebelumnya bahwa titik setimbang bebas penyakit adalah Ef = (Sh0 , Ih0 , Nh0 , Sv0 , Iv0 ) = b1 b1 b2 , 0, , , 0 , maka µh µh µv −µh −β1 0 0 −β2 0 β1 − (γh + δh + µh ) 0 0 β2 0 −δh −µh 0 0 J(Ef ) = β3 b2 µh 0 − 0 −µv 0 b1 µv β3 b2 µh 0 0 0 −(δv + µv ) b1 µv Untuk mempermudah mencari persamaan karakteristiknya, maka pada J(Ef ) akan diubah ke dalam
54 bentuk matriks segitiga atas dengan cara OBE sebagai berikut. −µh −β1 0 0 −β2 0 β1 − (γh + δh + µh ) 0 0 β2 0 −δh −µh 0 0 J(Ef ) = β3 b2 µh 0 0 −µv 0 − b1 µv β3 b2 µh 0 0 0 −(δv + µv ) b1 µv ∼ B3 + ∼ B4 + ∼ B5 −
δh B2 (β1 − (γh + δh + µh )) β3 b2 µh b1 µv
(β1 − (γh + δh + µh )) β3 b2 µh b1 µv
(β1 − (γh + δh + µh )) " J(Ef ) =
0 0 0 −µv 0
−µh 0 0 0 0
B2 B2
−β1 β1 − (γh + δh + µh ) 0 0 0
−β2 β2 0 0 −(δv + µv )(β1 − (γh + δh + µh ))b1 µv − β2 β3 b2 µh b1 µv (β1 − (γh + δh + µh ))
0 0 −µh 0 0
Selanjutnya dicari persamaan karakteristik dari matriks Jacobian tersebut dengan menggunakan | J(Ef ) − λI |= 0
55 Sehingga "
−µh 0 0 0 0
−β1 β1 − (γh + δh + µh ) 0 0 0
0 0 −µh 0 0
−β2 β2 0 0 −(δv + µv )(β1 − (γh + δh + µh ))b1 µv − β2 β3 b2 µh b1 µv (β1 − (γh + δh + µh ))
Maka
0 0 0 −µv − λ 0
−µh − λ 0 0 0 0
−β1 β1 − (γh + δh + µh ) − λ 0 0 0
0 0 0 −µv 0
"
−λ
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
0 0 0 0 1
0 0 −µh − λ 0 0
−β2 β2 0 0 −(δv + µv )(β1 − (γh + δh + µh ))b1 µv − β2 β3 b2 µh b1 µv (β1 − (γh + δh + µh )) − λ
=0
Dari matriks Jacobian tersebut maka akan diperoleh persamaan karakteristik sebagai berikut. (−µh − λ) (β1 − (γh + δh + µh ) − λ) (−µh − λ)(−µv − λ) −(δv + µv )(β1 − (γh + δh + µh )b1 µv ) − β2 β3 b2 µh −λ =0 b1 µv (β1 − (γh + δh + µh )) Sehingga akan diperoleh nilai karakteristiknya sebagai berikut.
eigen
dari
λ1 = −µh < 0 λ2 = β1 − (γh + δh + µh ) = (γh + δh + µh )( untuk
β1 − 1) (γh + δh + µh )
β1 < 1 maka (γh + δh + µh ) λ2 = β1 − (γh + δh + µh ) < 0
akar
# =0
56 λ3 = −µh < 0 λ4 = −µv < 0 −(δv + µv ) (β1 − (γh + δh + µh )) b1 µv − β2 β3 b2 µh λ5 = b1 µv (β1 − (γh + δh + µh )) −(δv + µv ) (β1 − (γh + δh + µh )) b1 µv − β2 β3 b2 µh = b1 µv (β1 − (γh + δh + µh )) β2 β3 b2 µh + b1 µv (δv + µv ) (β1 − (γh + δh + µh )) = − b1 µv (β1 − (γh + δh + µh )) β2 β3 b2 µh = − − (δv + µv ) b1 µv (β1 − (γh + δh + µh )) β2 β3 b2 µh −1 = (δv + µv ) − b1 µv (δv + µv )(β1 − (γh + δh + µh )) = (δv + µv ) β2 β3 b2 µh (− b1 µv (δv + µv )(γh + δh + µh ) (γh +δβh1+µh ) − 1 β1 − 1 (γh +δh +µh ) ) − β1 − 1 (γh +δh +µh ) =
(δv + µv )
β1 (γh +δh +µh )
−
−1
β2 β3 b2 µh β1 − −1 b1 µv (δv + µv )(γh + δh + µh ) (γh + δh + µh ) (δv + µv ) = − β1 − 1 (γh +δh +µh ) β2 β3 b2 µh β1 + −1 b1 µv (δv + µv )(γh + δh + µh ) (γh + δh + µh ) (δv + µv ) = − β1 − 1 (γh +δh +µh ) β2 β3 b2 µh β1 + −1 b1 µv (δv + µv )(γh + δh + µh ) (γh + δh + µh )
57 untuk R0 = λ5 =
β2 β3 b2 µh b1 µv (δv +µv )(γh +δh +µh )
+
β1 (γh +δh +µh )
< 1 maka
−(δv + µv ) (β1 − (γh + δh + µh )) b1 µv − β2 β3 b2 µh <0 b1 µv (β1 − (γh + δh + µh ))
Karena nilai eigen (λ1 , λ2 , λ3 , λ4 , dan λ5 ) bernilai negatif pada bagian realnya maka berdasarkan akarakar karakteristik b1 b1 b2 (nilai eigen λ) maka titik setimbang Ef = , 0, , , 0 µh µh µv stabil lokal asimtotis dan titik kesetimbangan bebas penyakit bersifat stabil untuk R0 < 1. 4.6.2 Kestabilan Lokal Titik Setimbang Endemik Telah diketahui sebelumnya bahwa titik setimbang endemik E1 = (Sh∗ , Ih∗ , Nh∗ , Sv∗ , Iv∗ ), dalam hal ini Ih∗ selalu positif dengan b1 − (γh + δh + µh )Ih∗ µh b2 (b1 − δh Ih∗ ) Sv∗ = µh β3 Ih∗ + µv (b1 − δh Ih∗ ) µh b2 β3 Ih∗ Iv∗ = (µv + δv )(µh β3 Ih∗ + µv (b1 − δh Ih∗ )) b1 − δh Ih∗ Nh∗ = µh Sh∗ =
Pada titik setimbang E1 = (Sh∗ , Ih∗ , Nh∗ , Sv∗ , Iv∗ ) matrik Jacobiannya adalah β1 I ∗ β2 I ∗ ∗ β1 S ∗ β1 S ∗ I ∗ β 2 S ∗ Iv J(E1 ) =
− N ∗h − N ∗v − µh h h ∗ ∗ β 1 Ih β2 I v + N ∗ N∗ h
h
0 0 0
∗ β1 Sh N∗ h
− N ∗h h
− γh − δh − µh −δh ∗ β3 S v − N ∗ h ∗ β3 S v N∗ h
h h + h (N ∗ )2 (N ∗ )2 h∗ ∗ h∗ ∗ β1 S h Ih β2 Sh Iv − − (N ∗ )2 (N ∗ )2 h h −µh β S∗ I ∗ − 3 ∗v 2h (N ) h β S∗ I ∗ − 3 ∗v 2h (N ) h
58 ∗ ∗
β S Ih − 2Nh ∗
0
h ∗ β2 S h N∗ h
0 0 ∗ β 3 Ih − N ∗ − µv
0
0
h ∗ β 3 Ih N∗ h
−δv − µv
Untuk mempermudah mencari persamaan karakteristiknya, maka pada J(E1 ) akan diubah ke dalam bentuk matriks segitiga atas dengan cara OBE sebagai berikut: ∗
∗
∗
β1 I h β2 I v − N ∗ − N ∗ − µh h h ∗ β1 I h N∗ h
J(E1 ) =
β1 Sh − N ∗ ∗ β1 S h N∗ h
β I∗
2 v + N ∗ h 0
h ∗ β3 S v N∗ h
0
∗ ∗
β S Ih − 2Nh ∗
0
h ∗ β2 S h N∗ h
0 0
0
∗
β 3 Ih − N ∗ − µv
0
h ∗ β 3 Ih N∗ h
∼ B2 + −β
−β1 Ih∗ −β2 Iv∗ Nh∗
− γh − δh − µh −δh ∗ β3 S v − N ∗
0
h
∗ ∗ ∗ I∗ β1 S h h + β2 Sh Iv (N ∗ )2 (N ∗ )2 h∗ ∗ h∗ ∗ β S I β S I − 1 ∗h 2h − 2 ∗h 2v (N ) (N ) h h −µh ∗ I∗ β3 S v h − (N ∗ )2 h β S∗ I ∗ − 3 ∗v 2h (N ) h
−δv − µv
∗ ∗ ∗ 1 Ih −β2 Iv −µh Nh Nh∗
B1
sehingga J(E1 ) =
∗ −β I ∗ −µ N ∗ −β1 Ih 2 v h h N∗ h
0 0 0 0
∗ β1 S h − N ∗ h ∗ ∗ +β I ∗ +µ N ∗ ) −β1 µh Sh +(γh +δh +µh )(β1 Ih 2 v h h ∗ −µ N ∗ −β1 I ∗ −β2 Iv h h h −δh ∗ β3 Sv − N ∗ h ∗ β3 Sv N∗ h
59 ∗ I ∗ +β S ∗ I ∗ β1 Sh 2 h v h (N ∗ )2 h ∗ I ∗ +β ∗ ∗ β1 S h S 2 h h I v µh ∗ −µ N ∗ )(N ∗ ) (−β1 I ∗ −β2 Iv h h h h −µh ∗ I∗ β3 Sv h (N ∗ )2 h β S∗ I ∗ − 3 ∗v 2h (N ) h
∗
β2 Sh − N ∗
0
h ∗µ −β2 Sh h ∗ −µ N ∗ −β1 I ∗ −β2 Iv h h h
0 0
0
∗ −µ N ∗ −β3 Ih v h N∗ h ∗ β3 I h N∗ h
0
−(δv + µv )
Misalkan −β1 Ih∗ − β2 Iv∗ − µh Nh∗ Nh∗ −β1 µh Sh∗ + (γh + δh + µh )(β1 Ih∗ + β2 Iv∗ + µh Nh∗ ) B= −β1 Ih∗ − β2 Iv∗ − µh Nh∗ A=
Untuk mempermudah perhitungan selanjutnya, maka W = −β1 µh Sh∗ + (γh + δh + µh )(β1 Ih∗ + β2 Iv∗ + µh Nh∗ ) Sehingga B dapat ditulis kembali sebagai berikut B=
W − β2 Iv∗ − µh Nh∗
−β1 Ih∗
Sehingga matriknya dapat dtuliskan kembali sebagai berikut: ∗ β1 S ∗ β1 S ∗ I ∗ +β2 S ∗ Iv A
J(E1 ) =
− N ∗h h
0
B
0
−δh
0
β3 S v − N ∗
0
h h h (N ∗ )2 h ∗ I ∗ +β ∗ I∗ µ β1 Sh S 2 h v h h ∗ −µ N ∗ )(N ∗ ) (−β1 I ∗ −β2 Iv h h h h −µh ∗ I∗ β3 Sv h (N ∗ )2 h β S∗ I ∗ − 3 ∗v 2h (N ) h
∗
h ∗ β3 S v N∗ h
∗
0 0 0
∗ −µ N ∗ −β3 Ih v h N∗ h∗ β3 I h N∗ h
∼ B3 +
δh B B2
β2 Sh − N ∗
h ∗ −β2 Sh µh ∗ −µ N ∗ −β1 I ∗ −β2 Iv h h h
0
0 −(δv + µv )
60
∼ B4 +
B2
B
∼ B5 −
∗ β3 Sv N∗ h
∗ β3 Sv N∗ h
B2
B
Sehingga menghasilkan matriks sebagai berikut ∗ I ∗ +β S ∗ I ∗ ∗ β1 Sh β1 Sh 2 h v h A
J(E1 ) =
− N∗ h
0
B
0 0 0
0 0 0
(N ∗ )2 h ∗ I ∗ +β ∗ ∗ β1 Sh 2 Sh Iv µh h ∗ ∗ −µ N ∗ )(N ∗ ) (−β1 I −β2 Iv h h h h
C P Q
∗
0 0 0
∗ −µ N ∗ −β3 Ih v h N∗ h∗ β3 I h N∗ h
β2 Sh − N ∗
h ∗ −β2 Sh µh ∗ −µ N ∗ −β1 I ∗ −β2 Iv h h h
R S
T
dengan (δh (β1 Sh∗ Ih∗ + β2 Sh∗ Iv∗ µh )) W β3 Sv∗ Ih∗ W + β3 Sv∗ (β1 Sh∗ Ih∗ + β2 Sh∗ Iv∗ µh ) (Nh∗ )2 W −β3 Sv∗ Ih∗ W − β3 Sv∗ (β1 Sh∗ Ih∗ + β2 Sh∗ Iv∗ µh ) (Nh∗ )2 W −δh β2 Sh∗ µh W −β3 β2 Sh∗ Sv∗ µh Nh∗ W −Nh∗ (δv + µv )W + β2 β3 Sv∗ Sh∗ µh Nh∗ W
C = −µh + P = Q= R= S= T =
Selanjutnya kembali dilakukan OBE sebagai berikut ∼ B4 + −P C B3
61 ∼ B5 + −Q C B3 Sehingga didapatkan matriks sebagai berikut ∗ ∗ I ∗ +β S ∗ I ∗ β1 Sh β1 Sh 2 h v h A
J(E1 ) =
− N∗ h
0
B
0 0 0
0 0 0
(N ∗ )2 h ∗ I ∗ +β ∗ ∗ β1 Sh 2 Sh Iv µh h ∗ −µ N ∗ )(N ∗ ) (−β1 I ∗ −β2 Iv h h h h
C 0 0
∗
β2 Sh − N ∗
0
h ∗ −β2 Sh µh ∗ −µ N ∗ −β1 I ∗ −β2 Iv h h h
0 0
R
∗ −µ N ∗ −β3 Ih v h N∗ h∗ β3 I h N∗ h
U
V
dengan −Nh∗ β3 β2 Sh∗ Sv∗ µh (−µh W + δh (β1 Sh∗ Ih∗ + β2 Sh∗ Iv∗ µh )) W (Nh∗ )2 −µh W + δh (β1 Sh∗ Ih∗ + β2 Sh∗ Iv∗ µh ) δh β2 Sh∗ µh (β3 Sv∗ Ih∗ W + β3 Sv∗ (β1 Sh∗ Ih∗ + β2 Sh∗ Iv∗ µh )) + W (Nh∗ )2 −µh W + δh (β1 Sh∗ Ih∗ + β2 Sh∗ Iv∗ µh ) N ∗ (−µh W + δh (β1 Sh∗ Ih∗ + β2 Sh∗ Iv∗ µh )) (−(δv + µv )Nh∗ W ) V = h W (Nh∗ )2 −µh W + δh (β1 Sh∗ Ih∗ + β2 Sh∗ Iv∗ µh ) N ∗ (−µh W + δh (β1 Sh∗ Ih∗ + β2 Sh∗ Iv∗ µh )) (β3 β2 Sh∗ Sv∗ µh ) + h W (Nh∗ )2 −µh W + δh (β1 Sh∗ Ih∗ + β2 Sh∗ Iv∗ µh ) δh β2 Sh∗ µh (−β3 Sv∗ Ih∗ W − β3 Sv∗ (β1 Sh∗ Ih∗ + β2 Sh∗ Iv∗ µh )) + W (Nh∗ )2 −µh W + δh (β1 Sh∗ Ih∗ + β2 Sh∗ Iv∗ µh ) U=
Selanjutnya kembali dilakukan OBE sebagai berikut −
∼ B5 +
∗ β3 Ih N∗ h
C
B4
Sehingga didapatkan matriks sebagai berikut ∗ ∗ I ∗ +β S ∗ I ∗ β1 Sh β1 Sh 2 h v h A
J(E1 ) =
− N∗ h
0
B
0 0 0
0 0 0
(N ∗ )2 h ∗ I ∗ +β ∗ ∗ β1 Sh 2 Sh Iv µh h ∗ −µ N ∗ )(N ∗ ) (−β1 I ∗ −β2 Iv h h h h
C 0 0
62 ∗
0
β2 Sh − N ∗
R
h ∗ −β2 Sh µh ∗ −µ N ∗ −β1 I ∗ −β2 Iv h h h
0 0
∗ −µ N ∗ −β3 Ih v h N∗ h
U
0
E
dengan Nh∗ L (−(δv + µv )Nh∗ W ) (−β3 Ih∗ − µv Nh∗ ) W (Nh∗ )2 L(−β3 Ih∗ − µv Nh∗ ) N ∗ L(β3 β2 Sh∗ Sv∗ µh )(−β3 Ih∗ − µv Nh∗ ) + h W (Nh∗ )2 L(−β3 Ih∗ − µv Nh∗ ) δh β2 Sh∗ µh L(−β3 Ih∗ − µv Nh∗ ) + W (Nh∗ )2 L(−β3 Ih∗ − µv Nh∗ ) β3 Ih∗ Nh∗ Lβ3 β2 Sh∗ Sv∗ µh + W (Nh∗ )2 L(−β3 Ih∗ − µv Nh∗ )
E=
dalam persamaan E tersebut dimisalkan kembali: L = (−µh W + δh (β1 Sh∗ Ih∗ + β2 Sh∗ Iv∗ µh )) Sehingga didapatkan bentuk matriks segitiga atas sebagai berikut:
∗
A
β1 Sh − N ∗ h
J(E1 ) =
0
B
0 0 0
0 0 0 ∗
0 0 0 D 0
β2 Sh − N ∗
h ∗ −β2 Sh µh ∗ −µ N ∗ −β1 I ∗ −β2 Iv h h h
R U E
∗ I ∗ +β S ∗ I ∗ β1 Sh 2 h v h (N ∗ )2 h ∗ I ∗ +β ∗ ∗ β1 Sh S 2 h h I v µh ∗ −µ N ∗ )(N ∗ ) (−β1 I ∗ −β2 Iv h h h h
C 0 0
dengan pemisalan sebagai berikut
63 −β1 Ih∗ − β2 Iv∗ − µh Nh∗ Nh∗ W B= −β1 Ih∗ − β2 Iv∗ − µh Nh∗ (δh (β1 Sh∗ Ih∗ + β2 Sh∗ Iv∗ µh )) C = −µh + W −β3 Ih∗ − µv Nh∗ D= Nh∗ (β3 β2 Sh∗ Sv∗ µh ) δh β2 Sh∗ µh + E = −(δv + µv ) + W Nh∗ W (Nh∗ )2 β32 Ih∗ β2 Sh∗ Sv∗ µh − W Nh∗ (β3 Ih∗ + µv Nh∗ ) A=
Selanjutnya dicari persamaan karakteristik dari matriks Jacobian tersebut dengan menggunakan | J(E1 ) − λI |= 0 Sehingga
∗
A
β1 S h − N ∗ h
0
B
0 0 0
0 0 0
C 0 0
∗
β2 Sh − N ∗
0 0 0 D 0
∗ I ∗ +β S ∗ I ∗ β1 S h 2 h v h (N ∗ )2 h ∗ I ∗ +β ∗ ∗ β1 S h S 2 h I v µh h ∗ −µ N ∗ )(N ∗ ) (−β1 I ∗ −β2 Iv h h h h
h ∗ −β2 Sh µh ∗ −µ N ∗ −β1 I ∗ −β2 Iv h h h
R U E
"
−λ
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
0 0 0 0 1
Maka
∗
A−λ
β1 Sh − N ∗ h
0
B−λ
0 0 0
0 0 0
∗ I ∗ +β S ∗ I ∗ β1 Sh 2 h v h (N ∗ )2 h ∗ I ∗ +β ∗ ∗ S β1 S h 2 h I v µh h ∗ −µ N ∗ )(N ∗ ) (−β1 I ∗ −β2 Iv h h h h
C−λ 0 0
# =0
64 ∗
β2 Sh − N ∗
0 0 0 D−λ 0
h ∗ −β2 Sh µh ∗ −µ N ∗ −β1 I ∗ −β2 Iv h h h
R U E−λ
=0
Dari matriks Jacobian tersebut maka akan diperoleh persamaan karakteristik sebagai berikut. (A − λ)(B − λ)(C − λ)(D − λ)(E − λ) = 0 Sehingga akan diperoleh nilai karakteristiknya sebagai berikut.
eigen
dari
akar
Untuk nilai A −β1 Ih∗ − β2 Iv∗ − µh Nh∗ Nh∗ (β1 Ih∗ + β2 Iv∗ + µh Nh∗ ) =− Nh∗
A=
Maka λ1 = A < 0 Untuk nilai B B=−
W (β1 Ih∗ + β2 Iv∗ + µh Nh∗ )
dengan W = −β1 µh Sh∗ + (γh + δh + µh )(β1 Ih∗ + β2 Iv∗ + µh Nh∗ ) dimana ((γh + δh + µh )(β1 Ih∗ + β2 Iv∗ + µh Nh∗ )) > −β1 µh Sh∗ Sehingga W > 0, karena (β1 Ih∗ + β2 Iv∗ + µh Nh∗ ) > 0, maka diperoleh λ2 = B < 0
65
Untuk nilai C (δh (β1 Sh∗ Ih∗ + β2 Sh∗ Iv∗ µh )) C = −µh + W (µh W − δh β1 Sh∗ Ih∗ − δh β2 Sh∗ Iv∗ µh ) =− W dengan (µh W − δh β1 Sh∗ Ih∗ − δh β2 Sh∗ Iv∗ µh ) > 0 Karena W > 0, maka λ3 = C < 0 Untuk nilai D −β3 Ih∗ − µv Nh∗ Nh∗ (β3 Ih∗ + µv Nh∗ ) =− Nh∗
D=
Maka λ4 = D < 0 Untuk nilai E (β3 β2 Sh∗ Sv∗ µh ) δh β2 Sh∗ µh E = −(δv + µv ) + + W Nh∗ W (Nh∗ )2 β32 Ih∗ β2 Sh∗ Sv∗ µh − W Nh∗ (β3 Ih∗ + µv Nh∗ ) ! (δv + µv )(Nh∗ )2 W (β3 Ih∗ + µv Nh∗ ) − (Nh∗ (β3 β2 Sh∗ Sv∗ µh )) E=− W (Nh∗ )2 (β3 Ih∗ + µv Nh∗ ) ∗ 2 ∗ Nh β3 Ih β2 Sh∗ Sv∗ µh − ((δh β2 Sh∗ µh )(β3 Ih∗ + µv Nh∗ )) − W (Nh∗ )2 (β3 Ih∗ + µv Nh∗ )
66 dimana (δv + µv )(Nh∗ )2 W (β3 Ih∗ + µv Nh∗ ) + Nh∗ β32 Ih∗ β2 Sh∗ Sv∗ µh > (Nh∗ (β3 β2 Sh∗ Sv∗ µh )) + ((δh β2 Sh∗ µh )(β3 Ih∗ + µv Nh∗ )) dan juga W (Nh∗ )2 (β3 Ih∗ + µv Nh∗ ) > 0 Maka λ5 = E < 0 Karena nilai eigen (λ1 , λ2 , λ3 , λ4 , dan λ5 ) bernilai negatif pada bagian realnya maka berdasarkan akar akar karakteristik (nilai eigen λ) maka titik setimbang E1 = (Sh∗ , Ih∗ , Nh∗ , Sv∗ , Iv∗ ) stabil lokal asimtotis. 4.7
Analisa Bifurkasi Dalam hal ini, menggunakan titik kesetimbangan endemik g(Ih ) untuk mencari persamaan R0 yang optimum untuk membuat kurva bifurkasinya sehingga untuk R0 yang lebih kecil dari nilai optimum tidak terjadi penyebaran penyakit menular. Diketahui g(Ih ) sebagai berikut g(Ih ) = AIh∗2 + BIh∗ + C = 0 dengan A = K1 K2 (µh β3 (β1 − δh ) + µv δh (δh − β1 )), B = K1 K2 (µv b1 β1 − 2µv b1 δh + µh b1 β3 ) + K2 b1 β1 (µv δh − µh β3 ) + K1 µh b2 β2 β3 , C = µv b21 K1 K2 (1 − R0 ), µh b2 β2 β3 β1 R0 = + µv b1 K1 K2 K1 Dalam hal ini koefisien A selalu bernilai positif. Untuk koefisien C bergantung pada nilai R0 . Jika R0 < 1, maka
67 C > 0,sedangkan jika R0 > 1, maka C < 0. Jika R0 = 1, maka C = 0. Sehingga diperoleh g(Ih ) = AIh∗2 + BIh∗ + C = AIh∗2 + BIh∗ + C ∗ (1 − R0 ), dengan C ∗ = µv b21 K1 K2 = AIh∗2 + BIh∗ + C ∗ (1 − 1) = AIh∗2 + BIh∗ Karena g(Ih ) = 0, maka 0 = AIh∗2 + BIh∗ 0 = Ih∗ (AIh∗ + B) Karena Ih∗ 6= 0, maka AIh∗ + B = 0, sehingga AIh∗ + B = 0 AIh∗ = −B −B Ih∗ = A Dalam hal ini Ih∗ bernilai positif jika dan hanya jika A > 0 dan B < 0 atau A < 0 dan B > 0. Karena A selalu positif, maka hanya berlaku untuk A > 0 dan B < 0 Pada persamaan g(Ih ) terjadi bifurkasi mundur pada R0 = 1 jika hanya jika A > 0 dan B < 0 sehingga B 2 − 4AC > 0. Karena R0 = 1 maka C = 0 sehingga B 2 > 0. Pada saat terjadi bifurkasi mundur ada dua titik kesetimbangan endemik dan g(Ih∗ ) = 0 sehingga diperoleh √ −B − B 2 − 4AC ∗ Ih1 = √2A −B + B 2 − 4AC Ih∗2 = 2A
68 Dengan Ih∗1 dan Ih∗2 masing-masing bersifat stabil dan tidak stabil. Dengan demikian, persamaan g(Ih ) mempunyai dua penyelesaian positif, yang berhubungan dengan dua keseimbangan endemik jika hanya jika C > 0 atau R0 < 1 dan B < 0, A > 0, dan B 2 > 4AC. Dengan adanya bifurkasi mundur berakibat bahwa penyakit menular semakin sulit diberantas. Jika C > 0 dan B > 0 atau B 2 < 4AC, maka persamaan g(Ih ) menghasilkan penyelesaian tak real dan tidak ada keseimbangan endemiknya. Selanjutnya dari nilai R0 dan Ih akan disimulasikan dengan menggunakan MATLAB R2009a yang menghasilkan kurva bifurkasi dengan sumbu (x,y) yang merupakan (R0 ,Ih ). Di bawah ini merupakan kurva bifurkasi hasil dari simulasi. Adapun nilai parameter yang digunakan seperti yang terdapat pada tabel berikut ini : Tabel 4.1: Input Parameter Parameter Nilai Parameter b1 3 µh b2 12 µv β1 4 γh β2 16 δh β3 1 δv
Nilai 1 4 2 78 2
69
0.1 0.09 0.08
Populasi (Ih)
0.07 0.06
I1 I2
0.05 0.04 0.03 0.02 0.01 0
0
1 2 3 Bilangan Reproduksi Dasar (R0)
4
5
Gambar 4.2: Kurva Bifurkasi Mundur
Pada Gambar 4.2 menunjukkan bahwa telah terjadi bifurkasi mundur untuk R0 < 1 sehingga diperoleh tiga titik tetap, yang terdiri dari titik setimbang bebas penyakit, titik setimbang endemik stabil, dan titik setimbang endemik tidak stabil. Pada saat R0 = 1 terjadi bifurkasi Transkritikal, sedangkan pada titik 0.0834 terjadi bifurkasi Sadle Node. Pada saat R0 < 0.0834 tidak terjadi penyebaran penyakit (bebas penyakit), sedangkan 0.0834 < R0 < 1 ada dua titik kesetimbangan endemik, satu titik bersifat stabil (warna biru) dan yang lain tidak stabil (warna merah) dan juga satu titik kesetimbangan bebas penyakit bersifat stabil dan untuk R0 > 1 titik endemik stabil sangat besar sehingga penularan (endemik) sulit diatasi dan terdapat titik setimbang bebas
70 penyakit yang tidak stabil. Pada Gambar 4.2 menunjukkan kurva bifurkasi mundur dengan R0 = 1 atau R0 < 1 dengan nilai A > 0, B < 0, C > 0 untuk R0 < 1 dan A > 0, B < 0, C = 0 untuk R0 = 1. Selanjutnya dengan menggunakan nilai parameter yang digunakan seperti yang terdapat pada tabel berikut ini : Tabel 4.2: Input Parameter Parameter Nilai Parameter b1 1.7 µh b2 2 µv β1 1 γh β2 4 δh β3 0.48 δv
Nilai 0.6 0.8 0.4 1.4 0.8
Gambar 4.3 menunjukkan kurva bifurkasi maju dengan R0 > 1. Terdapat satu titik setimbang endemik dan titik setimbang bebas penyakit yang tidak stabil. Pada saat R0 = 1 terjadi bifurkasi Transkritikal. Pada saat R0 < 1 tidak terjadi penyebaran penyakit (bebas penyakit), sedangkan R0 > 1 terdapat titik endemik stabil. Pada Gambar 4.3 menunjukkan bahwa telah terjadi bifurkasi maju dengan R0 > 1 dengan mengakibatkan nilai A > 0, B > 0, C < 0. Dalam hal ini perubahan bifurkasi terjadi karena perubahan nilai R0 yang mengakibatkan nilai A, B, dan C berubah sehingga nilai titik puncaknya pun berubah.
71
3.5
3
Populasi (Ih)
2.5
2
I1 I2
1.5
1
0.5
0
0
1
2 3 4 Bilangan Reproduksi Dasar (R0)
5
Gambar 4.3: Kurva Bifurkasi Maju
4.8
Model Kendali Optimal Dalam model pada persamaan (4.1) sampai (4.5) diperluas dengan memasukkan jumlah populasi yang bergantung pada angka kematian dalam populasi vektor dan host (inang), yang didefinisikan sebagai berikut [1] µh = µ1 + µ2 Nh µv = µ3 + µ4 Nv dengan µ1 ≥ 0 dan µ3 ≥ 0 yang merupakan kematian sebelum terinfeksi dan µ2 ≥ 0 dan µ4 ≥ 0 adalah kematian setelah terinfeksi. Laju penambahan populasi individu dan vektor yang rentan terhadap penyakit (susceptible) diperbaharui
72 berdasarkan berikut [1] b1 → b1 + αh Nh b2 → b2 Nv dengan αh ≥ 0 adalah konstanta dan pengendaliannya sebagai berikut u1 : Pengendalian pada populasi individu yang rentan (susceptible(Sh )) yang terinfeksi oleh populasi individu (infected (Ih )). u2 : Pengendalian pada populasi yang rentan (susceptible(Sh )) yang terinfeksi oleh populasi vektor (infected (Iv )) dan pada populasi vektor yang rentan (susceptible (Sv )) yang terinfeksi oleh populasi individu (infected (Ih )). u3 : Pengendalian untuk populasi vektor. Sehingga diperoleh model interaksi dinamis sebagai berikut: 1. Besarnya laju populasi individu yang rentan terhadap penyakit (susceptible) dipengaruhi oleh angka kelahiran populasi individu yang rentan,penambahan jumlah populasi individu, pengendalian yang dilakukan pada populasi individu yang rentan yang telah terinfeksi oleh populasi individu yang terinfeksi dan populasi vektor yang terinfeksi sedangkan populasi akan menurun dengan adanya beberapa kejadian penularan penyakit, seperti populasi individu yang terinfeksi penyakit karena melakukan kontak langsung dengan individu yang terinfeksi, populasi individu yg terinfeksi penyakit karena melakukan kontak dengan vektor yang terinfeksi, Selain itu berkurang karena kematian alami dari populasi individu yang rentan dan kematian karena
73 terinfeksi pada jumlah populasi individu yang rentan. dSh = b1 + αh Nh − β1 Sh Ih (1 − u1 ) dt − β2 Sh Iv (1 − u2 ) − (µ1 + µ2 Nh )Sh 2. Besarnya laju populasi individu yang terinfeksi penyakit (infected) akan bertambah saat terdapat populasi individu yang terinfeksi penyakit akibat kontak langsung dengan individu yang terinfeksi maupun dengan vektor yang terinfeksi dan populasi akan menurun karena adanya pengendalian yang dilakukan pada populasi individu yang terinfeksi maupun vektor yang terinfeksi, adanya kejadian individu yang terinfeksi namun telah sembuh dan juga menurun karena kematian yang disebabkan karena terinfeksi penyakit maupun secara alami pada populasi individu yang terinfeksi dan juga pada jumlah populasi individu yang terinfeksi. dIh = β1 Sh Ih (1 − u1 ) + β2 Sh Iv (1 − u2 ) dt − γh Ih − δh Ih − (µ1 + µ2 Nh )Ih 3. Besarnya laju populasi individu yang sembuh dari penyakit (recovered) bergantung pada jumlah individu yang terinfeksi namun telah sembuh dan populasi akan berkurang saat terdapat kejadian kematian alami populasi individu yang sembuh dan kematian karena penyakit pada jumlah populasi individu yang sembuh. dRh = γh Ih − (µ1 + µ2 Nh )Rh dt 4. Besarnya laju populasi vektor yang rentan terhadap penyakit (susceptible) bergantung pada angka kelahiran
74 populasi vektor yang rentan, pengendalian yang dilakukan pada populasi vektor rentan yang telah terinfeksi oleh populasi individu yang terinfeksi dan populasi akan berkurang saat terjadi pengendalian pada angka kelahiran populasi vektor, saat angka populasi individu yang terinfeksi yang menjadikan vektor rentan terinfeksi, kematian alami populasi vektor, kematian karena penyakit depada jumlah populasi vektor, dan juga kematian populasi vektor yang dikenakan pengendalian. dSv = b2 Nv (1 − u3 ) − β3 Sv Ih (1 − u2 ) dt − (µ3 + µ4 Nv )Sv − r0 u3 Sv 5. Besarnya laju populasi vektor yang terinfeksi penyakit (infected) bergantung pada angka populasi individu yang terinfeksi yang menjadikan vektor rentan terinfeksi dan populasi akan berkurang akibat pengendalian yang dilakukan pada populasi vektor yang telah terinfeksi oleh populasi individu yang terinfeksi, angka kematian dari vektor yang terinfeksi yang disebabkan penyakit, angka kematian alami populasi vektor yang terinfeksi, kematian karena penyakit pada jumlah populasi vektor yang terinfeksi dan juga kematian populasi vektor yang terinfeksi yang dikenakan pengendalian. dIv = β3 Sv Ih (1 − u2 ) − δv Iv dt − (µ3 + µ4 Nv )Iv − r0 u3 Iv Dari penjelasan diatas, maka dapat dirumuskan model kendali optimal untuk penyakit vektor borne adalah sebagai
75 berikut [1]. dSh dt
= b1 + αh Nh − β1 Sh Ih (1 − u1 ) −
dIh dt
= β1 Sh Ih (1 − u1 ) + β2 Sh Iv (1 − u2 ) −
dRh dt dSv dt
β2 Sh Iv (1 − u2 ) − (µ1 + µ2 Nh )Sh
γh Ih − δh Ih − (µ1 + µ2 Nh )Ih
= γh Ih − (µ1 + µ2 Nh )Rh
(4.61)
= b2 Nv (1 − u3 ) − β3 Sv Ih (1 − u2 ) − (µ3 + µ4 Nv )Sv − r0 u3 Sv
dIv dt
= β3 Sv Ih (1 − u2 ) − δv Iv − (µ3 + µ4 Nv )Iv − r0 u3 Iv
dengan kondisi awal: Sh (0) > 0, Ih (0) > 0, Rh (0) > 0, Sv (0) > 0, Iv (0) > 0 Fungsi tujuan untuk sistem keadaan diberikan oleh Z T 1 2 2 2 J(u1 , u2 , u3 ) = A1 Ih + A2 Nv + (B1 u1 + B2 u2 + B3 u3 ) dt 2 0 Fungsi tujuan tersebut harus memenuhi sistem keadaan yang diberikan oleh A1 : bobot populasi inang (host) yang terinfeksi. A2 : bobot pada total populasi vektor. B1 : bobot untuk pemeriksaan donor darah. B2 : bobot untuk perlindungan personal (pengurangan vektor dan kontak manusia). B3 : pengendalian vektor. Fungsi tujuan dalam masalah kendali optimal ini
76 adalah untuk meminimalkan jumlah populasi inang (host) yang terinfeksi, jumlah vektor, dan meminimalkan biaya pengendalian u1 ,u2 , dan u3 . Biaya tersebut diperuntukan pada kendali yang diantaranya adalah [1]. 1. yang berasal dari sistem penyeleksian donor. 2. yang berasal dari biaya vaksinasi, obat nyamuk, dan keperluan lain yang berkaitan dengan upaya pengendalian. 3. yang berasal dari membasmi nyamuk.
penggunaan
pestisida
untuk
Tujuan yang akan dicari adalah untuk menemukan fungsi kendali sehingga J(u∗1 , u∗2 , u∗3 ) =
min (u1 ,u2 ,u3 )∈U
J(u1 , u2 , u3 )
harus memenuhi sistem pada model sistem persamaan (4.61) dengan U = ((u1 , u2 , u3 ) | ui (t) ∈ [0, 1]) 4.9
Penyelesaian Kendali Optimal Dalam penyelesaian kendali optimal, harus menentukan Lagrangian dan fungsi Hamiltonian. Dengan memperhatikan model kendali optimal, Lagrangian dari masalah kendali diberikan oleh 1 L = A1 Ih + A2 Nv + (B1 u21 + B2 u22 + B3 u23 ) 2
Selanjutnya akan dicari nilai minimal Lagrangian tersebut sehingga dapat didefinisikan hamiltonian H untuk masalah
77 kendali sebagai berikut H = L(Ih , Nv , u1 , u2 , u3 ) + λ1
dIh dRh dSh + λ2 + λ3 dt dt dt
dIv dSv + λ5 dt dt 1 = A1 Ih + A2 Nv + (B1 u21 + B2 u22 + B3 u23 ) + λ1 (b1 2 + αh Nh − β1 Sh Ih (1 − u1 ) − β2 Sh Iv (1 − u2 ) + λ4
− (µ1 + µ2 Nh )Sh ) + λ2 (β1 Sh Ih (1 − u1 ) + β2 Sh Iv (1 − u2 ) − γh Ih − δh Ih − (µ1 + µ2 Nh )Ih ) + λ3 (γh Ih − (µ1 + µ2 Nh )Rh ) + λ4 (b2 Nv (1 − u3 ) − β3 Sv Ih (1 − u2 ) − (µ3 + µ4 Nv )Sv − r0 u3 Sv ) + λ5 (β3 Sv Ih (1 − u2 ) − δv Iv − (µ3 + µ4 Nv )Iv − r0 u3 Iv ) 1 = A1 Ih + A2 Nv + (B1 u21 + B2 u22 + B3 u23 ) + λ1 b1 2 + λ1 αh Nh − λ1 β1 Sh Ih + λ1 β1 Sh Ih u1 − λ1 β2 Sh Iv + λ1 β2 Sh Iv u2 − λ1 µ1 Sh − λ1 µ2 Nh Sh + λ2 β1 Sh Ih − λ2 β1 Sh Ih u1 + λ2 β2 Sh Iv − λ2 β2 Sh Iv u2 − λ2 γh Ih − λ2 δh Ih − λ2 µ1 Ih − λ2 µ2 Nh Ih + λ3 γh Ih − λ3 µ1 Rh − λ3 µ2 Nh Rh + λ4 b2 Nv − λ4 b2 Nv u3 − λ4 β3 Sv Ih + λ4 β3 Sv Ih u2 − λ4 µ3 Sv − λ4 µ4 Nv Sv − λ4 r0 u3 Sv + λ5 β3 Sv Ih − λ5 β3 Sv Ih u2 − λ5 δv Iv − λ5 µ3 Iv − λ5 µ4 Nv Iv − λ5 r0 u3 Iv 1 1 1 = A1 Ih + A2 (Sv + Iv ) + B1 u21 + B2 u22 + B3 u23 2 2 2 + λ1 b1 + λ1 αh (Sh + Rh + Ih ) − λ1 β1 Sh Ih + λ1 β1 Sh Ih u1 − λ1 β2 Sh Iv + λ1 β2 Sh Iv u2 − λ1 µ1 Sh − λ1 µ2 (Sh + Rh + Ih )Sh + λ2 β1 Sh Ih − λ2 β1 Sh Ih u1 + λ2 β2 Sh Iv − λ2 β2 Sh Iv u2 − λ2 γh Ih − λ2 δh Ih − λ2 µ1 Ih − λ2 µ2 (Sh + Rh + Ih )Ih
78
+λ3 γh Ih − λ3 µ1 Rh − λ3 µ2 (Sh + Rh + Ih )Rh + λ4 b2 (Sv + Iv ) −λ4 b2 (Sv + Iv )u3 − λ4 β3 Sv Ih + λ4 β3 Sv Ih u2 − λ4 µ3 Sv −λ4 µ4 (Sv + Iv )Sv − λ4 r0 u3 Sv + λ5 β3 Sv Ih − λ5 β3 Sv Ih u2 H
=
−λ5 δv Iv − λ5 µ3 Iv − λ5 µ4 (Sv + Iv )Iv − λ5 r0 u3 Iv 1 1 1 A1 Ih + A2 Sv + A2 Iv + B1 u21 + B2 u22 + B3 u23 + λ1 b1 2 2 2 +λ1 αh Sh + λ1 αh Rh + λ1 αh Ih − λ1 β1 Sh Ih +λ1 β1 Sh Ih u1 − λ1 β2 Sh Iv + λ1 β2 Sh Iv u2 − λ1 µ1 Sh −λ1 µ2 Sh2 − λ1 µ2 Rh Sh − λ1 µ2 Ih Sh + λ2 β1 Sh Ih −λ2 β1 Sh Ih u1 + λ2 β2 Sh Iv − λ2 β2 Sh Iv u2 − λ2 γh Ih −λ2 δh Ih − λ2 µ1 Ih − λ2 µ2 Sh Ih − λ2 µ2 Rh Ih − λ2 µ2 Ih2 +λ3 γh Ih − λ3 µ1 Rh − λ3 µ2 Sh Rh − λ3 µ2 Rh2 + −λ3 µ2 Ih Rh +λ4 b2 Sv + λ4 b2 Iv − λ4 b2 Sv u3 − λ4 b2 Iv u3 − λ4 β3 Sv Ih +λ4 β3 Sv Ih u2 − λ4 µ3 Sv − λ4 µ4 Sv2 − λ4 µ4 Iv Sv −λ4 r0 u3 Sv + λ5 β3 Sv Ih − λ5 β3 Sv Ih u2 − λ5 δv Iv −λ5 µ3 Iv − λ5 µ4 Sv Iv − λ5 µ4 Iv2 − λ5 r0 u3 Iv
Diberikan Sh∗ , Ih∗ , Rh∗ , Sv∗ , dan Iv∗ adalah persamaan state yang optimal dengan variabel kendali optimal (u∗1 , u∗2 , u∗3 ) untuk kendali optimal pada sistem persamaan (4.61). Kemudian terdapat variabel adjoint λi , untuk i = 1, 2, 3, 4, 5. Dengan menggunakan Hamiltonian pada (4.62) maka variabel adjoint dapat dicari dengan cara berikut ∂H 0 λ1 = − ∂Sh = −λ1 αh + λ1 β1 Ih − λ1 β1 Ih u1 + λ1 β2 Iv − λ1 β2 Iv u2 − λ2 β1 Ih + λ2 β1 Ih u1 − λ2 β2 Iv + λ2 β2 Iv u2 + λ1 µ1 + 2λ1 µ2 Sh + λ1 µ2 Ih + λ1 µ2 Rh + λ2 µ2 Ih + λ3 µ2 Rh = −λ1 αh + (λ1 − λ2 )(β1 (1 − u1 )Ih + β2 (1 − u2 )Iv )
(4.67)
79 + λ1 µ1 + λ1 µ2 Sh + λ1 µ2 Sh + λ1 µ2 Ih + λ1 µ2 Rh + λ2 µ2 Ih + λ3 µ2 Rh 0
λ1 = −λ1 αh + (λ1 − λ2 )(β1 (1 − u1 )Ih + β2 (1 − u2 )Iv ) + λ1 µ1 + λ1 µ2 Sh + λ1 µ2 Ih + λ1 µ2 Rh + λ1 µ2 Sh + λ2 µ2 Ih + λ3 µ2 Rh = −λ1 αh + (λ1 − λ2 )(β1 (1 − u1 )Ih + β2 (1 − u2 )Iv ) + λ1 µ1 + λ1 µ2 Nh + λ1 µ2 Sh + λ2 µ2 Ih + λ3 µ2 Rh = −λ1 αh + (λ1 − λ2 )(β1 (1 − u1 )Ih + β2 (1 − u2 )Iv ) + (µ1 + µ2 Nh )λ1 + λ1 µ2 Sh + λ2 µ2 Ih + λ3 µ2 Rh ∂H λ2 = − ∂Ih = −A1 − λ1 αh + λ1 β1 Sh − λ1 β1 Sh u1 + µ2 Sh λ1 0
− λ2 β1 Sh + λ2 β1 Sh u1 + λ2 γh + λ2 δh + λ2 µ1 + λ2 µ2 Sh + λ2 µ2 Rh + 2λ2 µ2 Ih − λ3 γh + λ3 µ2 Rh + λ4 β3 Sv − λ4 β3 Sv u2 − λ5 β3 Sv + λ5 β3 Sv u2 = −αh λ1 − A1 + λ1 β1 Sh − λ1 β1 Sh u1 − λ2 β1 Sh + λ2 β1 Sh u1 + µ2 Sh λ1 + λ2 γh + λ2 δh + λ2 µ1 + λ2 µ2 Sh + λ2 µ2 Rh + µ2 Ih λ2 + µ2 Ih λ2 − λ3 γh + λ3 µ2 Rh + λ4 β3 Sv − λ4 β3 Sv u2 − λ5 β3 Sv + λ5 β3 Sv u2 = −αh λ1 − A1 + λ1 β1 Sh − λ1 β1 Sh u1 − λ2 β1 Sh + λ2 β1 Sh u1 + µ2 Sh λ1 + λ2 γh + λ2 δh + λ2 µ1 + λ2 µ2 Nh + µ2 Ih λ2 − λ3 γh + λ3 µ2 Rh + λ4 β3 Sv − λ4 β3 Sv u2 − λ5 β3 Sv + λ5 β3 Sv u2 = −αh λ1 − A1 + β1 (λ1 − λ2 )(1 − u1 )Sh + µ2 λ1 Sh + (γh + δh )λ2 + (µ1 + µ2 Nh )λ2 + µ2 λ2 Ih − γh λ3 + µ2 λ3 Rh + β3 (λ4 − λ5 )(1 − u2 )Sv
80 0
∂H ∂Rh = −λ1 αh + µ2 Sh λ1 + µ2 Ih λ2 + µ1 λ3 + Sh µ2 λ3
λ3 = −
+ 2µ2 Rh λ3 + µ2 Ih λ3 = −λ1 αh + µ2 Sh λ1 + µ2 Ih λ2 + µ1 λ3 + Sh µ2 λ3 + µ2 Rh λ3 + µ2 Rh λ3 + µ2 Ih λ3 = −λ1 αh + µ2 Sh λ1 + µ2 Ih λ2 + µ1 λ3 + Sh µ2 λ3 + µ2 Rh λ3 + µ2 Ih λ3 + µ2 Rh λ3 = −λ1 αh + µ2 Sh λ1 + µ2 Ih λ2 + µ1 λ3 + (Sh + Rh + Ih )µ2 λ3 = +µ2 Rh λ3 − λ1 αh + µ2 Sh λ1 + µ2 Ih λ2 + µ1 λ3 + Nh µ2 λ3 = +µ2 Rh λ3 − λ1 αh + µ2 Sh λ1 + µ2 Ih λ2 + (µ1 + µ2 Nh )λ3 + µ2 Rh λ3 ∂H 0 λ4 = − ∂Sv = −A2 − b2 λ4 + b2 u3 λ4 + β3 Ih λ4 − β3 Ih u2 λ4 + µ3 λ4 + 2µ4 Sv λ4 + µ4 Iv λ4 + r0 u3 λ4 − β3 Ih λ5 + β3 Ih u2 λ5 + µ4 Iv λ5 = −A2 − b2 λ4 + b2 u3 λ4 + β3 Ih λ4 − β3 Ih u2 λ4 − β3 Ih λ5 + β3 Ih u2 λ5 + µ3 λ4 + µ4 Sv λ4 + µ4 Iv λ4 + µ4 Sv λ4 + r0 u3 λ4 + µ4 Iv λ5 = −A2 − b2 λ4 + b2 u3 λ4 + β3 (λ4 − λ5 )(1 − u2 )Ih + µ3 λ4 + µ4 λ4 (Sv + Iv ) + µ4 Sv λ4 + r0 u3 λ4 + µ4 Iv λ5 = −A2 − b2 λ4 + b2 u3 λ4 + β3 (λ4 − λ5 )(1 − u2 )Ih + µ3 λ4 + µ4 λ4 Nv + µ4 Sv λ4 + r0 u3 λ4 + µ4 Iv λ5
81 = −A2 − b2 (1 − u3 ) + β3 (λ4 − λ5 )(1 − u2 )Ih + (µ3 + µ4 Nv )λ4 + µ4 Sv λ4 + r0 u3 λ4 + µ4 Iv λ5 ∂H 0 λ5 = − ∂Iv = −A2 + β2 Sh λ1 − β2 Sh u2 λ1 − β2 Sh λ2 + β2 Sh u2 λ2 − b2 λ4 + b2 u3 λ4 + µ4 Sv λ4 + δv λ5 + µ3 λ5 + µ4 Sv λ5 + 2µ4 Iv λ5 + r0 u3 λ5 = −A2 + β2 Sh λ1 − β2 Sh u2 λ1 − β2 Sh λ2 + β2 Sh u2 λ2 − b2 λ4 + b2 u3 λ4 + µ4 Sv λ4 + δv λ5 + µ3 λ5 + µ4 Sv λ5 + µ4 Iv λ5 + µ4 Iv λ5 + r0 u3 λ5 = −A2 + β2 Sh λ1 − β2 Sh u2 λ1 − β2 Sh λ2 + β2 Sh u2 λ2 − b2 λ4 + b2 u3 λ4 + µ4 Sv λ4 + δv λ5 + µ3 λ5 + µ4 (Sv + Iv )λ5 + µ4 Iv λ5 + r0 u3 λ5 = −A2 + β2 (λ1 − λ2 )(1 − u2 )Sh − b2 λ4 (1 − u3 ) + µ4 λ4 Sv + λ5 δv + (µ3 + µ4 Nv )λ5 + µ4 λ5 Iv + r0 u3 λ5 dengan kondisi transversality(kondisi batas) λi (T ) = 0,
i = 1, 2, ..., 5
Karena Sh (t) = Sh∗ (t), Ih (t) = Ih∗ (t), Rh (t) = Rh∗ (t), Sv (t) = Sv∗ (t), dan Iv (t) = Iv∗ (t) dan berdasarkan prinsip optimal ∂H ∂H ∂H = 0, ∂u = 0, dan ∂u = 0, sehingga didapat: ∂u 1 2 3 diperoleh karakteristik pengendalian optimal u∗1 , u∗2 , dan u∗3 yang diberikan sebagai berikut: 1. Untuk kendali optimal u∗1 ∂H = B1 u1 + λ1 β1 Sh Ih − λ2 β1 Sh Ih = 0 ∂u1 B1 u1 = −λ1 β1 Sh Ih + λ2 β1 Sh Ih
82 B1 u1 = (λ2 − λ1 )β1 Sh Ih (λ2 − λ1 )β1 Sh∗ Ih∗ u∗1 = B1 Karena 0 ≤ u∗1 ≤ 1, maka dapat dtulis (λ2 − λ1 )β1 Sh∗ Ih∗ ∗ u1 = max min ,1 ,0 B1 2. Untuk kendali optimal u∗2 ∂H = B2 u2 + λ1 β2 Sh Iv − λ2 β2 Sh Iv + λ4 β3 Sv Ih ∂u2 − λ5 β3 Sv Iv =0 B2 u2 = −λ1 β2 Sh Iv + λ2 β2 Sh Iv − λ4 β3 Sv Ih + λ5 β3 Sv Iv B2 u2 = β2 (λ2 − λ1 )Sh Iv + β3 (λ5 − λ4 )Sv Ih β2 (λ2 − λ1 )Sh∗ Iv∗ + β3 (λ5 − λ4 )Sv∗ Ih∗ u∗2 = B2 Karena 0 ≤ u∗2 ≤ 1, maka dapat ditulis β2 (λ2 − λ1 )Sh∗ Iv∗ + β3 (λ5 − λ4 )Sv∗ Ih∗ ∗ u2 = max min ,1 ,0 B2 3. Untuk kendali optimal u∗3 ∂H = B3 u3 − λ4 b2 Sv − λ4 b2 Iv − λ4 r0 Sv ∂u3 − λ5 r0 Iv = 0 B3 u3 = λ4 b2 Sv + λ4 b2 Iv + λ4 r0 Sv + λ5 r0 Iv B3 u3 = b2 λ4 (Sv + Iv ) + r0 (λ4 Sv + λ5 Iv ) b2 λ4 Nv∗ + r0 (λ4 Sv∗ + λ5 Iv∗ ) u∗3 = B3
83 Karena 0 ≤ u∗3 ≤ 1, maka dapat dtulis u∗3
4.10
b2 λ4 Nv∗ + r0 (λ4 Sv∗ + λ5 Iv∗ ) = max min ,1 ,0 B3
Solusi Numerik
Pada bagian ini digunakan metode iterasi untuk mencari solusi numerik dari masalah kendali optimal. Algoritma numerik yang diberikan adalah metode beda hingga semi-implisit. Terlebih dahulu mendiskritkan interval [t0 , tf ] pada titik-titik ti = t0 + il, (i = 0, 1, · · · , n), dengan l adalah beda waktu sehingga tn = tf . Selanjutnya, didefinisikan variabel keadaan dan adjoint Sh (t), Ih (t), Rh (t), Sv (t), Iv (t), λ1 (t), λ2 (t), λ3 (t), λ4 (t), λ5 (t) dan kendali u∗1 , u∗2 dan u∗3 dalam Sh∗ , Ih∗ , Rh∗ , Sv∗ , Iv∗ , λ∗1 , λ∗2 , λ∗3 , λ∗4 , λ∗5 . Selanjutnya dikombinasikan pendekatan beda maju dan mundur. Dalam hal ini persamaan statenya diselesaikan dengan menggunakan metode beda hingga maju adalah sebagai berikut [1] Shi+1 − Shi = b1 + αh (Shi+1 + Ihi + Rhi ) − β1 Shi+1 Ihi (1 − ui1 ) l − β2 Shi+1 Ivi (1 − ui2 ) − [µ1 + µ2 (Shi+1 + Ihi + Rhi )]Shi+1 , Ihi+1 − Ihi = β1 Shi+1 Ihi+1 (1 − ui1 ) + β2 Shi+1 Ivi (1 − ui2 ) l − γh Ih (i + 1) − δh Ih (i + 1) − [µ1 + µ2 (Shi+1 + Ihi+1 + Rhi )]Ihi+1 , Rhi+1 − Rhi = γh Ih (i + 1) − [µ1 + µ2 (Shi+1 + Ihi+1 l
84 + Rhi+1 )]Rhi+1 , Svi+1 − Svi = b2 (Svi+1 + Ivi )(1 − ui3 ) − β3 Svi+1 Ihi+1 (1 − ui2 ) l − (µ3 + µ4 (Svi+1 + Ivi ))Svi+1 − r0 ui3 Ivi+1 , Ivi+1 − Ivi = β3 Svi+1 Ihi+1 (1 − ui2 ) − δv Iv (i + 1) l − (µ3 + µ4 (Svi+1 + Ivi+1 ))Ivi+1 − r0 ui3 Ivi+1 . Persamaan adjoint diselesaikan dengan menggunakan metode beda hingga mundur adalah sebagai berikut − λn−i−1 λn−i 1 1 = −αh λn−i−1 + (λn−i−1 − λn−i 1 1 2 )[β1 (1 l − ui1 )Ihi+1 + β2 (1 − ui2 )Ivi+1 ] Shi+1 + µ2 λn−i−1 + (µ1 + µ2 Nhi+1 )λn−i−1 1 1 i+1 i+1 + λn−i + µ2 λn−i 3 µ2 Rh , 2 Ih
λn−i − λn−i−1 2 2 = −αh λn−i−1 − A1 + β1 (λn−i−1 1 1 l − λn−i−1 )(1 − ui1 )Shi+1 + µ2 λn−i−1 Shi+1 2 1 + (γh + δh )λn−i−1 + (µ1 + µ2 Nhi+1 )λn−i−1 2 2 i+1 + µ2 λn−i−1 Ihi+1 − γh λn−i + µ2 λn−i 2 3 3 Rh i i+1 + β3 (λn−i − λn−i 4 5 )(1 − u2 )Sv ,
λn−i − λn−i−1 3 3 = −αh λn−i−1 + µ2 λn−i−1 Shi+1 1 1 l + µ2 λn−i−1 Ihi+1 + (µ1 + µ2 Nhi+1 )λn−i−1 2 3 + µ2 λn−i−1 Rhi+1 , 3 λn−i − λn−i−1 4 4 = −A2 − b2 λn−i−1 (1 − ui3 ) + β3 (λn−i−1 4 4 l
85 i i+1 − λn−i + (µ3 + µ4 Nvi+1 )λn−i−1 5 )(1 − u2 )Ih 4 i+1 ui3 + µ4 λn−i Svi+1 + r0 λn−i−1 + µ4 λn−i−1 5 Iv , 4 4
− λn−i−1 λn−i 5 5 )(1 − ui2 )Shi+1 − λn−i−1 = −A2 + β2 (λn−i−1 2 1 l Svi+1 (1 − u∗3 ) + µ4 λn−i−1 − b2 λn−i−1 4 4 + λn−i−1 δv + (µ3 + µ4 Nvi+1 )λn−i−1 5 5 ui3 Ivi+1 + r0 λn−i−1 + µ4 λn−i−1 5 5 Algoritma menjelaskan metode pendekatan untuk mendapatkan kendali optimal adalah sebagai berikut: Algoritma pada program Langkah 1: Sh (0) = Sh0 , Ih (0) = Ih0 , Rh (0) = Rh0 , Sv (0) = Sv0 , Iv (0) = Iv0 , λi (tf ) = 0(i = 1, ..., 5), u1 (0) = u2 (0) = u3 (0) = 0 Langkah 2: for i=1,...,n-1,do:
Shi+1
−1 − l −αh + β1 Ihi (1 − ui1 ) = 2lµ2 i l β2 Iv (1 − ui2 ) + µ1 + µ2 (Ihi + Rhi ) − 2lµ2 1 + {[1 + l(−αh + β1 Ihi (1 − ui1 ) 2lµ2 + β2 Ivi (1 − ui2 ) + µ1 + µ2 (Ihi + Rhi ))]2 1
+ 4lµ2 [Shi + lb1 + lαh (Ihi + Rhi )]} 2 , −1 − l −β1 Shi+1 (1 − ui1 ) + γh + δh i+1 Ih = 2lµ2
86 l µ1 + µ2 (Shi+1 + Rhi ) 1 − + {[1 2lµ2 2lµ2 + l(−β1 Shi+1 (1 − ui1 ) + γh + δh + µ1 + µ2 (Shi+1 + Rhi ))]2 + 4lµ2 [Ihi 1
+ lβ2 Shi+1 Ivi (1 − ui2 )]} 2 , −1 − l µ1 + µ2 (Shi+1 + Ihi+1 ) i+1 Rh = 2lµ2 2 1 + { l + lµ1 + lµ2 (Shi+1 + Ihi+1 ) 2lµ2 1 + 4lµ2 Rhi + lγh Ihi+1 } 2 , −1 − l[−b2 (1 − ui3 ) + β3 Ihi+1 (1 − ui2 ) 2lµ4 i l µ3 + µ4 Iv + r0 ui3 1 − + {[1 2lµ4 2lµ4 + l(−b2 (1 − ui3 ) + β3 Ihi+1 (1 − ui2 ) + µ3
Svi+1 =
1
+ µ4 Ivi + r0 ui3 )]2 + 4lµ4 (Svi + lb2 Ivi (1 − ui3 )]} 2 , −1 − l(δv + µ3 + µ4 Svi+1 + r0 ui3 ) 2lµ4 2 1 + { l + l(δv + µ3 + µ4 Svi+1 + r0 ui3 ) 2lµ4
Ivi+1 =
1
+ 4lµ4 (Ivi + lβ3 Svi+1 Ivi+1 (1 − ui2 ))} 2 , i i+1 λn−i−1 = {λn−i + l[λn−i 1 1 2 β1 (1 − u1 )Ih i i+1 + λn−i 2 β2 (1 − u2 )Iv
− µ2 Ihi+1 λn−i − µ2 Rhi+1 λn−i 2 3 ]} /{1 + l[β1 (1 − ui1 )Ihi+1 + β2 (1 − ui2 )Ivi+1 + µ2 Shi+1 + µ1 + µ2 Nhi+1 − αh ]}, λn−i−1 = {λn−i + l[αh λn−i−1 2 2 1
87 + A1 − β1 λn−i−1 (1 − ui1 )Shi+1 − µ2 Shi+1 λn−i−1 1 1 i+1 − µ2 λn−i + γh λn−i 3 Rh 3 i i+1 − λn−i − β3 (λn−i 5 )(1 − u2 )Sv } 4
/{1 + l[γh + δh + µ1 + µ2 Nhi+1 + µ2 Ihi+1 − β1 (1 − ui1 )Shi+1 ]}, − µ2 Shi+1 λn−i−1 + l[αh λn−i−1 = {λn−i λn−i−1 1 1 3 3 ]}/{1 + l[µ1 − µ2 Ihi+1 λn−i−1 2 + µ2 Nhi+1 + µ2 Rhi+1 ]}, i+1 + l[A2 + β3 (1 − ui2 )λn−i = {λn−i λn−i−1 5 Ih 4 4 i i+1 + µ3 − µ4 Ivi+1 λn−i 5 ]}/{1 + l[β3 (1 − u2 )Ih
− µ4 Nvi+1 + µ4 Svi+1 + r0 ui3 − b2 (1 − ui3 )]}, Svi+1 (1 − ui3 ) − µ4 λn−i−1 + l[A2 + b2 λn−i−1 = {λn−i λn−i−1 4 4 5 5 )(1 − ui2 )Shi+1 ]} − λn−i−1 − β2 (λn−i−1 2 1 /{1 + l[δv + µ3 + µ4 Nvi+1 + µ4 Ivi+1 + r0 ui3 ]}, R1i+1 R2i+1
)β1 Shi+1 Ihi+1 − λn−i−1 (λn−i−1 1 2 = , B1 = {(λn−i−1 − λn−i−1 )β2 Shi+1 Ivi+1 2 1 + (λn−i−1 − λn−i−1 )β3 Svi+1 Ihi+1 }/B2 , 5 4
R3i+1 = {b2 λn−i−1 Nvi+1 + r0 (λn−i−1 Svi+1 4 4 + λn−i−1 Ivi+1 )}/B3 , 5 ui+1 = min 1, max(R1i+1 , 0) , 1 ui+1 = min 1, max(R2i+1 , 0) , 2 ui+1 = min 1, max(R3i+1 , 0) , 3 end for Langkah 3:
88 for i=1,...,n-1,tulis Sh∗ (ti ) = Sh∗ , Ih∗ (ti ) = Ih∗ , Rh∗ (ti ) = Rh∗ , Sv∗ (ti ) = Sv∗ , Iv∗ (ti ) = Iv∗ , u∗1 (ti ) = ui1 , u∗2 (ti ) = ui2 , u∗3 (ti ) = ui3 end for Untuk perkembangan penyakit sesudah dikendalikan, akan dilakukan simulasi model. Adapun nilai parameter yang digunakan seperti yang terdapat pada tabel berikut ini : Tabel 4.3: Input Parameter Parameter b1 αh b2 µ1 µ2 µ3 µ4 µh
Nilai 25/hari 0.03/hari 0.4/hari 0.4/hari 0.00002/hari 0.15/hari 2.8 × 10(−3) /hari 10/hari
Parameter µv δh δv γh β1 β2 β3 r0
Nilai 0.02/hari 0.03/hari 0.04/hari 0.00037/hari 0.0004 0.0006 0.009 0.002
Pada Gambar 4.4 menunjukkan adanya populasi individu yang terinfeksi dengan menggunakan kendali. Dalam grafik simulasi tersebut terlihat bahwa populasi individu yang terinfeksi semakin menurun ketika kendali dilakukan. Dengan demikian, jumlah orang yang terinfeksi setelah kendali semakin kecil. Pada Gambar 4.5 menunjukkan adanya jumlah populasi vektor dengan menggunakan kendali. Dalam grafik simulasi tersebut terlihat bahwa jumlah populasi vektor semakin menurun ketika kendali dilakukan. Dengan demikian, jumlah populasi vektor setelah kendali semakin kecil.
89
Grafik Ih 50 45 40 35
h
I (t)
30 25 20 15 10 5 0
0
20
40
60
80
100
Time(days)
Gambar 4.4: Simulasi dari populasi individu yang terinfeksi dengan kendali
90
Grafik Nv 100 90 80 70
v
N (t)
60 50 40 30 20 10 0
0
20
40
60
80
100
Time(days)
Gambar 4.5: Simulasi dari jumlah populasi vektor dengan kendali
BAB V PENUTUP
Pada bab ini diberikan kesimpulan sebagai hasil dari analisa model yang telah diperoleh dan saran sebagai pertimbangan dalam pengembangan atau penelitian lebih lanjut. 5.1
Kesimpulan
Berdasarkan analisis dan pembahasan yang telah disajikan pada bab sebelumnya, dapat disimpulkan beberapa hal sebagai berikut : 1. Diperoleh titik kesetimbangan dan kestabilan lokal sebagai berikut:. Titik kesetimbangan bebas penyakit b1 b1 b2 0 0 0 0 0 Ef = Sh , Ih , Nh , Sv , Iv = , 0, , , 0 µh µh µv dan Titik kesetimbangan endemik E1 = (Sh∗ , Ih∗ , Nh∗ , Sv∗ , Iv∗ ) dengan b1 − (γh + δh + µh )Ih∗ b1 − δh Ih∗ , Nh∗ = µh µh ∗) b (b − δ I 2 1 h h Sv∗ = µh β3 Ih∗ + µv (b1 − δh Ih∗ ) µh b2 β3 Ih∗ Iv∗ = (µv + δv )(µh β3 Ih∗ + µv (b1 − δh Ih∗ ))
Sh∗ =
91
92 g(Ih ) = AIh∗2 + BIh∗ + C = 0 dimana A = K1 K2 (µh β3 (β1 − δh ) + µv δh (δh − β1 )), B = K1 K2 (µv b1 β1 − 2µv b1 δh + µh b1 β3 ) + K2 b1 β1 (µv δh − µh β3 ) + K1 µh b2 β2 β3 , C = µv b21 K1 K2 (1 − R0 ), sehingga Ih∗1 Ih∗2
√
B 2 − 4AC √2A −B + B 2 − 4AC = 2A
=
−B −
dan bilangan reproduksi dasar adalah R0 =
β2 β3 b2 µh β1 + (µv b1 )(δv + µv )(γh + δh + µh ) (γh + δh + µh )
Titik kesetimbangan bebas penyakit bersifat stabil untuk R0 < 1 dan titik kesetimbangan endemik bersifat stabil untuk R0 > 1. 2. Perubahan jenis kurva bifurkasi (maju atau mundur) dipengaruhi oleh perubahan nilai R0 yang mempengaruhi nilai A,B, dan C sehingga nilai titik puncaknya pun berubah. Untuk bifurkasi maju terjadi pada saat titik puncak dari sistem persamaan g(Ih ) B = 1. Dalam hal ini terdapat yaitu pada saat Ih = − 2A tiga kemungkinan untuk titik kesetimbangannya, yaitu a. Terdapat dua titik endemik yaitu titik endemik bersifat stabil dan titik endemik bersifat tak stabil. b. Terdapat satu titik endemik dan bersifat stabil.
93 c. Tidak ada titik keseimbangan endemik sehingga menyebabkan timbulnya titik bebas penyakit. 3. Diperoleh kendali optimal sebagai berikut: (λ2 − λ1 )β1 Sh∗ Ih∗ u∗1 = max min ,1 ,0 B1 β2 (λ2 − λ1 )Sh∗ Iv∗ + β3 (λ5 − λ4 )Sv∗ Ih∗ ∗ u2 = max min ,1 ,0 B2 b2 λ4 Nv∗ + r0 (λ4 Sv∗ + λ5 Iv∗ ) ∗ ,1 ,0 u3 = max min B3 4. Hasil perhitungan numerik menunjukkan keefektifan kendali yang ada bisa mengurangi populasi individu yang terinfeksi dan jumlah populasi vektor. 5.2
Saran
Pada Tugas Akhir ini tidak dibahas mengenai analisa kestabilan global untuk dapat mengetahui perilaku dinamis secara global, maka untuk selanjutnya bisa dilakukan analisa kestabilan global. Selain itu, dapat dapat dicari masalah kestabilan dan kendali optimal untuk penyakit yang lain.
Halaman ini sengaja dikosongkan.
LAMPIRAN A Kurva Bifurkasi Mundur
clear all; clc; close all; b_1=3; b_2=12; beta_1=4; beta_2=16; beta_3=1; mu_h=1; mu_v=4; gamma_h=2; delta_h=78; delta_v=2; K1 = gamma_h+delta_h+mu_h; K2 = mu_v+delta_v; A=K1*K2*(mu_h*beta_3*(beta_1-delta_h)+... (mu_v*delta_h)*(delta_h-beta_1)); A B=K1*K2*(mu_v*b_1*beta_1-2*mu_v*b_1*delta_h+... mu_h*b_1*beta_3)+(K2*b_1*beta_1)*... (mu_v*delta_h-mu_h*beta_3)+... K1*mu_h*b_2*beta_2*beta_3; B
97
98 r_0=(-B^2+4*A*(mu_v*b_1*b_1*K1*K2))/... (4*A*(mu_v*b_1*b_1*K1*K2)); r_0 R_0=r_0:(1-r_0)/100:1; for i=1:1:101 C=mu_v*b_1*b_1*K1*K2*(1-R_0(i)); I_1(i)=(-B+sqrt(B^2-4*A*C))/(2*A); I_2(i)=(-B-sqrt(B^2-4*A*C))/(2*A); end plot(R_0,I_2,’r’,R_0,I_1,’b’, ’LineWidth’,2); R_a=1:(3-1)/100:3; for i=1:1:101 C=mu_v*b_1*b_1*K1*K2*(1-R_a(i)); I_3(i)=(-B+sqrt(B^2-4*A*C))/(2*A);
end hold on plot(R_a,I_3,’b’, ’LineWidth’,2); hold off xlabel(’Bilangan Reproduksi Dasar (R0)’) ylabel(’Populasi (Ih)’) legend(’I_1’,’I_2’,20,’location’,’eastoutside’); grid on axis([r_0-1,R_a(101)+2,0,I_3(101)]);
LAMPIRAN B Kurva Bifurkasi Maju
clear all; clc; close all; b_1=1.7; b_2=2; beta_1=1; beta_2=4; beta_3=0.48; mu_h=0.6; mu_v=0.8; gamma_h=0.4; delta_h=1.4; delta_v=0.8; K1 = gamma_h+delta_h+mu_h; K2 = mu_v+delta_v; A=K1*K2*(mu_h*beta_3*(beta_1-delta_h)+... (mu_v*delta_h)*(delta_h-beta_1)); A B=K1*K2*(mu_v*b_1*beta_1-2*mu_v*b_1*delta_h+... mu_h*b_1*beta_3)+(K2*b_1*beta_1)*... (mu_v*delta_h-mu_h*beta_3)+... K1*mu_h*b_2*beta_2*beta_3; B
99
100 r_0=(-B^2+4*A*(mu_v*b_1*b_1*K1*K2))/... (4*A*(mu_v*b_1*b_1*K1*K2)); r_0 R_0=r_0:(1-r_0)/100:1; for i=1:1:101 C=mu_v*b_1*b_1*K1*K2*(1-R_0(i)); I_1(i)=(-B+sqrt(B^2-4*A*C))/(2*A); I_2(i)=(-B-sqrt(B^2-4*A*C))/(2*A); end plot(R_0,I_2,’r’,R_0,I_1,’b’, ’LineWidth’,2); R_a=1:(3-1)/100:3; for i=1:1:101 C=mu_v*b_1*b_1*K1*K2*(1-R_a(i)); I_3(i)=(-B+sqrt(B^2-4*A*C))/(2*A);
end hold on plot(R_a,I_3,’b’, ’LineWidth’,2); hold off xlabel(’Bilangan Reproduksi Dasar (R0)’) ylabel(’Populasi (Ih)’) legend(’I_1’,’I_2’,20,’location’,’eastoutside’); grid on axis([r_0-1,R_a(101)+2,0,I_3(101)]);
LAMPIRAN C Simulasi dari populasi individu yang terinfeksi dan Simulasi dari jumlah populasi vektor dengan kendali
clear all; clc; close all; b_1=25; alpha_h=0.03; b_2=0.4; mu_1=0.4; mu_2=0.00002; mu_3=0.15; mu_4=2.8*10^(-3); mu_h=10; mu_v=0.02; delta_h=0.03; delta_v=0.04; gama_h=0.00037; beta_1=0.0004; beta_2=0.0006; beta_3=0.009; r_0=0.002; n=100; S_h0=50; I_h0=50; R_h0=50; 101
102 S_v0=50; I_v0=50; l=1; A_1=50; A_2=50; B_1=10; B_2=30; B_3=30; u_1(1)=0; u_2(1)=0; u_3(1)=0; lamda_1(n) = 0; lamda_2(n) = 0; lamda_3(n) = 0; lamda_4(n) = 0; lamda_5(n) = 0; S_h=zeros(n,1); I_h=zeros(n,1); R_h=zeros(n,1); I_v=zeros(n,1); S_v=zeros(n,1); S_h(1)=S_h0; I_h(1)=I_h0; R_h(1)=R_h0; S_v(1)=S_v0; I_v(1)=I_v0; N_h(1)=S_h(1)+I_h(1)+R_h(1); N_v(1)=S_v(1)+I_v(1); %% Dengan Kendali N=n-1; for i=1:N-1
103 S_h(i+1)=(-1-l*((-alpha_h)+(beta_1*I_h(i)*... (1-u_1(i))))-(l*(((beta_2*I_v(i))*... (1-u_2(i)))+mu_1+(mu_2*(I_h(i)+... R_h(i)))))+((1+(l*(-alpha_h+... (beta_1*I_h(i)*(1-u_1(i)))+(beta_2*... I_v(i)*(1-u_2(i)))+mu_1+(mu_2*... (I_h(i)+R_h(i))))))^2+(4*l*mu_2)*... (S_h(i)+(l*b_1)+((l*alpha_h)*... (I_h(i)+R_h(i)))))^(1/2))/(2*l*mu_2); I_h(i+1)=(-1-(l*(-beta_1*S_h(i+1)*(1-u_1(i))+... gama_h+delta_h))-(l*(mu_1+mu_2*... (S_h(i+1)+R_h(i))))+((1+(l*(-beta_1*... S_h(i+1)*(1-u_1(i))+gama_h+delta_h+... mu_1+(mu_2*(S_h(i+1)+R_h(i))))))^2+... (4*l*mu_2)*(I_h(i)+(l*beta_2*S_h(i+1)*... I_v(i))*(1-u_2(i))))^(1/2))/(2*l*mu_2); R_h(i+1)=(-1-(l*(mu_1+mu_2*(S_h(i+1)+I_h(i+... 1))))+((1+l*mu_1+(l*mu_2)*(S_h(i+1)+... I_h(i+1)))^2+(4*l*mu_2)*(R_h(i)+... (l*gama_h*I_h(i+1))))^(1/2))/(2*l*mu_2); S_v(i+1)=(-1-l*(-b_2*(1-u_3(i))+(beta_3*I_h(i+1))*... (1-u_2(i)))-l*(mu_3+(mu_4*I_v(i))+(r_0*... u_3(i)))+((1+l*(-b_2*(1-u_3(i))+(beta_3*... I_h(i+1))*(1-u_2(i))+mu_3+(mu_4*I_v(i))+... r_0*u_3(i)))^2+(4*l*mu_4)*(S_v(i)+(l*b_2*... I_v(i))*(1-u_3(i))))^(1/2))/(2*l*mu_4); I_v(i+1)=(-1-l*(delta_v+mu_3+(mu_4*S_v(i+1))+... (r_0*u_3(i)))+((1+l*(delta_v+mu_3+(mu_4*... S_v(i+1))+(r_0*u_3(i))))^2+(4*l*mu_4)*... (I_v(i)+(l*beta_3*S_v(i+1)*I_v(i+1))*... (1-u_2(i))))^(1/2))/(2*l*mu_4); N_h(i+1)=S_h(i+1)+I_h(i+1)+R_h(i+1); N_v(i+1)=S_v(i+1)+I_v(i+1);
104
lamda_1(n-i-1)=(lamda_1(n-i)+(l*(lamda_2(n-i)*... beta_1*(1-u_1(i))*I_h(i+1))+... (lamda_2(n-i)*beta_2*(1-u_2(i))*... I_v(i+1))-(mu_2*I_h(i+1)*lamda_2(n-i))-... (mu_2*R_h(i+1)*lamda_3(n-i))))/... (1+l*((beta_1*(1-u_1(i))*I_h(i+1))+... (beta_2*(1-u_2(i))*I_v(i+1))+... (mu_2*S_h(i+1))+mu_1+... (mu_2*N_h(i+1))-alpha_h)); lamda_2(n-i-1)=(lamda_2(n-i)+l*((alpha_h*lamda_1... (n-i-1))+A_1-(beta_1*lamda_1(n-i-1)*... (1-u_1(i))*S_h(i+1))-(mu_2*S_h(i+1)*... lamda_1(n-i-1))+(gama_h*lamda_3(n-i))-... (mu_2*lamda_3(n-i)*R_h(i+1))-(beta_3*... (lamda_4(n-i)-lamda_5(n-i))*(1-u_2(i))*... S_v(i+1))))/(1+l*(gama_h+delta_h+mu_1+... (mu_2*N_h(i+1))+(mu_2*I_h(i+1))-... (beta_1*(1-u_1(i))*S_h(i+1)))); lamda_3(n-i-1)=(lamda_3(n-i)+l*((alpha_h*lamda_1... (n-i-1))-(mu_2*S_h(i+1)*lamda_1... (n-i-1))-(mu_2*I_h(i+1)*lamda_2... (n-i-1))))/(1+l*(mu_1+... (mu_2*N_h(i+1))+(mu_2*R_h(i+1)))); lamda_4(n-i-1)=(lamda_4(n-i)+l*(A_2+(beta_3*... (1-u_2(i))*lamda_5(n-i)*I_h(i+1))-... (mu_4*I_v(i+1)*lamda_5(n-i))))/... (1+l*((beta_3*(1-u_2(i))*... I_h(i+1))+mu_3+(mu_4*N_v(i+1))+... (mu_4*S_v(i+1))+(r_0*u_3(i))-... (b_2*(1-u_3(i))))); lamda_5(n-i-1)=(lamda_5(n-i)+l*(A_2+(b_2*lamda_4... (n-i-1)*(1-u_3(i)))-(mu_4*lamda_4...
105 (n-i-1)*S_v(i+1))-(beta_2*(lamda_1... (n-i-1)-lamda_2(n-i-1))*(1-u_2(i))*... S_h(i+1))))/(1+l*(delta_v+mu_3+... (mu_4*N_v(i+1))+(mu_4*I_v(i+1))+... (r_0*u_3(i)))); R_1(i+1)=((lamda_2(n-i-1)-lamda_1(n-i-1))*... beta_1*S_h(i+1)*I_h(i+1))/B_1; R_2(i+1)=((lamda_2(n-i-1)-lamda_1(n-i-1))*... beta_2*S_h(i+1)*I_v(i+1)+(lamda_5... (n-i-1)-(n-i-1))*beta_3*S_v(i+1)*... I_h(i+1))/B_2; R_3(i+1)=((b_2*lamda_4(n-i-1)*N_v(i+1))+... r_0*(lamda_4(n-i-1)*S_v(i+1)+... lamda_5(n-i-1)*I_v(i+1)))/B_3; u_1(i+1)=min(1,max(R_1(i+1),0)); u_2(i+1)=min(1,max(R_2(i+1),0)); u_3(i+1)=min(1,max(R_3(i+1),0)); end u_1 u_2 u_3 figure(3); plot([0:l:(n-1)*l],I_h,’b’); hold on title(’Grafik I_h’); xlabel(’Time(days)’); ylabel(’I_h(t)’); figure(4);
106 plot([0:l:(n-2)*l],N_v,’b’); hold on title(’Grafik N_v’); xlabel(’Time(days)’); ylabel(’N_v(t)’);
DAFTAR PUSTAKA
[1] Lashari. A. A, Hattaf. K, Zaman. G, dan Li. X, (2013). ”Backward bifurcation and optimal control of a vector borne disease”. Applied Mathematical and Information Sciences, No.1, 301-309 [2] Zaman. G, Khan. M. A, Islam. S, Chohan. M. I, dan Jung. I. H, (2012). ”Modeling Dynamical Interaction between Leptospirosis Infected Vector and Human Population”. Applied Mathematical Sciences, Vol. 6, 2012, no. 26, 1287-1302. [3] Iskandar. A. (2011). ”Pemberantasan serangga dan binatang pengganggu”. APKTS Pusdiknakes. Depkes RI. Jakarta [4] Artikel Kesehatan. (2014). ”Vector-borne diseases jadi tema hari kesehatan dunia”. [5] World Health Organization. (2014). ”Vector BorneDiseases”. [6] Zainal. F.D. (2014). ”Analisis Kestabilan pada Model Transmisi Virus Hepatitis B yang Dipengaruhi oleh Migrasi”. Tugas Akhir S1 Jurusan Matematika ITS Surabaya. [7] Iswahyuni. N. (2013). ”Analisa Kualitatif pada Model Epidemik SIS dengan Pengobatan”. Tugas Akhir S1 Jurusan Matematika ITS Surabaya. 95
96 [8] Naidu. S. D. (2002). ”Optimal Control System”. USA: CRC Press LLC. [9] Resmi, F. (2014). ”Kendali Optimal pada Sistem PreyPredator dengan Pemberian Makan Alternatif pada Predator”. Tugas Akhir S1 Jurusan Matematika ITS Surabaya [10] Fitria. I. (2014). ”Kendali Optimal dalam Produksi Sumber Energi Terbarukan dan Tidak Terbarukan”. Tugas Akhir S1 Jurusan Matematika ITS Surabaya [11] Londok, J., (2011). Leptospirosis. http://jcgirlonthemove.blogspot.com/2011/03/vectorborne-disease.html, diakses pada tanggal 5 Desember 2014. [12] Driessche. P. V. D dan Watmough. J, (2002). ”Reproduction Numbers and Sub-threshold Endemic Equilibria for Compartmental Models of Disease Transmission”. Mathematical Biosciences, No.180, 29-48 [13] Bryson, A. E. dan Ho, Y. C., (1975). Applied Optimal Control. New York: Taylor & Francis Group.
LAMPIRAN D Biodata Penulis
Penulis bernama Charisma Juni Kumalasari, lahir di Nganjuk, 11 Juni 1994. Penulis merupakan anak terakhir dari tiga bersaudara. Penulis menempuh pendidikan formal dimulai dari TK NawaKartika (1997-1999), SDN Dawu 2 Ngawi (1999-2005), SMP Negeri 1 Ngawi (2005-2008), dan SMA Negeri 2 Ngawi (2008-2011). Setelah lulus dari SMA, pada tahun 2011 penulis melanjutkan studi ke jenjang S1 di Jurusan Matematika ITS Surabaya melalui jalur Undangan dengan NRP 1211 100 032. Di Jurusan Matematika, penulis mengambil Bidang Minat Pemodelan dan Simulasi Sistem. Selain aktif kuliah, penulis juga aktif berorganisasi di KM ITS melalui HIMATIKA ITS sebagai staf Depart. Hubungan Luar (2012-2014) dan UKM Taekwondo ITS sebagai Sekretaris Umum (2012-2013). Informasi lebih lanjut mengenai Tugas Akhir ini dapat ditujukan ke penulis melalui email:
[email protected]
