TUGAS BESAR

MODEL MATEMATIKA SEICR PADA PENYEBARAN

HEPATITIS B DENGAN PEMBERIAN VAKSINASI

ANGGOTA KELOMPOK:

Julinar (06111740000025)
Karohmatul Amalia MS. (06111740000027)
Siti Masriyah (06111740000073)
Ilham Dwi Prasetyo (06111740000106)

Dosen Pengampu:
Dr. Dra. Mardlijah, MT
NIP. 19670114 199102 2 00

Mata Kuliah :
Matematika Sistem (Kelas B)

DEPARTEMEN MATEMATIKA
FAKULTAS SAINS DAN ANALITIKA DATA
INSTITUT TEKNOLOGI SEPULUH NOPEMBER
SURABAYA
2020
 BAB I

PENDAHULUAN

1.1 Latar Belakang

Secara langsung sistem bias dikatakan sebagai bagian dari realita. Realita di luar
sistem dinamakan sekitar sistem. Interaksi antara sistem dan sekitar sistem
direalisasikan lewat besaran yang biasanya merupakan fungsi waktu. Besaran ini disebut
masukan (input) dan keluaran (output). Masukan dan keluaran sistem yang disajikan
oleh fungsi waktu bias merupakan waktu yang kontinu atau distrit. Hal ini berkaitan
dengan apa yang dinamakan sistem kontinu dan sistem diskrit. Mengkaji (menganalisis)
proses fisis atau mendisainnya dinamakan sistem fisis dalam hal ini hubungan masukan
dan keluran sistem disajikan oleh suatu model matematika. Sering kali model
matematika ini berbentuk suatu persamaan differensial (untuk yang kontinu) dan
persamaan beda (untuk yang diskrit). Untuk sistem dengan masukan dan keluaran yang
disajikan oleh bentuk persamaan differensial biasa dinamakan sistem tergumpal
(lumped), bila tidak demikian dinamakan sistem terdistribusi. Pembentukan suatu model
matematika sering membutuhkan asumsi tentang sifat dasar proses fisis.

Pada makalah ini akan dibahas mengenai model dari sistem penyebaran virus
Hepatitis B menggunakan Model SEICR. Hepatitis merupakan salah satu jenis penyakit
menular. Tipe yang berbahaya dari hepatitis ialah Hepatitis B. Hal ini dikarenakan
penularan Hepatitis B terjadi sangat cepat melalui kontak tubuh, air liur, hubungan seks
dan sebagainya. Bahkan penyakit ini juga bisa diturunkan ke bayi yang baru lahir oleh
ibu yang terkena Hepatitis B. Penularan Hepatitis B dapat dimodelkan secara matematis
menggunakan pemodelan matematika. Model SEICR untuk Hepatitis B mengasumsikan
populasi menjadi lima sub-populasi yaitu: Susceptible (individu rentan), Exposed
(individu laten), Infected (individu tertular), Carrier (individu pembawa penyakit, dan
Recovery (individu sembuh). Adanya pembagian populasi ini bertujuan untuk
mengetahui tingkat penyebaran penyakit disetiap sub-populasi.
Model yang akan diteliti dalam penelitian ini ialah model SEICR yang merupakan
pengembangan dari model SACR (Suspectible, Acute, Infection, Carrier, Recovery).
Dimana dalam model ini ditambahkan sub populasi individu laten yaitu individu yang
 terkena virus namun tidak berpotensi untuk menularkan ke individu lain. . Kemudian
dalam model tersebut ditambahkan pula faktor kontrol vaksinasi dengan tujuan
mengontrol penyebaran Hepatitis B di suatu populasi dan menentukan tingkat vaksin
untuk menghentikan laju penyebaran.Perilaku dinamik yang diselidiki dari model
SEICR adalah mencari titik kesetimbangan, dilakukan pelinieran di sekitar titik
setimbang, mencari solusi sistem, menyelidiki sifat - sifat sistem yaitu kestabilan,
kekontrolan, dan keteramatan, membentuk feedback controller dan observernya, dan
melakukan simulasi gambar dari hasil sistem tersebut.

1.2 Rumusan Masalah

1. Bagaimana model sistem non-linier dari penyebaran Virus Hepatitis B dengan
vaksinasi?
2. Bagaimana titik setimbang dari model non-liner penyebaran Virus Hepatitis B
dengan vaksinasi?
3. Bagaimana proses linierisasi di titik setimbangnya?
4. Bagaimana solusi dari sistem penyebaran Virus Hepatitis B dengan vaksinasi?
5. Bagaimana sifat-sifat dari sistem penyebaran Virus Hepatitis B dengan vaksinasi?
6. Bagaimana feedback controller dan observer dari sistem penyebaran Virus
Hepatitis B dengan vaksinasi?
7. Bagaimana hasil simulasi yang didapat dari sistem penyebaran Virus Hepatitis B
dengan vaksinasi?

1.3 Tujuan
1. Mengetahui model non-linier dari penyebaran Virus Hepatitis B dengan vaksinasi.
2. Mengetahui titik setimbang dari model non-liner penyebaran Virus Hepatitis B
dengan vaksinasi.
3. Mengetahui proses linierisasi di titik setimbangnya.
4. Mengetahui solusi dari sistem penyebaran Virus Hepatitis B dengan vaksinasi.
5. Mengetahui sifat-sifat dari sistem penyebaran Virus Hepatitis B dengan vaksinasi.
6. Mengetahui feedback controller dan observer dari sistem penyebaran Virus
Hepatitis B dengan vaksinasi.
7. Mengetahui hasil simulasi yang didapat dari sistem penyebaran Virus Hepatitis B
dengan vaksinasi.
 BAB II

KAJIAN TEORI

2.1 Sistem Persamaan Diferensial Non Linier

Persamaan diferensial non linear adalah persamaan diferensial biasa yang tak linear.
Suatu persamaan diferensial dikatakan non linear jika persamaan diferensial tersebut
paling sedikit memenuhi satu dari krteria berikut ini (Ross, 2010:5) :
1. Memuat variable tak bebas dari turunan-turunanya berpangkat selain satu.
2. Terdapat perkalian dari variable tak bebas dan atau turunan-turunannya.
3. Terdapat fungsi transcendental dari variable tak bebas dan turunan-turunannya.

Diberikan sistem persamaan diferensial 𝑥̇ = 𝑓(𝑡, 𝑥) dengan :

𝑥1 (𝑡) 𝑓1 (𝑡, 𝑥1 , 𝑥2 , … , 𝑥𝑛 )

𝑥 = 𝑥2 (𝑡) 𝑑𝑎𝑛 𝑓(𝑡, 𝑥) = 𝑓2 (𝑡, 𝑥1 , 𝑥2 , … , 𝑥𝑛 ) (2.1)

⋮ ⋮
[𝑥𝑛 (𝑡)] [𝑓𝑛 (𝑡, 𝑥1 , 𝑥2 , … , 𝑥𝑛 )]

Sistem persamaan diatas adalah sistem persamaan diferensial non-linier dengan fungsi
non-linier dalam 𝑥1 , 𝑥2 , … , 𝑥𝑛 (Finizio & Ladas, 1971).

2.2 Nilai Eigen dan Vektor Eigen

Jika 𝐴 merupakan matriks 𝑛 × 𝑛, maka vektor tak nol 𝑥 dalam ℝ𝑛 disebut vektor eigen
dari 𝐴 jika 𝐴𝑥 merupakan kelipatan skalar dari 𝑥, yaitu :

𝐴𝑥 = 𝜆𝑥 (2.2)

Untuk skalar 𝜆. Skalar 𝜆 disebut nilai eigen dari 𝐴 dan 𝑥 dinamakan vektor yang
bersesuaian dengan skalar 𝜆.
Persamaan (2.2) dapat ditulis sebagai berikut :
𝐴𝑥 = 𝜆𝐼𝑥

(𝜆𝐼 − 𝐴)𝑥 = 0 (2.3)

Dengan I adalah matriks identitas. Persamaan (2.3) memiliki solusi tak nol jika dan
hanya jika

|𝜆𝐼 − 𝐴| = 0 (2.4)
 Persamaan (2.4) merupakan persamaan karakteristik dari matriks 𝐴 dan skalar yang
memenuhi persamaan (2.4) adalah nilai eigen dari 𝐴.

|𝜆𝐼 − 𝐴| = 𝜆𝑛 + 𝑐1 𝜆𝑛−1 + 𝑐2 𝜆𝑛−2 + ⋯ + 𝑐𝑛 (2.5)

Sehingga persamaan karakteristik dari 𝐴 menjadi :

𝜆𝑛 + 𝑐1 𝜆𝑛−1 + 𝑐2 𝜆𝑛−2 + ⋯ + 𝑐𝑛 = 0 (2.6)

Dengan 𝑐𝑖 ∈ ℝ, 𝑖 = 1,2,3, … , 𝑛

2.3 Titik Kesetimbangan

Titik kesetimbangan merupakan sebuah titik tetap yang tidak berubah. Diberikan
persamaan 𝑥̇ = 𝑓(𝑥), titik 𝑥̅ ∈ ℝ𝑛 disebut titik kesetimbangan apabila memenuhi
turunan pertama sama dengan nol, dituliskan dengan 𝑓(𝑥̅ ) = 0.

2.4 Kestabilan Titik Setimbang

Kondisi stabil dari titik kesetimbangan memiliki tiga kondisi kestabilan, yaitu: stabil,
stabil asimtotik, dan tidak stabil. Apabila diberikan persamaan diferensial tingkat 1,
yaitu 𝑥̇ (𝑡) = 𝑓(𝑥(𝑡)) dengan 𝑥 ∈ ℝ𝑛 . Penyelesaian keadaan awal 𝑥(0) = 𝑥0
dinotasikan oleh 𝑥(𝑡, 𝑥0 ) (Side, Sanusi, & Setiawan, 2016).
1. Suatu titik kesetimbangan 𝑥̅ dikatakan stabil untuk setiap 𝜖 > 0 terdapat 𝛿 > 0 dan
𝑡𝛿 sedemikian sehingga bila ‖𝑥𝑡𝛿 − 𝑥̅ ‖ < 𝛿 maka ‖𝑥(𝑡, 𝑥𝑡𝛿 ) − 𝑥̅ ‖ < 𝜖 untuk semua
𝑡 > 𝑡𝛿 .
2. 𝑥̅ dikatakan stabil asimtotik bila ia stabil dan bila ada 𝛿1 > 0 sehingga :
lim ‖𝑥(𝑡, 𝑥𝑡𝛿 ) − 𝑥̅ ‖ = 0 bila ‖𝑥𝑡𝛿 − 𝑥̅ ‖ < 𝛿1 .
𝑡→∞

3. Suatu titik kesetimbangan dikatakan tak stabil bila ia tidak stabil.

Secara umum, stabil berarti penyelesaiannya sangat dekat ke titik setimbang. Sedangkan
stabil asimtotik berarti penyelesaiannya konvergen ke titik setimbang (Susanto, 2012).

2.5 Linierisasi
Linearisasi merupakan suatu proses yang digunakan untuk membentuk suatu sistem
persamaan diferensial non linear menjadi sistem persamaan diferensial linear. Melalui
linearisasi tersebut dapat juga digunakan untuk menganalisis kestabilan. Pencarian hasil
linearisasi dari sistem persamaan diferensial non-linear dapat menggunakan matriks
Jacobian.
 Misalkan diberikan sistem persamaan diferensial non-linier berikut :
𝑥̇ = 𝑓(𝑥) (2.7)

Maka, linierisasai dari sistem persamaan (2.7) sebagai berikut :

𝑥̇ = 𝐽(𝑓(𝑥)) (2.8)

Dengan 𝐽 merupakan matriks Jacobian yang dapat dituliskan sebagai berikut :

𝜕𝑓1 𝜕𝑓1 𝑥1
⋯
𝑥1 𝑥𝑛
𝑥̇ = 𝐽(𝑓(𝑥))𝑥̇ = ⋮ ⋱ ⋮ 𝑥2

𝜕𝑓𝑚 𝜕𝑓𝑚 ⋮
⋯ [𝑥 ]
[ 𝑥1 𝑥𝑛 ] 𝑛

Berdasarkan matriks Jacobian 𝐽(𝑓(𝑥)) tersebut dapat dilakukan analisis kestabilan

disekitar titik kesetimbangan dengan nilai eigen (𝜆) (Olsder, 2003).

1. Apabila semua bagian real nilai eigen matriks 𝐽(𝑓(𝑥)) bernilai negatif, maka titik
ekuilirium dari sistem non-linier stabil asimtotik lokal.
2. Apabila terdapat paling sedikit satu nilai eigen matriks 𝐽(𝑓(𝑥)) yang bagian realnya
positif, maka titik kesetimbangan dari sistem non-linier tidak stabil.

2.6 Penyelesaian dari Persamaan Diferensial Linier

Persamaan diferensial varian-waktu berbentuk :

𝑥̇ (𝑡) = 𝐴𝑥(𝑡) + 𝐵𝑢(𝑡)

𝑦(𝑡) = 𝐶𝑥(𝑡) + 𝐷𝑢(𝑡)

Missal 𝑥̇ (𝑡) = 𝐴𝑥(𝑡), dengan 𝐴 adalah matriks konstan yang bersesuaian dan kondisi
awal 𝑥(𝑡0 ) = 𝑥0 , maka solusinya :

𝑥(𝑡) = 𝑒 𝐴(𝑡−𝑡0 ) ∙ 𝑥0

Dimana 𝑒 𝐴(𝑡−𝑡0 ) adalah matriks transisi. 𝑒 𝐴(𝑡−𝑡0 ) dapat dituliskan sebagai 𝑒 𝐴𝑡 . Matriks
𝑒 𝐴𝑡 merupakan matriks yang bergantung pada matriks 𝐴, dimana :

1. Jika matriks 𝐴 adalah matriks diagonal

Jika matriks 𝐴 adalah matriks diagonal, maka nilai eigennya adalah diagonalnya.
 Misalkan :
𝑎 0 0
𝐴 = [0 𝑏 0]
0 0 𝑐
Maka nilai eigennya adalah 𝜆1 = 𝑎, 𝜆2 = 𝑏 𝑑𝑎𝑛 𝜆3 = 𝑐. Sehingga dapat
didefinisikan :
𝜆1 0 ⋯ 0 𝑒 𝜆1 𝑡 0 ⋯ 0
0 𝜆2 ⋯ 0 𝜆2 𝑡
𝐴=[ ] dan 𝑒 𝐴𝑡
=[ 0 𝑒 ⋯ 0 ]
⋮ ⋮ ⋱ ⋮ ⋮ ⋮ ⋱ ⋮
0 0 ⋯ 𝜆𝑛 0 0 ⋯ 𝑒 𝜆𝑛𝑡
Maka penyelesaian atau solusi dari sistem adalah
𝑥(𝑡) = 𝑒 𝐴𝑡 ∙ 𝑥0

2. Jika matriks 𝐴 bukan matriks diagonal, tetapi dapat didiagonalkan

Ciri-ciri matriks 𝐴 yang dapat didiagonalkan adalah :
 𝐴 merupakan matriks persegi
 Mempunyai 𝑛 vektor eigen (vektor karakteristik yang bebas linier)
 Mempunyai 𝑛 nilai eigen (nilai karakteristik) yang berbeda yang
menghasilkan 𝑛 vektor karakteristik yang bebas linier.

Berikut adalah rumus untuk mendiagonalkan matriks 𝐴 :

𝐷 = 𝑇 −1 𝐴 𝑇

Dimana

𝑇 = [𝑣1 𝑣2 ⋯ 𝑣𝑛 ]

Dengan 𝑣1 , 𝑣2 , 𝑣𝑛 adalah vektor eigen. Sehingga matriks transisinya adalah :

𝑒 𝐴𝑡 = 𝑇 𝑒 𝐷𝑡 𝑇 −1

Maka penyelesaian atau solusi dari sistem adalah

𝑥(𝑡) = 𝑒 𝐷𝑡 ∙ 𝑥0

2.7 Sifat-Sifat Sistem

1. Kestabilan
Ada beberapa konsep kestabilan untuk persamaan differensial. Kestabilan ini
dibedakan menurut kestabilan sistem autonomus (berkaitan dengan vektor keadaan)
dan kestabilan yang dikaitkan dengan masukan dan keluaran sistem (kestabilan
didefinisikan dari segi masukan dan keluaran).
 Kestabilan dari segi nilai karakteristik
Diberikan persamaan differensial tingkat satu 𝑥̇ (𝑡) = 𝑓(𝑥(𝑡)) dengan 𝑥 ∈ 𝑅 𝑛 ,
penyelesaan dengan keadaan awal 𝑥(0) = 𝑥0 dinotasikan oleh 𝑥(𝑡, 𝑥0 ).
 Vektor 𝑥̅ yang memenuhi 𝑓(𝑥) = 0 disebut suatu titik setimbang.
 Suatu titik setimbang 𝑥 dikatakan stabil bila untuk setiap 𝜖 > 0 ada 𝛿 > 0 dan

𝑡𝛿 sedemikian hingga ||𝑥𝑡𝛿 − 𝑥̅ || < 𝛿 maka ||𝑥(𝑡, 𝑥𝑡𝛿 ) − 𝑥̅ || < 𝜖 untuk semua

𝑡 > 𝑡𝛿
 Suatu titik setimbang 𝑥̅ dikatakan stabil asimtotik bila ia stabil dan bila ada 𝛿1 > 0

sedemikian hingga lim |𝑥(𝑡, 𝑥𝑡𝛿 ) − 𝑥̅ | = 0 bila ||𝑥𝑡𝛿 − 𝑥̅ || < 𝛿1

𝑡→∞

 Suatu titik setimbang dikatakan takstabil bila ia tidak stabil.

Tanda || . || berarti norm, biasanya digunakan norm Euclidean. Secara intuisi

stabil berarti penyelesaian sangat dekat ketitik setimbang didalam suatu sekitar.
Sedangkan stabil asimtotik berarti penyelesaian konvergen ke titik setimbang
(asalkan titik awal adalah cukup dekat ke titik setimbang). Takstabil artinya selalu
ada penyelesaian yang dimulai dari manapun dekatnya dengan titik setimbang tapi
akhirnya menjauh dari titik setimbang.

Diberikan persamaan differensial 𝑥̇ = 𝐴𝑥dengan matriks 𝐴 berukuran n×n

dan mempunyai nilai karakteristik yang berbeda 𝜆1 , … , 𝜆𝑘 (𝑘 ≤ 𝑛)
 Titik asal 𝑥̅ = 0 adalah stabil asimtotik bila dan hanya bila 𝑅𝑒𝜆𝑖 < 0 untuk semua
𝑖 = 1, . . . , 𝑘.
 Titik asal 𝑥̅ = 0 adalah stabil bila dan hanya bila 𝑅𝑒𝜆𝑖 ≤ 0 untuk semua 𝑖 = 1,
. . . , 𝑘 dan untuk semua 𝜆𝑖 dengan 𝑅𝑒𝜆𝑖 = 0 multisiplisistas aljabar sama dengan
mutiplisistas geometrinya.
 Titik asal 𝑥̅ = 0 adalah takstabil bila dan hanya bila 𝑅𝑒𝜆𝑖 > 0 untuk beberapa
semua 𝑖 = 1, . . . , 𝑘 atau ada 𝜆𝑖 dengan 𝑅𝑒𝜆𝑖 = 0 dan multisiplisistas aljabar lebih
besar dari mutiplisitas geometrinya.
Kriteria Routh-Hurwitz
Kriteria Routh Hurwitz merupakan salah satu teknik yang digunakan untuk
menganalisis kestabilan. Analisis kestabilan dapat diperoleh dengan menentukan
nilai eigen dari persamaan karakteristik 𝑓(𝜆) = 𝜆𝑛 + 𝑐1 𝜆𝑛−1 + 𝑐2 𝜆𝑛−2 + ⋯ + 𝑐𝑛
Namun seringkali akar-akar dari persamaan karakteristik tersebut sulit ditentukan.
Sehingga Kriteria Routh-Hurwitz dapat digunakan sebagai salah satu cara untuk
menentukan nilai eigen tersebut (Nurfitriaa, Kiftiah, & Yudhi , 2019).
Diberikan suatu sistem persamaan karakteristik dalam bentuk polinomial berikut :
𝑐0 𝜆𝑛 + 𝑐1 𝜆𝑛−1 + 𝑐2 𝜆𝑛−2 + ⋯ + 𝑐𝑛−1 𝜆 + 𝑐𝑛 = 0
Sistem Persamaan tersebut akan memiliki kestabilan apabila keseluruhan dari
determinan Hurwitz adalah positif. Dimana dari Persamaan diatas dapat dibentuk
determinan Hurwitz (𝑑𝑒𝑡(𝐻)) sebagai berikut:
𝑐1 𝑐0 0 … 0
𝑐3 𝑐2 𝑐1 … 0
|
𝑑𝑒𝑡(𝐻) = 𝑐 𝑐4 𝑐3 … 0|
| .5 . . . .|
. . . . .
𝑐2𝑛−1 𝑐2𝑛−2 𝑐2𝑛−3 … 𝑐𝑛
Nilai-nilai untuk koefisien dengan indeks lebih besar dari n atau dengan indeks
negatif diganti dengan nol. Kondisi stabilitas akan terpenuhi jika:
∆𝑘 > 0, untuk 𝑘 = 1, 2, 3, … , 𝑛
dimana
∆1 = |𝑐1 |
𝑐1 𝑐0
∆2 = |𝑐 𝑐 |
3 2

𝑐1 𝑐0 0
∆3 = |𝑐3 𝑐2 𝑐1 |
𝑐5 𝑐4 𝑐3
Untuk 𝑘 = 𝑛, maka ∆𝑛 = det(𝐻)
Selain menggunakan determinan Hurwitz alternatif penyelesaian dalam
menganalisis kestabilan ialah menggunakan matriks Routh. Kriteria dari kestabilan
dilihat berdasarkan banyaknya perubahan tanda pada kolom pertama. Matriks
Routh dari Persamaan pertama adalah sebagai berikut:
 𝜆𝑛 𝑐0 𝑐2 𝑐4 𝑐6 …
𝜆𝑛−1 𝑐1 𝑐3 𝑐5 𝑐7 …
𝜆𝑛−2 𝑎 1 𝑎2 𝑎3 𝑎4 …
𝜆𝑛−3 𝑏 1 𝑏2 𝑏3 𝑏4 …
. = . .
. . .
𝜆 2 𝑑1 𝑑2
1 𝑒1
𝜆
[ 𝜆0 ] [ 𝑓1 ]
Dengan koefisien-koefisien:
𝑐1 𝑐2 −𝑐0 𝑐3 𝑎1 𝑐3 −𝑐1 𝑎2
𝑎1 = 𝑏1 =
𝑐1 𝑎1
𝑐1 𝑐4 −𝑐0 𝑐5 𝑎1 𝑐5 −𝑐1 𝑎3
𝑎2 = 𝑏2 =
𝑐1 𝑎1
𝑐1 𝑐6 −𝑐0 𝑐7 𝑎1 𝑐7 −𝑐1 𝑎4
𝑎3 = 𝑏3 =
𝑐1 𝑎1

dan seterusnya. Hal yang menjadi syarat kestabilan matriks Routh ialah berdasarkan
syarat perlu dan cukup sebagai berikut (Ningsih & Winarko , 2016)
1. Semua koefisien persamaan karakteristik positif.
2. Semua suku pada kolom pertama matriks Routh bertanda positif.

2. Kekontrolan
Diberikan Sistem: 𝑥̇ (𝑡) = 𝐴𝑥(𝑡) + 𝐵𝑢(𝑡) 𝑑𝑎𝑛 𝑦(𝑡) = 𝐶𝑥(𝑡) + 𝐷𝑢(𝑡). Sistem
linear tersebut dikatakan terkontrol bila untuk setiap kedaan sebarang 𝑥(0) = 𝑥0 ada
masukan (𝑡) yang tidak dibatasi mentransfer keadaan 𝑥0 kesebarang keadaan akhir
𝑥(𝑡1 ) = 𝑥1 dengan waktu akhir 𝑡1 hingga.
Artinya : bila diberikan sebarang keadaan awal 𝑥(0) dan sebarang keadaan akhir
𝑥(𝑡1 ) akan selalu ada pengontrol 𝑢(𝑡) yang akan mentransfer keadaan awal 𝑥(0) ke
keadaan akhir yang diinginkan 𝑥(𝑡1 ) dalam waktu yang berhingga 𝑡1 .
Penyelesaian sistem tersebut adalah:
𝑡
𝑥(𝑡) = 𝑒 𝑥0 + ∫ 𝑒 𝐴(𝑡−𝜏) 𝐵𝑢(𝜏)𝑑𝜏
𝐴𝑡
0

Bila sistem terkontrol, yaitu ada masukan 𝑢(𝑡) yang mentransfer 𝑥0 ke 𝑥1 dalam
waktu berhingga 𝑡1 . Dalam hal ini 𝑥1 diberikan oleh
𝑡1
𝑥1 = 𝑒 𝐴𝑡1 𝑥0 + ∫ 𝑒 𝐴(𝑡1 −𝜏) 𝐵𝑢(𝜏)𝑑𝜏
0
 Teorema
Syarat perlu dan cukup sistem 𝑥̇ (𝑡) = 𝐴𝑥(𝑡) + 𝐵𝑢(𝑡) 𝑑𝑎𝑛 𝑦(𝑡) = 𝐶𝑥(𝑡) + 𝐷𝑢(𝑡)
terkontrol adalah:
𝑡 𝑇
1. 𝜔(0, 𝑡1 ) = ∫0 1 𝑒 −𝐴𝜏 𝐵𝐵 𝑇 𝑒 −𝐴 𝜏 𝑑𝜏 non singular
2. 𝑀𝑐 = (𝐵 |𝐴𝐵|𝐴2 𝐵|… … …|𝐴𝑛−1 𝐵) mempunyai rank sama dengan n (ukuran
matriks A)

3. Keteramatan
Diberikan Sistem: 𝑥̇ (𝑡) = 𝐴𝑥(𝑡) + 𝐵𝑢(𝑡) 𝑑𝑎𝑛 𝑦(𝑡) = 𝐶𝑥(𝑡) + 𝐷𝑢(𝑡).
Bila setiap keadaan awal 𝑥(0) = 𝑥0 secara tunggal dapat diamati dari setiap
pengukuran keluaran sistem dari waktu 𝑡0 = 0 ke 𝑡 = 𝑡1 , maka sistem dikatakan
"teramati".
Artinya : sebarang keadaan awal 𝑥0 lewat sebarang pengukuran keluaran 𝑦(𝑡)
diamati pada interval waktu 0 ≤ 𝑡 ≤ 𝑡1.

Keluaran sistem 𝑥̇ (𝑡) = 𝐴𝑥(𝑡) + 𝐵𝑢(𝑡) 𝑑𝑎𝑛 𝑦(𝑡) = 𝐶𝑥(𝑡) + 𝐷𝑢(𝑡) diberikan
oleh:
𝑡
𝑦(𝑡) = 𝐶𝑒 𝑥0 + 𝐶 ∫ 𝑒 𝐴(𝑡−𝜏) 𝐵𝑢(𝜏)𝑑𝜏 + 𝐷𝑢(𝑡)
𝐴𝑡
0

Bila diukur keluaran 𝑦(𝑡) pada 𝑡 = 0, maka diperoleh

𝑦(0) = 𝐶𝑥0 + 𝐷𝑢(0)
Bila diukur keluaran 𝑦(𝑡) pada 𝑡𝑠 dengan 0 ≤ 𝑡𝑠 ≤ 𝑡1 , maka diperoleh
𝑡
𝑦(𝑡𝑠 ) = 𝐶𝑒 𝐴𝑡𝑠 𝑥0 + 𝐶 ∫ 𝑒 𝐴(𝑡𝑠 −𝜏) 𝐵𝑢(𝜏)𝑑𝜏 + 𝐷𝑢(𝑡𝑠 )
0

𝑦(𝑡𝑠 ) = 𝐶𝑥(𝑡𝑠 ) + 𝐷𝑢(𝑡𝑠 )

Bila keadaan awal dapat diamati maka keadaan ini akan muncul pada pengukuran
luaran 𝑦(𝑡𝑠 ) yaitu : 𝑦(𝑡𝑠 ) = 𝐶𝑥(0) + 𝐷𝑢(𝑡𝑠 )
Teorema
Syarat perlu dan cukup sistem 𝑥̇ (𝑡) = 𝐴𝑥(𝑡) + 𝐵𝑢(𝑡) 𝑑𝑎𝑛 𝑦(𝑡) = 𝐶𝑥(𝑡) + 𝐷𝑢(𝑡)
teramati adalah:
𝑡 𝑇
1. Matriks : 𝑚(0, 𝑡) = ∫0 𝑒 𝐴 𝜏 𝐶 𝑇 𝐶𝑒 𝐴𝜏 𝑑𝜏 nonsingular
 𝐶
𝐶𝐴
𝐶𝐴2
2. Matriks keteramatan 𝑀0 = .. mempunyai rank sama dengan n

.
(𝐶𝐴(𝑛−1) )
(ukuran matriks A)

2.8 Feedback Controller

Sifat-sifat sistem dapat dipengaruhi oleh suatu kompensator. Kondisi pada
matriks 𝐴 dan 𝐵 akan diberikan sedemikian hingga matriks 𝐴 + 𝐵𝐹 stabil asimtotik bila
matriks F dipilih sesuai dengan apa yang diinginkan. Dimana matriks F disini
merupakan Feedback Controller.
Sistem 𝑥̇ (𝑡) = 𝐴𝑥(𝑡) + 𝐵𝑢(𝑡) bisa distabilkan apabila ada suatu matriks F
ukuran 𝑚 × 𝑛 sedemikian hingga semua nilai karakteristik 𝜆 dari matriks 𝐴 + 𝐵𝐹,
memiliki ℝ(𝜆) < 0.
Teorema
Sistem 𝑥̇ (𝑡) = 𝐴𝑥(𝑡) + 𝐵𝑢(𝑡) terkontrol bila dan hanya bila untuk setiap polynomial
𝑝(𝜆) = 𝜆𝑛 + 𝑝1 𝜆𝑛−1 + 𝑝2 𝜆𝑛−2 + ⋯ + 𝑝𝑛
Dengan koefisien 𝑝𝑖 real, ada matriks 𝐹 ukuran 𝑚 × 𝑛 sedemikian hingga
det[𝜆𝐼 − (𝐴 + 𝐵𝐹)] = 𝑝(𝜆).

2.9 Observer
Untuk membangun suatu sistem real terkadang membutuhkan biaya yang cukup mahal,
ataupun beberapa sistem memiliki kesulitan dalam melakukan pengukuran. Maka dari
itu diperlukan suatu sistem pembantu yang dinamakan pengamat (observer) yang
mempunyai masukan 𝑢(𝑡) dan keluaran 𝑦(𝑡) dari sistem riil dan keluaran 𝑥̂(𝑡) suatu
pendekatan dari keadaan 𝑥(𝑡) dari sistem riil. Suatu pengamat dari sistem 𝑥̇ (𝑡) =
𝐴𝑥(𝑡) + 𝐵𝑢(𝑡), 𝑦(𝑡) = 𝐶𝑥(𝑡) diasumsikan berbentuk:
𝑧̇ (𝑡) = 𝑃𝑧(𝑡) + 𝑄𝑢(𝑡) + 𝐾𝑦(𝑡)
𝑥̂(𝑡) = 𝑆𝑧(𝑡) + 𝑇𝑢(𝑡) + 𝑅𝑦(𝑡)
Vektor 𝑧(𝑡) adalah keadaan dari pengamat. Sedangka matriks 𝑃, 𝑄, 𝐾, 𝑇, 𝑆, dan 𝑅 adalah
matriks –matriks yang dapat ditentukan.
 BAB III
PEMBAHASAN
3.1 Model Sistem Non-Linier Penyebaran Virus Hepatitis B
Model SEICR penyebaran Hepatitis B dengan pemberian vaksin membagi populasi menjadi
lima sub-populasi sehingga model juga akan terdiri dari lima variabel yaitu:

𝑥 = jumlah populasi individu rentan (Susceptible)

𝑒 = jumlah individu yang terinfeksi virus namun belum menularkan (Exposed)

𝑦 = jumlah individu yang terinfeksi akut (Infected)

𝑐 = jumlah individu pembawa Hepatitis B kronis (Carrier)

𝑣 = jumlah individu yang mengalami pemulihan (Recovery)

Sedangkan parameter yang digunakan model SEICR penyebaran Hepatitis B dengan

pemberian vaksin adalah:

𝛽 = koefisien penularan penyakit

𝜇 = laju kematian dan kelahiran

𝜔 = proporsi dari kegagalan imunisasi

𝜎 = tingkat individu exposed menjadi terinfeksi akut

𝛿 = hilangnya tingkat kekebalan

𝑟1 = tingkat individu yang meninggalkan kelas infeksi akut

𝑟2 = tingkat perubahan individu pembawa yang mengalami pemulihan

𝛼 = pembawa infeksi akut dari garis keturunan

𝜀 = proporsi penularan vertical

𝑝 = proporsi individu suspectible yang tervaksinasi

𝑞 = proporsi individu terinfeksi akut menjadi individu pembawa

Skema berikut menjelaskan tentang hubungan antar sub-populasi yang terjadi pada
penyebaran Hepatitis B.

Berikut adalah penjelasan dari skema diatas:

a. Perpindahan populasi recovered menjadi susceptible

Individu yang tergolong populasi kelas susceptible akan terjadi karena tubuh kehilangan
imun saat virus memasuki tubuh. Populasi yang masuk dalam individu susceptible ialah
individu yang memiliki imun lemah dan populasi yang pernah mengalami pemulihan
setelah sebelumnya terserang Hepatitis B. Persamaan dirumuskan sebagai berikut:
𝑑𝑥
= 𝜇𝜔(1 − 𝜀𝑐) − 𝛽𝑥(𝑦 + 𝛼𝑐) − (1 − 𝜔)𝑝𝑥 + 𝛿𝑣 − 𝜇𝑥
𝑑𝑡
dengan 𝛿𝑣 merupakan tingkat perubahan populasi yang mengalami pemulihan
(recovery) menjadi populasi suspectible.

b. Perpindahan populasi susceptible menjadi populasi exposed

Populasi susceptible akan berubah menjadi individu exposed akan terjadi melalui
penularan secara horizontal. Penularan tersebut terjadi dikarenakan kontak langsung
dengan pembawa infeksi virus. Dimana dirumusukan dengan 𝛽(𝑦 + 𝛼𝑐)𝑥 dengan
𝛽 merupakan tingkat kontak dan 𝛼 sebagai pembawa infeksi. Persamaannya
dirumuskan sebagai berikut:
𝑑𝑒
= 𝛽𝑥(𝑦 + 𝛼𝑐) + 𝜇𝜔𝜀𝑐 − 𝜎𝑒 − 𝜇𝑒
𝑑𝑡

c. Perpindahan populasi yang telah terinfeksi menjadi terinfeksi akut

Perpindahan populasi exposed menjadi populasi terinfeksi akut bergantung pada
besarnya 𝜎 yaitu tingkat individu yang tertular menjadi terinfeksi. Populasi yang telah
tertular Hepatitis B dapat masuk ke dalam kelas populasi terinfeksi akut dalam kurun
waktu tertentu. Persamaannya dirumuskan sebagai berikut:
𝑑𝑦
= 𝜎𝑒 − 𝑟1 𝑦 − 𝜇𝑦
𝑑𝑡
 Dengan 𝑟1 merupakan tingkat individu yang meninggalkan kelas infeksi akut.

d. Perpindahan populasi terinfeksi akut menjadi populasi pembawa virus

Adanya populasi pembawa virus terjadi melalui penularan vertikal dari populasi yang
tergolong terinfeksi akut. Sehingga memunculkan populasi pembawa virus (carriers).
Kasus ini terjadi pada penularan dari ibu hamil yang menularkan virus ke bayinya.
Perumusan penularan ini dilambangkan dengan 𝑞𝑟1 𝑦, dimana 𝑞 merupakan tingkat
individu terinfeksi akut menjadi individu pembawa virus dan 𝑟1 merupakan tingkat
individu yang meninggalkan kelas infeksi akut. Persamaan dari perpindahan populasi
ini dirumusukan sebagai berikut:
𝑑𝑐
= 𝑞𝑟1 𝑦 − 𝑟2 𝑐 − 𝜇𝑐
𝑑𝑡

e. Perpindahan populasi pembawa virus ke populasi (recovery).

Pada populasi (recovery) yang termasuk dalam populasi ini diantaranya: bayi yang
berhasil tervaksinasi (𝜇(1 − 𝜔)), populasi dari susceptible yang sukses
tervaksinasi ((1 − 𝜔)𝑝𝑥), individu terinfeksi akut yang mengalami pemulihan
((1 − 𝑞)𝑟1 𝑦), dan individu pembawa virus yang berhasil pulih dari hepatitis B (𝑟2 𝑐).
Persamaannya dirumuskan sebagai berikut:
𝑑𝑣
= (1 − 𝜔)𝑝𝑥 + (1 − 𝑞)𝑟1 𝑦 + 𝑟2 𝑐 − 𝛿𝑣 − 𝜇𝑣 + 𝜇(1 − 𝑤)
𝑑𝑡
Secara singkat, model epidemik persamaan diferensial yang didapatkan sebagai berikut:

𝑑𝑥
𝑆(𝑡) = = 𝜇𝜔(1 − 𝜀𝑐) − 𝛽𝑥(𝑦 + 𝛼𝑐) − (1 − 𝜔)𝑝𝑥 + 𝛿𝑣 − 𝜇𝑥
𝑑𝑡
𝑑𝑒
𝐸(𝑡) = = 𝛽𝑥(𝑦 + 𝛼𝑐) + 𝜇𝜔𝜀𝑐 − 𝜎𝑒 − 𝜇𝑒
𝑑𝑡
𝑑𝑦
𝐼(𝑡) = = 𝜎𝑒 − 𝑟1 𝑦 − 𝜇𝑦
𝑑𝑡
𝑑𝑐
𝐶(𝑡) = = 𝑞𝑟1 𝑦 − 𝑟2 𝑐 − 𝜇𝑐
𝑑𝑡
𝑑𝑣
𝑅(𝑡) = = (1 − 𝜔)𝑝𝑥 + (1 − 𝑞)𝑟1 𝑦 + 𝑟2 𝑐 − 𝛿𝑣 − 𝜇𝑣 + 𝜇(1 − 𝑤)
𝑑𝑡

Jika mengasumsikan bahwa total populasi yaitu 𝑥 + 𝑒 + 𝑦 + 𝑐 + 𝑣 = 1, maka diperoleh

model SEICR penyebaran Hepatitis B dengan pemberian vaksin adalah sebagai berikut:
𝑥̇ = 𝜇𝜔(1 − 𝜀𝑐) − 𝛽𝑥(𝑦 + 𝛼𝑐) − (1 − 𝜔)𝑝𝑥 + 𝛿[1 − (𝑥 + 𝑒 + 𝑦 + 𝑐)] − 𝜇𝑥
𝑒̇ = 𝛽𝑥(𝑦 + 𝛼𝑐) + 𝜇𝜔𝜀𝑐 − 𝜎𝑒 − 𝜇𝑒
𝑦̇ = 𝜎𝑒 − 𝑟1 𝑦 − 𝜇𝑦
𝑐̇ = 𝑞𝑟1 𝑦 − 𝑟2 𝑐 − 𝜇𝑐

Dengan diberikan nilai awal parameter nya adalah sebagai berikut:

 Parameter Nilai

𝛽 0.85

𝜎 6

𝑟1 4

𝑟2 0.025

𝜇 1/70

𝜀 0.8

𝑝 0.1

𝛼 0.5

𝜔 0.1

𝑞 0.1

𝛿 1/22

3.2 Penentuan Titik Setimbang

Suatu model dapat diperoleh titik kesetimbangannya apabila turunan dari tiap variabelnya
sama dengan nol. Kemudian ukuran populasi ditetepakan sebagai 𝑁 = 𝑥 + 𝑒 + 𝑦 + 𝑐 +
𝑣 = 1. Sehingga diperoleh:

𝑥̇ = 𝜇𝜔(1 − 𝜀𝑐) − 𝛽𝑥(𝑦 + 𝛼𝑐) − (1 − 𝜔)𝑝𝑥 + 𝛿[1 − (𝑥 + 𝑒 + 𝑦 + 𝑐)] − 𝜇𝑥 = 0

𝑒̇ = 𝛽𝑥(𝑦 + 𝛼𝑐) + 𝜇𝜔𝜀𝑐 − 𝜎𝑒 − 𝜇𝑒 = 0

𝑦̇ = 𝜎𝑒 − 𝑟1 𝑦 − 𝜇𝑦 = 0

𝑐̇ = 𝑞𝑟1 𝑦 − 𝑟2 𝑐 − 𝜇𝑐 = 0 …(3.1)

Titik kesetimbangan model ada dua, yaitu titik kesetimbangan bebas penyakit dan titik
kesetimbangan endemik. Titik kesetimbangan yang pertama ialah titik kesetimbangan bebas
penyakit (𝑇1 ). Titik kesetimbangan bebas penyakit ialah keadaan dimana tidak terjadi
penyebaran penyakit, sehingga populasi 𝑒 = 0, 𝑦 = 0, dan 𝑐 = 0.

Selanjutnya nilai 𝑥0 diperoleh dari persamaan (3.1) sebagai berikut :

 𝑥̇ = 0

0 = 𝜇𝜔(1 − 𝜀𝑐) − 𝛽𝑥(𝑦 + 𝛼𝑐) − (1 − 𝜔)𝑝𝑥 + 𝛿[1 − (𝑥 + 𝑒 + 𝑦 + 𝑐)] − 𝜇𝑥

𝛽𝑥(𝑦 + 𝛼𝑐) + (1 − 𝜔)𝑝𝑥 + 𝛿𝑥 + 𝜇𝑥 = 𝜇𝜔(1 − 𝜀𝑐) + 𝛿 − 𝛿𝑒 − 𝛿𝑦 − 𝛿𝑐

𝑥(𝛽(𝑦 + 𝛼𝑐) + (1 − 𝜔)𝑝 + 𝛿 + 𝜇) = 𝜇𝜔(1 − 𝜀𝑐) + 𝛿(1 − 𝑒 − 𝑦 − 𝑐)

𝜇𝜔(1 − 𝜀𝑐) + 𝛿(1 − 𝑒 − 𝑦 − 𝑐)

𝑥=
𝛽(𝑦 + 𝛼𝑐) + (1 − 𝜔)𝑝 + 𝛿 + 𝜇

Karena 𝑒 = 0, 𝑦 = 0, dan 𝑐 = 0, maka diperoleh

𝜇𝜔 + 𝛿
𝑥0 = (3.2)
(1 − 𝜔)𝑝 + 𝛿 + 𝜇

Sehingga diperoleh titik kesetimbangan bebas penyakit adalah 𝑇1 = (𝑥0 , 0,0,0) =

𝜇𝜔+𝛿
((1−𝜔)𝑝+𝛿+𝜇 , 0,0,0).

Selanjutnya akan diselidiki titik kesetimbangan endemik (𝑇2 ) yang diperoleh dari
persamaan (3.1). Titik ini digunakan untuk menunjukkan adanya penyebaran penyakit
dalam populasi sedemikian hingga 𝑒 ≠ 0, 𝑦 ≠ 0, 𝑐 ≠ 0.

Pertama akan dicari 𝑒 ∗ sebagai berikut :

𝑦̇ = 0

0 = 𝜎𝑒 − 𝑟1 𝑦 − 𝜇𝑦

0 = 𝜎𝑒 − 𝑦(𝑟1 + 𝜇)

𝜎𝑒 = 𝑦(𝑟1 + 𝜇)

(𝑟1 + 𝜇)
𝑒∗ = 𝑦∗ (3.3)
𝜎

Kemudian, akan dicari nilai 𝑐 ∗

𝑐̇ = 0

0 = 𝑞𝑟1 𝑦 − 𝑟2 𝑐 − 𝜇𝑐

0 = 𝑞𝑟1 𝑦 − 𝑐(𝑟2 + 𝜇)
 𝑐(𝑟2 + 𝜇) = 𝑞𝑟1 𝑦

𝑞𝑟1
𝑐∗ = 𝑦∗ (3.4)
(𝑟2 + 𝜇)

Selanjutnya, akan diselidiki nilai 𝑥 ∗ sebagai berikut

𝑒̇ = 0

0 = 𝛽𝑥(𝑦 + 𝛼𝑐) + 𝜇𝜔𝜀𝑐 − 𝜎𝑒 − 𝜇𝑒

𝛽𝑥(𝑦 + 𝛼𝑐) = 𝑒(𝜎 + 𝜇) − 𝜇𝜔𝜀𝑐

𝑒(𝜎 + 𝜇) − 𝜇𝜔𝜀𝑐
𝑥∗ =
𝛽(𝑦 + 𝛼𝑐)

𝑦(𝑟1 + 𝜇) 𝑞𝑟1 𝑦
(𝜎 + 𝜇) − 𝜇𝜔𝜀
𝜎 (𝑟2 + 𝜇)
𝑥∗ =
𝛼𝑞𝑟1 𝑦
𝛽 (𝑦 + )
(𝑟2 + 𝜇)

𝑦(𝑟1 + 𝜇)(𝜎 + 𝜇)(𝑟2 + 𝜇) − 𝜇𝜔𝜀(𝑞𝑟1 𝑦)𝜎

𝜎(𝑟2 + 𝜇)
𝑥∗ =
(𝑟 + 𝜇) + 𝛼𝑞𝑟1
𝛽𝑦 ( 2 𝑟 + 𝜇 )
2

(𝑟1 + 𝜇)(𝜎 + 𝜇)(𝑟2 + 𝜇) − 𝜇𝜔𝜀𝑞𝑟1 𝜎

𝑥∗ = (3.5)
𝛽𝜎(𝑟2 + 𝜇 + 𝛼𝑞𝑟1 )

Terakhir, akan dicari 𝑦 ∗ dengan menyatakan 𝑥̇ = 0. Kemudian mensubstitusikan

persamaan-persamaan yang telah diperoleh sebelumnya.

𝑥̇ = 0

0 = 𝜇𝜔(1 − 𝜀𝑐) − 𝛽𝑥(𝑦 + 𝛼𝑐) − (1 − 𝜔)𝑝𝑥 + 𝛿[1 − (𝑥 + 𝑒 + 𝑦 + 𝑐)] − 𝜇𝑥

−𝛿𝑦 − 𝛽𝑥𝑦 = (1 − 𝜔)𝑝𝑥 + 𝛿𝑥 + 𝛿𝑒𝛿 + 𝛿𝑐 + 𝜇𝑥 + 𝛽𝑥𝛼𝑐 − 𝜇𝜔(1 − 𝜀𝑐)

𝑦(−𝛿 − 𝛽𝑥) = ((1 − 𝜔)𝑝 + 𝛿 + 𝜇)𝑥 ∗ − 𝜇𝜔(1 − 𝜀𝑐) + 𝛿𝑒 − 𝛿 + 𝛿𝑐 + 𝛽𝑥𝛼𝑐

𝜀𝑞𝑟1 𝑦 𝛿𝑦(𝑟1 + 𝜇) 𝛿𝑞𝑟1 𝑦 𝛽𝑥 ∗ 𝛼𝑞𝑟1 𝑦

𝑦(−𝛿 − 𝛽𝑥) = ((1 − 𝜔)𝑝 + 𝛿 + 𝜇)𝑥 ∗ − 𝜇𝜔 (1 − )+ −𝛿+ +
𝑟2 + 𝜇 𝜎 𝑟2 + 𝜇 𝑟2 + 𝜇

𝜇𝜔𝜀𝑞𝑟1 𝑦 𝛿𝑦(𝑟1 + 𝜇) 𝛿𝑞𝑟1 𝑦 𝛽𝑥 ∗ 𝛼𝑞𝑟1 𝑦

𝑦(−𝛿 − 𝛽𝑥) = ((1 − 𝜔)𝑝 + 𝛿 + 𝜇)𝑥 ∗ − 𝜇𝜔 + + −𝛿+ +
𝑟2 + 𝜇 𝜎 𝑟2 + 𝜇 𝑟2 + 𝜇
 𝜇𝜔𝜀𝑞𝑟1 𝑦 𝛿𝑦(𝑟1 + 𝜇) 𝛿𝑞𝑟1 𝑦 𝛽𝑥 ∗ 𝛼𝑞𝑟1 𝑦
𝑦(−𝛿 − 𝛽𝑥) − − − − = ((1 − 𝜔)𝑝 + 𝛿 + 𝜇)𝑥 ∗ − (𝜇𝜔 + 𝛿)
𝑟2 + 𝜇 𝜎 𝑟2 + 𝜇 𝑟2 + 𝜇

𝜇𝜔𝜀𝑞𝑟1 𝛿(𝑟1 + 𝜇) 𝛿𝑞𝑟1 𝛽𝑥 ∗ 𝛼𝑞𝑟1

𝑦 [(−𝛿 − 𝛽𝑥) − − − − ] = ((1 − 𝜔)𝑝 + 𝛿 + 𝜇)𝑥 ∗ − (𝜇𝜔 + 𝛿)
𝑟2 + 𝜇 𝜎 𝑟2 + 𝜇 𝑟2 + 𝜇

𝛼𝑞𝑟1 𝑟1 + 𝜇 𝑞𝑟1 𝜇𝜔𝜀𝑞𝑟1

𝑦 [−𝛽𝑥 ∗ (1 + )−𝛿( +1+ )− ] = ((1 − 𝜔)𝑝 + 𝛿 + 𝜇)𝑥 ∗ − (𝜇𝜔 + 𝛿)
𝑟2 + 𝜇 𝜎 𝑟2 + 𝜇 𝑟2 + 𝜇

∗
((1 − 𝜔)𝑝 + 𝛿 + 𝜇)𝑥 ∗ − (𝜇𝜔 + 𝛿)
𝑦 = 𝛼𝑞𝑟 𝑟 +𝜇 𝑞𝑟 𝜇𝜔𝜀𝑞𝑟 (3.6)
−𝛽𝑥 ∗ (1 + 𝑟 + 1𝜇 ) − 𝛿 ( 1 𝜎 + 1 + 𝑟 +1 𝜇 ) − 𝑟 + 𝜇1
2 2 2

Bilangan Reproduksi Dasar

Perhitungan bilangan reproduksi dasara (𝑅0 ) dilakukan dengan menggunakan
metode next generation matrix dari sistem persamaan model yang telah terbentuk. Pada
sistem ini yang termasuk kelas terinfeksi ialah individu pada kelas Exposed, terinfeksi
akut, dan carrier. Sehingga diperoleh bahwa :
𝛽𝑥(𝑦 + 𝛼𝑐) −𝜇𝜔𝜀𝑐 + 𝜎𝑒 + 𝜇𝑒
𝜑=[ 0 ] dan 𝜓 = [ −𝜎𝑒 + 𝑟1 𝑦 + 𝜇𝑦 ]
0 −𝑞𝑟1 𝑦 + 𝑟2 𝑐 + 𝜇𝑐
Selanjutnya diselidiki matriks Jacobian dari 𝜑 dan 𝜓 untuk mendapatkan nilai
𝐹 dan 𝑉 sebagai berikut :
𝜕𝜑(1,1) 𝜕𝜑(1,2) 𝜕𝜑(1,3)
𝜕𝑒 𝜕𝑦 𝜕𝑒 0 𝛽𝑥 𝛽𝑥𝛼
𝜕𝜑(2,1) 𝜕𝜑(2,2) 𝜕𝜑(2,3)
𝐹= = [0 0 0 ]
𝜕𝑒 𝜕𝑦 𝜕𝑐
𝜕𝜑(3,1) 𝜕𝜑(3,2) 𝜕𝜑(3,3) 0 0 0
[ 𝜕𝑒 𝜕𝑦 𝜕𝑐 ]
𝜕𝜓(1,1) 𝜕𝜓(1,2) 𝜕𝜓(1,3)
𝜕𝑒 𝜕𝑦 𝜕𝑒 𝜎+𝜇 0 −𝜇𝜔𝜀
𝜕𝜓(2,1) 𝜕𝜓(2,2) 𝜕𝜓(2,3)
𝑉= = [ −𝜎 𝑟1 + 𝜇 0 ]
𝜕𝑒 𝜕𝑦 𝜕𝑐
𝜕𝜓(3,1) 𝜕𝜓(3,2) 𝜕𝜓(3,3) 0 −𝑞𝑟1 𝑟2 + 𝜇
[ 𝜕𝑒 𝜕𝑦 𝜕𝑐 ]
Selanjutnya akan dicari 𝑉 −1 , yaitu
1
𝑉 −1 = 𝑎𝑑𝑗 𝑉
det 𝑉
 (𝑟1 + 𝜇)(𝑟2 + 𝜇) 𝑞𝑟1 𝜇𝜀𝜔 𝜇𝜔𝜀(𝑟1 + 𝜇)
1
= [ 𝜎(𝑟2 + 𝜇) (𝜎 + 𝜇)(𝑟2 + 𝜇) 𝜎𝜇𝜔𝜀 ]
det 𝑉
𝜎𝑞𝑟1 𝑞𝑟1 (𝜎 + 𝜇) (𝜎 + 𝜇)(𝑟1 + 𝜇)
Dengan det(𝑉) = (𝜎 + 𝜇)(𝑟1 + 𝜇)(𝑟2 + 𝜇) − 𝜎𝑞𝑟1 𝜇𝜔𝜀, didapatkan next generation
matriks sebagai berikut:
𝐾 = 𝐹𝑉 −1
0 𝛽𝑥 𝛽𝑥𝛼 (𝑟1 + 𝜇)(𝑟2 + 𝜇) 𝑞𝑟1 𝜇𝜀𝜔 𝜇𝜔𝜀(𝑟1 + 𝜇)
1
= [0 0 0 ][ 𝜎(𝑟2 + 𝜇) (𝜎 + 𝜇)(𝑟2 + 𝜇) 𝜎𝜇𝜔𝜀 ]
det 𝑉
0 0 0 𝜎𝑞𝑟1 𝑞𝑟1 (𝜎 + 𝜇) (𝜎 + 𝜇)(𝑟1 + 𝜇)
𝑎 𝑏 𝑐
= [0 0 0]
0 0 0
Dimana,
𝛽𝑥𝜎(𝑟2 + 𝜇 + 𝛽𝑥𝛼𝜎𝑞𝑟1 )
𝑎=
(𝛼 + 𝜇)(𝑟1 + 𝜇)(𝑟2 + 𝜇) − 𝜎𝑞𝑟1 𝜇𝜔𝜀
𝛽𝑥(𝜎 + 𝜇)(𝑟2 + 𝜇) + 𝛽𝑥𝛼𝜎𝑞𝑟1 (𝜎 + 𝜇)
𝑏=
(𝛼 + 𝜇)(𝑟1 + 𝜇)(𝑟2 + 𝜇) − 𝜎𝑞𝑟1 𝜇𝜔𝜀
𝛽𝑥(𝜎 + 𝜇)(𝑟1 + 𝜇) + 𝛽𝑥𝜎𝜇𝜔𝜀
𝑐=
(𝛼 + 𝜇)(𝑟1 + 𝜇)(𝑟2 + 𝜇) − 𝜎𝑞𝑟1 𝜇𝜔𝜀
Menurut Van den Driesche dan Watmough (2002), bilangan reproduksi dasar
diperoleh dari nilai eigen dominan next generation matrix. Sehingga 𝑅0 diperoleh
dengan:
|𝐾 − 𝜆𝐼| = 0
𝑎 𝑏 𝑐 𝜆 0 0
|[0 0 0] − [0 𝜆 0]| = 0
0 0 0 0 0 𝜆
2 (𝑎
𝜆 − 𝜆) = 0
𝜆1,2 = 0
𝜆3 = 𝑎
Sehingga dari Persamaan (4.7) didapatkan nilai eigen dominan yaitu:
𝛽𝑥𝜎(𝑟2 + 𝜇 + 𝛽𝑥𝛼𝜎𝑞𝑟1 )
𝑅0 =
(𝛼 + 𝜇)(𝑟1 + 𝜇)(𝑟2 + 𝜇) − 𝜎𝑞𝑟1 𝜇𝜔𝜀
3.3 Linierisasi di Titik Setimbang
Linearisasi dilakukan untuk mengubah suatu sistem non-linear menjadi sistem linear.
Pencarian hasil linearisasi dari sistem persamaan diferensial non-linear dapat menggunakan
matriks Jacobian. Sehingga jika dilihat dari model, didapat:

𝑓1 = 𝑥̇ = 𝜇𝜔(1 − 𝜀𝑐) − 𝛽𝑥(𝑦 + 𝛼𝑐) − (1 − 𝜔)𝑝𝑥 + 𝛿[1 − (𝑥 + 𝑒 + 𝑦 + 𝑐)] − 𝜇𝑥

𝑓2 = 𝑒̇ = 𝛽𝑥(𝑦 + 𝛼𝑐) + 𝜇𝜔𝜀𝑐 − 𝜎𝑒 − 𝜇𝑒

𝑓3 = 𝑦̇ = 𝜎𝑒 − 𝑟1 𝑦 − 𝜇𝑦

𝑓4 = 𝑐̇ = 𝑞𝑟1 𝑦 − 𝑟2 𝑐 − 𝜇𝑐

Matriks Jacobiannya adalah:

𝜕𝑓1 𝜕𝑓1 𝜕𝑓1 𝜕𝑓1

𝜕𝑥 𝜕𝑒 𝜕𝑦 𝜕𝑐
𝜕𝑓2 𝜕𝑓2 𝜕𝑓2 𝜕𝑓2
𝜕𝑥 𝜕𝑒 𝜕𝑦 𝜕𝑐
𝐽=
𝜕𝑓3 𝜕𝑓3 𝜕𝑓3 𝜕𝑓3
𝜕𝑥 𝜕𝑒 𝜕𝑦 𝜕𝑐
𝜕𝑓4 𝜕𝑓4 𝜕𝑓4 𝜕𝑓4
[ 𝜕𝑥 𝜕𝑒 𝜕𝑦 𝜕𝑐 ]

𝜕𝑓1 𝜕𝑓2
= −𝛽(𝑦 + 𝛼𝑐) − (1 − 𝜔)𝑝 − 𝛿 − 𝜇 = 𝛽(𝑦 + 𝛼𝑐)
𝜕𝑥 𝜕𝑥

𝜕𝑓1 𝜕𝑓2
= −𝛿 = −𝜎 − 𝜇
𝜕𝑒 𝜕𝑒

𝜕𝑓1 𝜕𝑓2
= −𝛽𝑥 − 𝛿 = 𝛽𝑥
𝜕𝑦 𝜕𝑦

𝜕𝑓1 𝜕𝑓2
= −𝜇𝜔𝜀 − 𝛽𝑥𝛼 − 𝛿 = 𝛽𝑥𝛼 + 𝜇𝜔𝜀
𝜕𝑐 𝜕𝑐

𝜕𝑓3 𝜕𝑓4
=0 =0
𝜕𝑥 𝜕𝑥

𝜕𝑓3 𝜕𝑓4
=𝜎 =0
𝜕𝑒 𝜕𝑒

𝜕𝑓3 𝜕𝑓4
= −𝑟1 − 𝜇 = 𝑞𝑟1
𝜕𝑦 𝜕𝑦
 𝜕𝑓3 𝜕𝑓4
=0 = −𝑟2 − 𝜇
𝜕𝑐 𝜕𝑐

−𝛽(𝑦 + 𝛼𝑐) − (1 − 𝜔)𝑝 − 𝛿 − 𝜇 −𝛿 −𝛽𝑥 − 𝛿 −𝜇𝜔𝜀 − 𝛽𝑥𝛼 − 𝛿

𝛽(𝑦 + 𝛼𝑐) −𝜎 − 𝜇 𝛽𝑥 𝛽𝑥𝛼 + 𝜇𝜔𝜀
𝐽=[ ]
0 𝜎 −𝑟1 − 𝜇 0
0 0 𝑞𝑟1 −𝑟2 − 𝜇

Sistem pelinearan di sekitar titik kesetimbangan 𝑇1 = (𝑥0 , 0, 0, 0)

𝑥̇ −(1 − 𝜔)𝑝 − 𝛿 − 𝜇 −𝛿 −𝛽𝑥0 − 𝛿 −𝜇𝜔𝜀 − 𝛽𝑥0 𝛼 − 𝛿 𝑥

𝑒̇ 0 −𝜎 − 𝜇 𝛽𝑥0 𝛽𝑥0 𝛼 + 𝜇𝜔𝜀 𝑒
[ ]=[ ] [𝑦 ]
𝑦̇ 0 𝜎 −𝑟1 − 𝜇 0
𝑐̇ 0 0 𝑞𝑟1 −𝑟2 − 𝜇 𝑐

Sistem pelinearan di sekitar titik kesetimbangan 𝑇2 = (𝑥 ∗ , 𝑒 ∗ , 𝑦 ∗ , 𝑐 ∗ )

𝑥̇ −𝛽(𝑦 ∗ + 𝛼𝑐 ∗ ) − ((1 − 𝜔)𝑝 + 𝛿) − 𝜇 −𝛿 −(𝛽𝑥 ∗ + 𝛿) −(𝜇𝜔𝜀 + 𝛽𝛼𝑥 ∗ + 𝛿) 𝑥

𝑒̇ 𝛽(𝑦 ∗ + 𝛼𝑐 ∗ ) −𝛿 − 𝜇 𝛽𝑥 ∗ 𝛽𝛼𝑥 ∗ + 𝜇𝜔𝜀 𝑒
[ ]=[ ] [𝑦]
𝑦̇ 0 𝛿 −𝑟1 − 𝜇 0
𝑐̇ 0 0 𝑞𝑟 1 −𝑟 2−𝜇
𝑐

3.4 Solusi Sistem Linier

1. Solusi untuk Sistem pelinearan di sekitar titik kesetimbangan 𝑇1 = (𝑥0 , 0, 0, 0) =
(0.3131, 0, 0, 0)

𝑥̇ −(1 − 𝜔)𝑝 − 𝛿 − 𝜇 −𝛿 −𝛽𝑥0 − 𝛿 −𝜇𝜔𝜖 − 𝛽𝑥0 𝛼 − 𝛿 𝑥

𝑒̇ 0 −𝜎 − 𝜇 𝛽𝑥0 𝛽𝑥0 𝛼 + 𝜇𝜔𝜀 𝑒
[ ]=[ ] [𝑦 ]
𝑦̇ 0 𝜎 −𝑟1 − 𝜇 0
𝑐̇ 0 0 𝑞𝑟1 −𝑟2 − 𝜇 𝑐

Misalkan

−(1 − 𝜔)𝑝 − 𝛿 − 𝜇 −𝛿 −𝛽𝑥0 − 𝛿 −𝜇𝜔𝜖 − 𝛽𝑥0 𝛼 − 𝛿

0 −𝜎 − 𝜇 𝛽𝑥0 𝛽𝑥0 𝛼 + 𝜇𝜔𝜀
𝐴=[ ]
0 𝜎 −𝑟1 − 𝜇 0
0 0 𝑞𝑟1 −𝑟2 − 𝜇

Substitusi nilai-nilainya, maka matriks 𝐴 menjadi :

1153 1 85 1 8 425 1
− − − 𝑥0 − − − 𝑥0 −
7700 22 100 22 7000 1000 22
421 85 425 8
0 − 𝑥 𝑥 +
70 100 0 1000 0 7000
𝐴=
281
0 6 − 0
70
4 185
0 0 −
[ 10 700 ]
 −0.1497 −0.0455 −0.3116 −0.1797
0 −6.0143 0.2661 0.1342
=[ ]
0 6 −4.0143 0
0 0 0.4 −0.0393
Pertama akan dicari nilai eigennya:

|𝜆𝐼 − 𝐴| = 0

𝜆 + 0.1497 0.0455 0.3116 0.1797

0 𝜆 + 6.0143 −0.2661 −0.1342
|𝜆𝐼 − 𝐴| = | |
0 −6 𝜆 + 4.0143 0
0 0 −0.4 𝜆 + 0.0393

0 = 𝜆4 + 10.2176𝜆3 + 24.4479𝜆2 + 3.9982𝜆 + 0.084

0 = (𝜆 + 0.1497)(𝜆 + 6.6104)(𝜆 + 3.4326)(𝜆 + 0.0249)

Sehingga didapat nilai eigennya adalah

𝜆1 = −0.1497

𝜆2 = −6.6104

𝜆3 = −3.4326

𝜆4 = −0.0249

Kemudian akan dicari vaktor eigennya :

[𝜆𝐼 − 𝐴]𝑣⃗ = 0

 Untuk 𝜆1 = −0.1497 ∶
−0.1497 + 0.1497 0.0455 0.3116 0.1797 𝑣1
0 −0.1497 + 6.0143 −0.2661 −0.1342 𝑣2
[ ] [𝑣3 ]
0 −6 −0.1497 + 4.0143 0
0 0 −0.4 −0.1497 + 0.0393 𝑣4
0 0.0455 0.3116 0.1797 𝑣1
0 5.8646 −0.2661 −0.1342 𝑣2
=[ ] [𝑣3 ] = 0
0 −6 3.8646 0
0 0 −0.4 −0.1104 𝑣4
Maka diperoleh bahwa vektor eigen untuk 𝜆1 = −0.1497 adalah
1
0
𝑉1 = [ ] 𝑠
0
0

 Untuk 𝜆2 = −6.6104:
−6.6104 + 0.1497 0.0455 0.3116 0.1797 𝑣1
0 −6.6104 + 6.0143 −0.2661 −0.1342 𝑣2
[ ] [𝑣3 ]
0 −6 −6.6104 + 4.0143 0
0 0 −0.4 −6.6104 + 0.0393 𝑣4
 −6.4607 0.0455 0.3116 0.1797 𝑣1
0 −0.5961 −0.2661 −0.1342 𝑣2
= [ ] [𝑣3 ] = 0
0 −6 −2.5961 0
0 0 −0.4 −6.5711 𝑣4
Maka diperoleh bahwa vektor eigen untuk 𝜆2 = −6.6104 adalah
0.0398
−0.3962
𝑉2 = [ ]𝑠
0.9159
−0.0557

 Untuk 𝜆3 = −3.4326
−3.4326 + 0.1497 0.0455 0.3116 0.1797 𝑣1
0 −3.4326 + 6.0143 −0.2661 −0.1342 𝑣2
[ ] [𝑣3 ]
0 −6 −3.4326 + 4.0143 0
0 0 −0.4 −3.4326 + 0.0393 𝑣4
−3.2829 0.0455 0.3116 0.1797 𝑣1
0 2.5817 −0.2661 −0.1342 𝑣2
=[ ] [𝑣3 ] = 0
0 −6 0.5817 0
0 0 −0.4 −3.3933 𝑣4
Maka diperoleh bahwa vektor eigen untuk 𝜆3 = −3.4326 adalah
−0.0884
−0.0955
𝑉3 = [ ]𝑠
−0.9847
0.1161

 Untuk 𝜆4 = −0.0249
−0.0249 + 0.1497 0.0455 0.3116 0.1797
0 −0.0249 + 6.0143 −0.2661 −0.1342
[ ]
0 −6 −0.0249 + 4.0143 0
0 0 −0.4 −0.0249 + 0.0393
0.1248 0.0455 0.3116 0.1797 𝑣1
0 5.9894 −0.2661 −0.1342 𝑣2
=[ ] [𝑣3 ] = 0
0 −6 3.9894 0
0 0 −0.4 0.0144 𝑣4
Maka diperoleh bahwa vektor eigen untuk 𝜆4 = −0.0249 adalah
−0.8382
0.0131
𝑉4 = [ ]𝑠
0.0197
0.5449
Dengan demikian diperoleh matriks fundamentalnya yaitu:

𝑒 −0.1497 𝑡 0.0398 𝑒 −6.6104 𝑡 −0.0884 𝑒 −3.4326 𝑡 −0.8382 𝑒 −0.0249 𝑡

𝐹=[ 0 −0.3962 𝑒 −6.6104 𝑡 −0.0955 𝑒 −3.4326 𝑡 0.0131 𝑒 −0.0249 𝑡 ]
0 0.9159 𝑒 −6.6104 𝑡 −0.9847 𝑒 −3.4326 𝑡 0.0197 𝑒 −0.0249 𝑡
0 −0.0557 𝑒 −6.6104 𝑡 0.1161 𝑒 −3.4326 𝑡 0.5449 𝑒 −0.0249 𝑡
Akan dicari matriks transisinya.
 1 0.0398 −0.0884 −0.8382
0 −0.3962 −0.0955 0.0131
𝑇 = [𝑉1 𝑉2 𝑉3 𝑉4 ] = [ ]
0 0.9159 −0.9847 0.0197
0 −0.0557 0.1161 0.5449
1 0.0778 0.0835 1.5334
0 −2.0572 0.2045 0.0421
𝑇 −1 =[ ]
0 −1.9095 −0.8214 0.0756
0 0.1966 0.1959 1.8234

−0.1497 0.0004 0.0002 −0.0000

0 −6.6104 −0.0003 0.0002
𝐷 = 𝑇 −1 𝐴 𝑇 = [ ]
0 −0.0006 −3.4326 0.0002
0 0.0008 0.0001 −0.0249

−0.1497 0 0 0
0 −6.6104 0 0
=[ ]
0 0 −3.4326 0
0 0 0 −0.0249

𝑒 −0.1497𝑡 0 0 0
−0.3962𝑡
0 𝑒 0 0
𝑒 𝐷𝑡 =[ ]
0 0 𝑒 −0.9847𝑡 0
0 0 0 𝑒 0.5449𝑡

Misalkan :

𝑎 = 𝑒 −0.1497𝑡 ;𝑏 = 𝑒 −0.3962𝑡 ; 𝑐 = 𝑒 −0.9847𝑡 ; 𝑑 = 𝑒 0.5449𝑡

Maka matriks 𝑒 𝐷𝑡 menjadi :

𝑒 𝐴𝑡 = 𝑇𝑒 𝐷𝑡 𝑇 −1

𝑎 − 0.8868 0.0778𝑎 − 0.0818𝑏 + 0.1687𝑐 + 0.187𝑑 − 0.4177

𝐴𝑡 −0.4786 0.815𝑏 + 0.1823𝑐 + 0.0025𝑑 + 0.7672
𝑒 =[
−0.0491 −1.884𝑏 + 1.8803𝑐 + 0.0038𝑑 + 0.2066
0.6053 0.1145𝑏 + 1.8802𝑐 + 0.0038𝑑 − 0.3679

0.0835𝑎 − 0.0818𝑏 + 0.1688𝑐 + 0.187𝑑 − 0.4224 1.5334𝑎 + 0.0016𝑏 − 0.0067𝑐 + 1.7347𝑑 − 0.2056
−0.081𝑏 + 0.0784𝑐 + 0.0025𝑑 + 0.1616 −0.0166𝑏 − 0.0072𝑐 + 0.0238𝑑 − 1.6627
]
0.1873𝑏 + 0.8088𝑐 + 0.0038𝑑 − 0.9833 0.0385𝑏 − 0.077 = 44𝑐 + 0.0359𝑑 − 0.1705
−0.0113𝑏 + 0.4476𝑐 + 0.106𝑑 − 0.0846 −0.0023𝑏 − 0.0411𝑐 + 0.9935𝑑 + 1.0148

Sehingga solusinya adalah

𝑥1 (𝑡) 𝑎 − 0.8868 0.0778𝑎 − 0.0818𝑏 + 0.1687𝑐 + 0.187𝑑 − 0.4177

𝑥2 (𝑡) −0.4786 0.815𝑏 + 0.1823𝑐 + 0.0025𝑑 + 0.7672
=[
𝑥3 (𝑡) −0.0491 −1.884𝑏 + 1.8803𝑐 + 0.0038𝑑 + 0.2066
[𝑥4 (𝑡)] 0.6053 0.1145𝑏 + 1.8802𝑐 + 0.0038𝑑 − 0.3679
0.0835𝑎 − 0.0818𝑏 + 0.1688𝑐 + 0.187𝑑 − 0.4224 1.5334𝑎 + 0.0016𝑏 − 0.0067𝑐 + 1.7347𝑑 − 0.2056 𝑥1 (0)
−0.081𝑏 + 0.0784𝑐 + 0.0025𝑑 + 0.1616 −0.0166𝑏 − 0.0072𝑐 + 0.0238𝑑 − 1.6627 𝑥 (0)
] 2
0.1873𝑏 + 0.8088𝑐 + 0.0038𝑑 − 0.9833 0.0385𝑏 − 0.077 = 44𝑐 + 0.0359𝑑 − 0.1705 𝑥3 (0)
−0.0113𝑏 + 0.4476𝑐 + 0.106𝑑 − 0.0846 −0.0023𝑏 − 0.0411𝑐 + 0.9935𝑑 + 1.0148 [𝑥4 (0)]

2. Solusi Sistem pelinearan di sekitar titik kesetimbangan 𝑇2 = (𝑥 ∗ , 𝑒 ∗ , 𝑦 ∗ , 𝑐 ∗ ) =

(0.7750, −0.0101, −0.0152, −0.1543)
𝑥̇ −𝛽(𝑦 ∗ + 𝛼𝑐 ∗ ) − ((1 − 𝜔)𝑝 + 𝛿) − 𝜇 −𝛿 −(𝛽𝑥 ∗ + 𝛿) −(𝜇𝜔𝜖 + 𝛽𝛼𝑥 ∗ + 𝛿)
𝑒̇ 𝛽(𝑦 ∗ + 𝛼𝑐 ∗ ) −𝜎 − 𝜇 𝛽𝑥 ∗ 𝛽𝛼𝑥 ∗ + 𝜇𝜔𝜀
[ ]=[ ]
𝑦̇ 0 𝜎 −𝑟1 − 𝜇 0
𝑐̇ 0 0 𝑞𝑟1 −𝑟2 − 𝜇
Maka diperoleh
−𝛽(𝑦 ∗ + 𝛼𝑐 ∗ ) − ((1 − 𝜔)𝑝 + 𝛿) − 𝜇 −𝛿 −(𝛽𝑥 ∗ + 𝛿) −(𝜇𝜔𝜖 + 𝛽𝛼𝑥 ∗ + 𝛿)
𝛽(𝑦 ∗ + 𝛼𝑐 ∗ ) −𝜎 − 𝜇 𝛽𝑥 ∗ 𝛽𝛼𝑥 ∗ + 𝜇𝜔𝜀
𝐴=[ ]
0 𝜎 −𝑟1 − 𝜇 0
0 0 𝑞𝑟1 −𝑟2 − 𝜇
−0.0712 −0.0455 −0.7042 −0.3760
−0.0785 −6.0143 0.6587 0.3305
=[ ]
0 6 −4.0143 0
0 0 0.4 −0.0393
Akan dicari nilai eigen dari A
|𝜆𝐼 − 𝐴| = 0
𝜆 + 0.0712 0.0455 0.7042 0.3760
0.0785 𝜆 + 6.0143 −0.6587 −0.3305
|𝜆𝐼 − 𝐴| = | |
0 6 𝜆 + 4.0143 0
0 0 −0.4 𝜆 + 0.0393
4 3 2
= 𝜆 + 10.1391𝜆 + 21.2984𝜆 + 1.1198𝜆 − 0.0844

= (𝜆 + 7.2047)(𝜆 + 2.8784)(𝜆 − 0.0417)(𝜆 + 0.0977)

Sehingga didapat nilai eigennya adalah:

𝜆1 = −7.2047

𝜆2 = −2.8784

𝜆3 = 0.0417

𝜆4 = −0.0977

Selanjutnya akan diselidiki vektor eigen dari tiap nilai eigennya, yaitu dengan cara

[𝜆𝐼 − 𝐴]𝑣⃗ = 0

 Untuk 𝜆1 = −7.2047, maka

−7.2047 + 0.0712 0.0455 0.7042 0.3760 𝑣1 0
0.0785 −7.2047 + 6.0143 −0.6587 −0.3305 𝑣2 0
[ ] [𝑣 ] = [ ]
0 −6 −7.2047 + 4.0143 0 3 0
0 0 −0.4 −7.2047 + 0.0393 𝑣4 0
 −7.1335 0.0455 0.7042 0.3760 𝑣1 0
0.0785 −1.1904 −0.6587 −0.3305 𝑣2 0
[ ] [𝑣 ] = [ ]
0 −6 −3.1904 0 3 0
0 0 −0.4000 −7.1654 𝑣4 0

Maka diperoleh bahwa vektor eigen untuk 𝜆1 = −7.2047 adalah

0.0812
−0.4674
𝑉1 = [ ]𝑠
0.8790
−0.0491

 Untuk 𝜆2 = −2.8784, maka diperoleh

−2.8784 + 0.0712 0.0455 0.7042 0.3760 𝑣1 0
0.0785 −2.8784 + 6.0143 −0.6587 −0.3305 𝑣2 0
[ ] [𝑣 ] = [ ]
0 −6 −2.8784 + 4.0143 0 3 0
0 0 −0.4 −2.8784 + 0.0393 𝑣4 0
−2.8072 0.0455 0.7042 0.3760 𝑣1 0
0.0785 3.1359 −0.6587 −0.3305 𝑣2 0
[ ] [𝑣 ] = [ ]
0 −6 1.1359 0 3 0
0 0 −0.4000 −2.8391 𝑣4 0
Maka diperoleh bahwa vektor eigen untuk 𝜆2 = −2.8784 adalah
−0.2230
−0.1796
𝑉2 = [ ]𝑠
−0.9488
0.1337

 Untuk 𝜆3 = 0.0417, maka didapat

0.0417 + 0.0712 0.0455 0.7042 0.3760 𝑣1 0
0.0785 0.0417 + 6.0143 −0.6587 −0.3305 𝑣2 0
[ ] [𝑣 ] = [ ]
0 −6 0.0417 + 4.0143 0 3 0
0 0 −0.4 0.0417 + 0.0393 𝑣4 0
0.1129 0.0455 0.7042 0.3760 𝑣1 0
0.0785 6.0560 −0.6587 −0.3305 𝑣2 0
[ ] [𝑣 ] = [ ]
0 −6 4.0560 0 3 0
0 0 −0.4000 0.0810 𝑣4 0
Maka diperoleh bahwa vektor eigen untuk 𝜆3 = 0.0417 adalah
−0.9764
0.0287
𝑉3 = [ ]𝑠
0.0425
0.2100

 Untuk 𝜆4 = −0.0977, maka

−0.0977 + 0.0712 0.0455 0.7042 0.3760 𝑣1 0
0.0785 −0.0977 + 6.0143 −0.6587 −0.3305 𝑣2 0
[ ] [𝑣 ] = [ ]
0 −6 −0.0977 + 4.0143 0 3 0
0 0 −0.4 −0.0977 + 0.0393 𝑣4 0
 −0.0265 0.0455 0.7042 0.3760 𝑣1 0
0.0785 5.9166 −0.6587 −0.3305 𝑣2 0
[ ] [𝑣 ] = [ ]
0 −6 3.9166 0 3 0
0 0 −0.4000 −0.0584 𝑣4 0
Maka diperoleh bahwa vektor eigen untuk 𝜆4 = −0.0977 adalah
0.9951
−0.0093
𝑉4 = [ ]𝑠
−0.0143
0.0979

Dengan demikian diperoleh matriks fundamentalnya yaitu:

0.0812 𝑒 −7.2047 𝑡 −0.2230 𝑒 −2.8784 𝑡 −0.9764 𝑒 0.0417 𝑡 0.9951 𝑒 −0.0977 𝑡

−7.2047 𝑡 −2.8784 𝑡 0.0417 𝑡 −0.0977 𝑡
𝐹 = [−0.4674 𝑒−7.2047 𝑡 −0.1796 𝑒 −2.8784 𝑡 0.0287 𝑒 0.0417 𝑡 −0.0093 𝑒 −0.0977 𝑡 ]
0.8790 𝑒 −0.9488 𝑒 0.0425 𝑒 −0.0143 𝑒
−7.2047 𝑡 −2.8784 𝑡 0.0417 𝑡
−0.0491 𝑒 0.1337 𝑒 0.2100 𝑒 0.0979 𝑒 −0.0977 𝑡

Akan dicari matriks transisinya

0.0812 −0.2230 −0.9764 0.9951
−0.4674 −0.1796 0.0287 −0.0093
𝑇 = [𝑉1 𝑉2 𝑉3 𝑉4 ] = [ ]
0.8790 −0.9488 0.0425 −0.0143
−0.0491 0.1337 0.2100 0.0979
−7.2047 0.0002 −0.0004 0.0003
−0.0000 −2.8784 −0.0002 0.0001
𝐷 = 𝑇 −1 𝐴 𝑇 = [ ]
−0.0007 0.0002 0.0417 −0.0000
−0.0007 0.0003 0.0000 −0.0977
−7.2047 0 0 0
0 −2.8784 0 0
=[ ]
0 0 0.0417 0
0 0 0 −0.0977

𝑒 −7.2047𝑡 0 0 0
0 𝑒 −2.8784𝑡 0 0
𝑒 𝐷𝑡 =[ 0.0417𝑡 ]
0 0 𝑒 0
0 0 0 𝑒 −0.0977𝑡

Misalkan :

𝑎 = 𝑒 −7.2047𝑡 ;𝑏 = 𝑒 −2.8784𝑡 ; 𝑐 = 𝑒 0.0417𝑡 ; 𝑑 = 𝑒 −0.0977𝑡

Maka matriks 𝑒 𝐷𝑡 menjadi :

𝑒 𝐴𝑡 = 𝑇𝑒 𝐷𝑡 𝑇 −1
 𝑒 𝐴𝑡 =
0.0013𝑎 + 0.0089𝑏 + −0.2971𝑐 + 0.6953𝑑 − 1.0485 −0.1263𝑎 + 0.199𝑏 − 0.4371𝑐 + 0.2343𝑑 − 0.3533
0.008𝑎 + 0.0072𝑏 − 0.0087𝑐 − 0.0065𝑑 − 0.2117 0.7272𝑎 + 0.2576𝑏 + 0.0125𝑐_0.0021𝑑 + 0.4609
[
−0.0151𝑎 + 0.038𝑏 − 0.0129𝑐 − 0.0099𝑑 − 0.014 −1.3765𝑎 + 1.3613𝑏 + 0.0185𝑐 − 0.0033𝑑 + 0.0968
−0.0008𝑎 + 0.0053𝑏 − 0.0639𝑐 + 0.0684𝑑 + 0.1324 0.0768𝑎 + 0.1918𝑏 + 0.0918𝑐 + 0.0132𝑑 − 0.9136

0.08252𝑎 + 0.1672𝑏 − 0.429𝑐 + 0.2365𝑑 + 0.0976 0.0057𝑎 − 0.0366𝑏 − 3.123𝑐 + 3.1532𝑑 − 0.7833
−0.1451𝑎 + 0.1347𝑏 + 0.0126𝑐 − 0.0022𝑑 − 0.149 −0.0333𝑎 − 0.0241𝑏 + 0.0917𝑐 + 0.0295𝑑 − 4.1422
]
0.2729𝑎 + 0.7117𝑏 + 0.0186𝑐 − 0.0034𝑑 − 1.0099 0.0626𝑎 − 0.1532𝑏 + 0.1359𝑐 − 0.0453𝑑 − 0.2746
−0.0152𝑎 + 0.1003𝑏 + 0.0922𝑐 + 0.2327𝑑 + 0.0932 −0.0035𝑎 − 0.0215𝑏 + 0.6716𝑐 + 0.3102𝑑 + 1.5905

Sehingga solusinya adalah

𝑥(𝑡) = 𝑒 𝐴𝑡 ∙ 𝑥(0)

𝑥1 (𝑡)
𝑥2 (𝑡)
=
𝑥3 (𝑡)
[𝑥4 (𝑡)]
0.0013𝑎 + 0.0089𝑏 + −0.2971𝑐 + 0.6953𝑑 − 1.0485 −0.1263𝑎 + 0.199𝑏 − 0.4371𝑐 + 0.2343𝑑 − 0.3533
0.008𝑎 + 0.0072𝑏 − 0.0087𝑐 − 0.0065𝑑 − 0.2117 0.7272𝑎 + 0.2576𝑏 + 0.0125𝑐_0.0021𝑑 + 0.4609
[
−0.0151𝑎 + 0.038𝑏 − 0.0129𝑐 − 0.0099𝑑 − 0.014 −1.3765𝑎 + 1.3613𝑏 + 0.0185𝑐 − 0.0033𝑑 + 0.0968
−0.0008𝑎 + 0.0053𝑏 − 0.0639𝑐 + 0.0684𝑑 + 0.1324 0.0768𝑎 + 0.1918𝑏 + 0.0918𝑐 + 0.0132𝑑 − 0.9136

0.08252𝑎 + 0.1672𝑏 − 0.429𝑐 + 0.2365𝑑 + 0.0976 0.0057𝑎 − 0.0366𝑏 − 3.123𝑐 + 3.1532𝑑 − 0.7833
−0.1451𝑎 + 0.1347𝑏 + 0.0126𝑐 − 0.0022𝑑 − 0.149 −0.0333𝑎 − 0.0241𝑏 + 0.0917𝑐 + 0.0295𝑑 − 4.1422
]
0.2729𝑎 + 0.7117𝑏 + 0.0186𝑐 − 0.0034𝑑 − 1.0099 0.0626𝑎 − 0.1532𝑏 + 0.1359𝑐 − 0.0453𝑑 − 0.2746
−0.0152𝑎 + 0.1003𝑏 + 0.0922𝑐 + 0.2327𝑑 + 0.0932 −0.0035𝑎 − 0.0215𝑏 + 0.6716𝑐 + 0.3102𝑑 + 1.5905

𝑥1 (0)
𝑥2 (0)
𝑥3 (0)
[𝑥4 (0)]

3.5 Sifat-Sifat Sistem

1. Kestabilan Sistem
Berdasarkan hasil perhitungan titik kesetimbangan, didapatkan dua titik kesetimbangan
dari sistem. Oleh karena itu, akan dilakukan analisis untuk masing-masing kestabilan.
a) Titik Kesetimbangan Bebas Penyakit
Analisis titik kesetimbangan bebas penyakit akan dilakukan dengan
mensubtitusikan 𝑇1 = (𝑥0 , 0,0,0) ke persamaan 4.9 (Persamaan Jacobiannya).
Sehingga didapatkan :
 −(1 − 𝜔)𝑝 − 𝛿 − 𝜇 −𝛿 −𝛽𝑥0 − 𝛿 −𝜇𝜔𝜀 − 𝛽𝑥0 𝛼 − 𝛿
0 −𝜎 − 𝜇 𝛽𝑥0 𝛽𝑥0 𝛼 + 𝜇𝜔𝜀
𝐽0 = [ ] (6.1)
0 𝜎 −𝑟1 − 𝜇 0
0 0 𝑞𝑟1 −𝑟2 − 𝜇

Selanjutnya akan dicari persamaan karakteristik dari matriks 𝐽0 dengan

menggunakan 𝑡[𝜆𝐼 − 𝐽0 ] = 0 , dimana :

𝜆 0 0 0 −(1 − 𝜔)𝑝 − 𝛿 − 𝜇 −𝛿 −𝛽𝑥0 − 𝛿 −𝜇𝜔𝜀 − 𝛽𝑥0 𝛼 − 𝛿

0 𝜆 0 0 0 −𝜎 − 𝜇 𝛽𝑥0 𝛽𝑥0 𝛼 + 𝜇𝜔𝜀
[𝜆𝐼 − 𝐽0 ] = [
0 0 𝜆 0] − [ 0 𝜎 −𝑟1 − 𝜇 0
]
0 0 0 𝜆 0 0 𝑞𝑟1 −𝑟2 − 𝜇

𝜆 + (1 − 𝜔)𝑝 + 𝛿 + 𝜇 −𝛿 −𝛽𝑥0 − 𝛿 −𝜇𝜔𝜀 − 𝛽𝑥0 𝛼 − 𝛿

0 𝜆 + (𝜎 + 𝜇) 𝛽𝑥0 𝛽𝑥0 𝛼 + 𝜇𝜔𝜀
=[ ]
0 𝜎 𝜆 + (𝑟1 + 𝜇) 0
0 0 𝑞𝑟1 𝜆 + (𝑟2 + 𝜇)

Sehingga diperoleh:

𝜆 + (1 − 𝜔)𝑝 + 𝛿 + 𝜇 −𝛿 −𝛽𝑥0 − 𝛿 −𝜇𝜔𝜀 − 𝛽𝑥0 𝛼 − 𝛿

0 𝜆 + (𝜎 + 𝜇) 𝛽𝑥0 𝛽𝑥0 𝛼 + 𝜇𝜔𝜀
𝑑𝑒𝑡[𝜆𝐼 − 𝐽0 ] = | |
0 𝜎 𝜆 + (𝑟1 + 𝜇) 0
0 0 𝑞𝑟1 𝜆 + (𝑟2 + 𝜇)

= ((1 − 𝜔)𝑝 + 𝛿 + 𝜇 + 𝜆)(𝜆3 + 𝑙1 𝜆2 + 𝑙2 𝜆 + 𝑙1 )

= 𝜆4 + (𝑙0 +𝑙1 )𝜆3 + (𝑙0 𝑙1 + 𝑙2 )𝜆2 + (𝑙0 𝑙2 + 𝑙3 )𝜆 + 𝑙0 𝑙3

=0

Dapat dituliskan juga sebagai berikut :

𝜆4 + (𝑙0 +𝑙1 )𝜆3 + (𝑙0 𝑙1 + 𝑙2 )𝜆2 + (𝑙0 𝑙2 + 𝑙3 )𝜆 + 𝑙0 𝑙3 = 0 (6.2)

Dengan :

𝑙0 = (1 − 𝜔)𝑝 + 𝛿 + 𝜇

𝑙1 = 𝛿 + 𝑟1 + 𝑟2 + 3𝜇

𝑙2 = (𝜎 + 𝜇)(𝑟1 + 𝜇) + (𝜎 + 𝜇)(𝑟2 + 𝜇) + (𝑟1 + 𝜇)(𝑟2 + 𝜇) − 𝛽𝑥0 𝜎

𝑙3 = (𝜎 + 𝜇)(𝑟1 + 𝜇)(𝑟2 + 𝜇) − (𝛽𝜎(𝑟2 + 𝜇)𝑥0 + 𝑞𝑟1 𝜎(𝛽𝛼𝑥0 + 𝜇𝜔𝜀))

 Penentuan kestabilan dari titik kesetimbangan bebas penyakit dapat ditentukan
dengan metode matrik Routh. Matriks Routh berbentuk sebagai berikut :

𝜆4 1 𝑙0 𝑙1 + 𝑙2 𝑙0 𝑙3
𝜆3 𝑙0 + 𝑙1 𝑙0 𝑙2 + 𝑙3 0
2 = 𝑎1 𝑎2 (6.3)
𝜆
1
𝜆 𝑏1
[𝜆0 ] [ 𝑑1 ]

Dimana :

(𝑙0 + 𝑙1 )(𝑙0 𝑙1 + 𝑙2 ) − (𝑙0 𝑙2 + 𝑙3 ) 𝑙1 (𝑙0 2 + 𝑙0 𝑙1 + 𝑙2 ) − 𝑙3

𝑎1 = =
𝑙0 𝑙1 𝑙0 𝑙1

(𝑙0 + 𝑙1 )(𝑙0 𝑙3 )
𝑎2 = = 𝑙0 𝑙3
𝑙0 +𝑙1

𝑙1 (𝑙0 2 + 𝑙0 𝑙1 + 𝑙2 ) − (𝑙0 + 𝑙1 )2 𝑙0 𝑙2
𝑏1 =
𝑙1 (𝑙0 2 + 𝑙0 𝑙1 + 𝑙2 )

𝑑1 = 𝑙0 𝑙3

Berdasarkan syarat dari matriks Routh dimana akan stabil apabila semua nilai pada
kolom pertama dari matriks bernilai positif. Jadi, titik kesetimbangan bebas penyakit
akan bersifat stabil bebas penyakit apabila 𝑙0 + 𝑙1 > 0, 𝑙0 𝑙3 > 0 dan jika
memenuhi 𝑅0 < 1.

Dengan mesubstitusikan nilai awal parameter-parameter yang sudah diketahui,

didapat 𝑅0 = 0.6055 < 1 untuk 𝑇1 = (0.3131, 0, 0, 0). Karena memenuhi 𝑅0 < 1,
maka untuk titik kesetimbangan bebas penyakit bersifat stabil.

b) Titik Kesetimbangan Endemik

Analisis kestabilan titik kesetimbangan endemik dilakukan dengan mensubtitusikan
𝑇2 = (𝑥 ∗ , 𝑒 ∗ , 𝑦 ∗ , 𝑐 ∗ ) ke matriks 𝐽0 . Sehingga dapat diperoleh:

−𝛽(𝑦 ∗ + 𝑎𝑐 ∗ ) − ((1 − 𝜔)𝑝 − 𝛿) − 𝜇 −𝛿 −(𝛽𝑥 ∗ + 𝛿) −(𝜇𝜔𝜀 + 𝛽𝛼𝑥 ∗ + 𝛿)

−𝛽(𝑦 ∗ + 𝑎𝑐 ∗ ) −𝜎 − 𝜇 −𝛽𝑥 ∗ 𝛽𝛼𝑥 ∗ + 𝜇𝜔𝜀
[ ]
0 𝜎 −𝑟1 − 𝜇 0
0 0 𝑞𝑟1 −𝑟2 − 𝜇
(6.4)
 Misalkan :
𝑎1 = −𝛽(𝑦 ∗ + 𝑎𝑐 ∗ ) − (1 − 𝜔)𝑝 − 𝛿 − 𝜇
𝑏1 = −𝛿
𝑐1 = −(𝛽𝑥 ∗ + 𝛿)
𝑑1 = −(𝜇𝜔𝜀 + 𝛽𝛼𝑥 ∗ + 𝛿)
𝑎2 = −𝛽(𝑦 ∗ + 𝑎𝑐 ∗ )
𝑏2 = −𝜎 − 𝜇
𝑐2 = −𝛽𝑥 ∗
𝑑2 = 𝛽𝛼𝑥 ∗ + 𝜇𝜔𝜀
𝑏3 = 𝜎
𝑐3 = −𝑟1 − 𝜇
𝑐4 = 𝑞𝑟1
𝑑4 = −𝑟2 − 𝜇
Maka :

𝑎1 𝑏1 𝑐1 𝑑1
𝑎 𝑏2 𝑐2 𝑑2
[ 2 ]
0 𝑏3 𝑐3 0
0 0 𝑐4 𝑑4

Kemudian akan dicari 𝑃∗ (𝜆) = |𝜆𝐼 − 𝐽∗ |.

𝜆 0 0 0 𝑎1 𝑏1 𝑐1 𝑑1 𝜆 − 𝑎1 −𝑏1 −𝑐1 −𝑑1

0 𝜆 0 0 𝑎 2 𝑏2 𝑐2 𝑑2 −𝑎2 𝜆 − 𝑏2 −𝑐2 −𝑑2
|[0 0 𝜆 0] − [ 0 𝑐3 ]| = |[ 𝜆 − 𝑐3 ]|
𝑏3 0 0 −𝑏3 0
0 0 0 𝜆 0 0 𝑐4 𝑑4 0 0 −𝑐4 𝜆 − 𝑑4

Sehingga diperoleh persamaan karakteristik :

𝑃∗ (𝜆) = 𝜆4 − 𝑏2 𝜆3 − 𝑐3 𝜆3 − 𝑎1 𝜆3 − 𝑑4 𝜆3 + 𝑎1 𝑏2 𝜆2 − 𝑎1 𝑐2 𝜆2 + 𝑎1 𝑑4 𝜆2 + 𝑏2 𝑐3 𝜆2
− 𝑏3 𝑐2 𝜆2 + 𝑏2 𝑑4 𝜆2 + 𝑐2 𝑑4 𝜆2 + 𝑎1 𝑏2 𝑐3 𝑑4 − 𝑎1 𝑏3 𝑐2 𝑑4 + 𝑎1 𝑏3 𝑐3 𝑑2
− 𝑎2 𝑏1 𝑐3 𝑑4 − 𝑎2 𝑏3 𝑐4 𝑑1 − 𝑎1 𝑏2 𝑐3 𝜆 − 𝑎1 𝑏3 𝑐2 𝜆 + 𝑎2 𝑏1 𝑐3 𝜆 − 𝑎2 𝑏3 𝑐1 𝜆
− 𝑎1 𝑏2 𝑑4 𝜆 + 𝑎2 𝑏1 𝑑4 𝜆 − 𝑎1 𝑐3 𝑑4 𝜆 − 𝑏2 𝑐3 𝑑4 𝜆 + 𝑏3 𝑐2 𝑑4 𝜆 − 𝑏3 𝑐4 𝑑2 𝜆
 = 𝜆4 − (𝑏2 + 𝑐3 + 𝑎1 )𝜆3
+ (−𝑎2 𝑏1 + 𝑎1 𝑏2 + 𝑎1 𝑐3 + 𝑎1 𝑑4 + 𝑏2 𝑐3 − 𝑏3 𝑐2
+ 𝑏2 𝑑4 +𝑐3 𝑑4 )𝜆2
+ (−𝑎1 𝑏2 𝑐3 +𝑎1 𝑏3 𝑐2 + 𝑎2 𝑏1 𝑐3 − 𝑎2 𝑏3 𝑐1 + 𝑎1 𝑏2 𝑑4 + 𝑎2 𝑏1 𝑑4
− 𝑎1 𝑐3 𝑑4 − 𝑏2 𝑐3 𝑑4 + 𝑏3 𝑐2 𝑑4 − 𝑏3 𝑐4 𝑑4 )𝜆 + 𝑎1 𝑏2 𝑐3 𝑑4
− 𝑎1 𝑏3 𝑐2 𝑑4 + 𝑎1 𝑏3 𝑐4 𝑑2 − 𝑎2 𝑏1 𝑐3 𝑑4 + 𝑎2 𝑏3 𝑐4 𝑑1 − 𝑎2 𝑏3 𝑐4 𝑑1

Misalkan :

𝑙1 = 𝑏2 + 𝑐3 + 𝑎1

𝑙2 = −𝑎2 𝑏1 + 𝑎1 𝑏2 + 𝑎1 𝑐3 + 𝑎1 𝑑4 + 𝑏2 𝑐3 − 𝑏3 𝑐2 + 𝑏2 𝑑4 +𝑐3 𝑑4

𝑙3 = −𝑎1 𝑏2 𝑐3 +𝑎1 𝑏3 𝑐2 + 𝑎2 𝑏1 𝑐3 − 𝑎2 𝑏3 𝑐1 + 𝑎1 𝑏2 𝑑4 + 𝑎2 𝑏1 𝑑4 − 𝑎1 𝑐3 𝑑4
− 𝑏2 𝑐3 𝑑4 + 𝑏3 𝑐2 𝑑4 − 𝑏3 𝑐4 𝑑4

𝑙4 = 𝑎1 𝑏2 𝑐3 𝑑4 − 𝑎1 𝑏3 𝑐2 𝑑4 + 𝑎1 𝑏3 𝑐4 𝑑2 − 𝑎2 𝑏1 𝑐3 𝑑4 + 𝑎2 𝑏3 𝑐4 𝑑1 − 𝑎2 𝑏3 𝑐4 𝑑1

Maka persamaan karakteristiknya menjadi :

𝜆4 + 𝑙1 𝜆3 + 𝑙2 𝜆2 + 𝑙3 𝜆 + 𝑙4 (6.5)

Titik kesetimbangan endemik akan stabil jika dan hanya jika 𝑅0 > 1 dan juga
apabila nilai eigen dari akar-akar bagian real dari persamaan karakteristiknya
bernilai negatif. Matriks Routh digunakan untuk mengetahui kestabilan tersebut
yaitu :

𝜆4 1 𝑙2 𝑙4
3 𝑙1 𝑙3 0
𝜆
𝜆2 = 𝑚1 𝑚2
𝜆1 𝑛1
[𝜆 ] [ 𝑒1
0 ]

Dengan :

𝑙1 𝑙2 − 𝑙0 𝑙3 𝑙1 𝑙2 − 𝑙3
𝑚1 = =
𝑙1 𝑙1

𝑚2 = 𝑙 4
 𝑙3 (𝑙1 𝑙2 − 𝑙3 ) − 𝑙1 2 𝑙4
𝑛1 =
𝑙1 𝑙2 − 𝑙3

𝑝1 = 𝑙4

Dengan demikian titik kesetimbangan endemik akan stabil apabila 𝑅0 > 1 dan
memenuhi tiga kondisi berikut :

1. 𝑙1 , 𝑙2 , 𝑙3 , 𝑙4 > 0
2. 𝑙1 𝑙2 − 𝑙3 > 0
3. 𝑙3 (𝑙1 𝑙2 − 𝑙3 ) − 𝑙1 2 𝑙4 > 0

Dengan mesubstitusikan nilai awal parameter-parameter yang sudah diketahui,

didapat 𝑅0 = 3.4679 > 1 untuk 𝑇2 = (0.7750, −0.0101, −0.0152, −0.1543).
Karena memenuhi 𝑅0 > 1, maka untuk titik kesetimbangan endemik bersifat stabil.

2. Keterkontrolan
Apabila diperkenalkan suatu pengontrol yaitu pengaruh terapi untuk menghalangi
tingkat penyebaran infeksi virus Hepatitis B. Pengontrol ini berlaku sebagi suatu matriks
input dari model yang dibahas pada tugas kali ini. Misalkan 𝜂 merupakan tingkat
individu yang tidak terinfeksi akibat melakukan terapi, maka bentuk state space-nya
setelah mensubstitusikan nilai parameter awlanya menjadi
a) Untuk model bebas penyakit, dengan 𝜂 = 0.1
−0.1497 −0.0455 −0.3116 −0.1797 0.1
0 −6.0143 0.2661 0.1342 0.1
𝑋̇ = [ ] 𝑋 + [ ] 𝑢(𝑡)
0 6 −4.0143 0 0
0 0 0.4 −0.0393 0
Akan diselidiki apakah sistem yang terbentuk merupakan sistem yang terkontrol
ataupun tidak. Suatu sistem disebut terkontrol apabila dapat dibentuk matriks
keterkontrolannya, 𝑀𝐶 , dimana 𝑟𝑎𝑛𝑘(𝑀𝐶 ) = 𝑛 = dim(𝐴). Karena dim(𝐴) = 4,
maka matriks keterkontrolan dari sistem yang diberikan adalah
𝑀𝐶 = [𝐵 𝐴𝐵 𝐴𝐵 2 𝐴𝐵 3 ]
0.1 −0.0195 −0.1567 1.6834
0.1 −0.6014 3.7768 −24.2840
=[ ]
0 0.6 −6.0172 46.8157
0 0 0.2400 −2.4163
 Jelas bahwa 𝑟𝑎𝑛𝑘(𝑀𝐶 ) = 4 = dim(𝐴), hal ini dibuktikan dengan tiap vektor
kolomnya saling bebas linear, sehingga dapat disimpulkan bahwa sistemnya
terkontrol.
b) Untuk model endemik penyakit, misalkan 𝜂 = 0.1 maka diperoleh
−0.0712 −0.0455 −0.7042 −0.3760 0.1
−0.0785 −6.0143 0.6587 0.3305 0.1
𝑋̇ = [ ] 𝑋 + [ ] 𝑢(𝑡)
0 6 −4.0143 0 0
0 0 0.4 −0.0393 0
Akan diselidiki sifat keterkontrolan dari sistem yang diberikan, akan dibentuk
matriks keterkontrolannya yaitu
0.1 −0.0117 −0.3940 4.0235
0.1 −0.6093 4.0605 −28.3055
𝑀𝐶 = [ ]
0 0.6 −6.0643 48.7069
0 0 0.24 −2.4351
Karena vektor-vektor kolom dari 𝑀𝐶 saling bebas linear, maka jelas bahwa
𝑟𝑎𝑛𝑘(𝑀𝐶 ) = 4, sehingga dapat dikatakan bahwa sistemnya terkontrol.

3. Keteramatan
Suatu sistem dikatakan teramati apabila terdapat matriks keteramatan 𝑀𝑂 =
[𝐶 𝐶𝐴 𝐶𝐴2 𝐶𝐴3 … 𝐶𝐴𝑛−1 ]𝑇 dan mempunyai 𝑟𝑎𝑛𝑘(𝑀𝑜 ) = dim(𝐴) = 𝑛.
Karena didapatkan dua titik kesetimbangan dari sistem, maka akan dibentuk matriks
keteramatan 𝑀𝑜 untuk masing-masing titik kesetimbangan
a) Titik Kesetimbangan Bebas Penyakit
Akan dibentuk matriks keteramatan 𝑀𝑜 untuk titik kesetimbangan bebas penyakit
𝜇𝜔+𝛿
yaitu 𝑇1 = (𝑥0 , 0,0,0) dimana 𝑥0 = (1−𝜔)𝑝+𝛿+𝜇. Diketahui bahwa:

−(1 − 𝜔)𝑝 − 𝛿 − 𝜇 −𝛿 −𝛽𝑥0 − 𝛿 −𝜇𝜔𝜀 − 𝛽𝑥0 𝛼 − 𝛿

0 −𝜎 − 𝜇 𝛽𝑥0 𝛽𝑥0 𝛼 + 𝜇𝜔𝜀
𝐴=[ ]
0 𝜎 −𝑟1 − 𝜇 0
0 0 𝑞𝑟1 −𝑟2 − 𝜇

dengan dim(𝐴) = 𝑛 = 4 dan 𝐶 = [1 0 0 0]

Dengan mensubstitusikan nilai awal parameter-parameter yang sudah diketahui

sebelumnya, didapat 𝑇1 = (0.3131, 0, 0, 0) dan
 −0.1497 −0.0455 −0.3116 −0.1797
0 −6.0143 0.2661 0.1342
𝐴=[ ]
0 6 −4.0143 0
0 0 0.4 −0.0393

Sehingga dapat dibentuk matriks keteramatan

𝐶 1 0 0 0
𝐶𝐴 −0.1497 −0.0455 −0.3116 −0.1797
𝑀𝑜 = [ 2 ] = [ ]
𝐶𝐴 0.0224 −1.5893 1.2135 0.0279
3 −0.0034 16.8386 −5.2901 −0.2184
𝐶𝐴

Terlihat jelas bahwa 𝑟𝑎𝑛𝑘(𝑀𝑜 ) = 4 = dim(𝐴), maka dapat disimpulkan bahwa

untuk titik kesetimbangan bebas penyakit, sistem teramati.

b) Titik Kesetimbangan Endemik

Akan dibentuk matriks keteramatan 𝑀𝑜 untuk titik kesetimbangan endemik yaitu
𝑇2 = (𝑥 ∗ , 𝑒 ∗ , 𝑦 ∗ , 𝑐 ∗ ) dimana
(𝑟1 + 𝜇)(𝜎 + 𝜇)(𝑟2 + 𝜇) − 𝜇𝜔𝜀𝑞𝑟1 𝜎
𝑥∗ =
𝛽𝜎(𝑟2 + 𝜇 + 𝛼𝑞𝑟1 )

∗
((1 − 𝜔)𝑝 + 𝛿 + 𝜇)𝑥 ∗ − (𝜇𝜔 + 𝛿)
𝑦 = 𝛼𝑞𝑟 𝑟 +𝜇 𝑞𝑟 𝜇𝜔𝜀𝑞𝑟
−𝛽𝑥 ∗ (1 + 𝑟 + 1𝜇 ) − 𝛿 ( 1 𝜎 + 1 + 𝑟 +1 𝜇 ) − 𝑟 + 𝜇1
2 2 2

(𝑟1 + 𝜇)
𝑒∗ = 𝑦∗
𝜎

𝑞𝑟1
𝑐∗ = 𝑦∗
(𝑟2 + 𝜇)
Diketahui bahwa:
−𝛽(𝑦∗ + 𝛼𝑐∗ ) − ((1 − 𝜔)𝑝 + 𝛿) − 𝜇 −𝛿 −(𝛽𝑥∗ + 𝛿) −(𝜇𝜔𝜀 + 𝛽𝛼𝑥∗ + 𝛿)
𝛽(𝑦∗ + 𝛼𝑐∗ ) −𝜎 − 𝜇 𝛽𝑥∗ 𝛽𝛼𝑥∗ + 𝜇𝜔𝜀
𝐴=[ ]
0 𝜎 −𝑟1 − 𝜇 0
0 0 𝑞𝑟1 −𝑟2 − 𝜇
dengan dim(𝐴) = 𝑛 = 4 dan 𝐶 = [1 0 0 0]

Dengan mensubstitusikan nilai awal parameter-parameter yang sudah diketahui

sebelumnya, didapat 𝑇2 = (0.7750, −0.0101, −0.0152, −0.1543) dan
 −0.0712 −0.0455 −0.7042 −0.3760
−0.0785 −6.0143 0.6587 0.3305
𝐴=[ ]
0 6 −4.0143 0
0 0 0.4 −0.0393

Sehingga dapat dibentuk matriks keteramatan

𝐶 1 0 0 0
𝐶𝐴 −0.0712 −0.0455 −0.7042 −0.3760
𝑀𝑜 = [ 2 ] = [ ]
𝐶𝐴 0.0086 −3.9486 2.6967 −0.0265
3 0.3093 39.9280 −13.4220 −1.3094
𝐶𝐴

Terlihat jelas bahwa 𝑟𝑎𝑛𝑘(𝑀𝑜 ) = 4 = dim(𝐴), maka dapat disimpulkan bahwa

untuk titik kesetimbangan endemik, sistem teramati.

3.6 Feedback Controller

Akan dikonstruksikan suatu kontrol umpan balik (feedback controller) dimana umpan
balik keadaannya adalah 𝑢(𝑡) = 𝐹𝑥(𝑡) sedemikian hingga sistem akan menjadi
𝑥̇ (𝑡) = (𝐴 + 𝐵𝐹)𝑥(𝑡)

Untuk mendapatkan sistem 𝑥̇ (𝑡) = (𝐴 + 𝐵𝐹)𝑥(𝑡) yang stabil dan terkontrol, maka akan
dicari matriks F sedemikian hingga semua nilai karakteristik 𝜆 dari matriks 𝐴 + 𝐵𝐹 adalah
harus 𝑅𝑒(𝜆) < 0 dan det[𝜆𝐼 − (𝐴 + 𝐵𝐹)] = 𝑝(𝜆) dimana 𝑝(𝜆) = 𝜆𝑛 + 𝑝1 𝜆𝑛−1 +
𝑝2 𝜆𝑛−2 + ⋯ + 𝑝𝑛 .

Dari perhitungan sebelumnya telah diperoleh model linear sistem bebas penyakit, yaitu

−0.1497 −0.0455 −0.3118 −0.1797 0.1

0 −6.0143 0.2663 0.1343 0.1
𝑥̇ (𝑡) = [ ] 𝑥(𝑡) + [ ] 𝑢(𝑡)
0 6 −4.0143 0 0
0 0 0.4 −0.0393 0
Misalkan matriks 𝐹 = [𝑓1 𝑓2 𝑓3 𝑓4 ], maka

−0.1497 −0.0455 −0.3118 −0.1797 0.1

0 −6.0143 0.2663 0.1343 0.1
𝐴 + 𝐵𝐹 = [ ] + [ ] [𝑓1 𝑓2 𝑓3 𝑓4 ]
0 6 −4.0143 0 0
0 0 0.4 −0.0393 0
0.1𝑓1 − 0.1497 𝑓2 − 0.0455 0.1𝑓3 − 0.3118 0.1𝑓4 − 0.1797
0.1𝑓1 0.1𝑓2 − 6.0143 0.1𝑓3 + 0.2663 0.1𝑓4 + 0.1343
=[ ]
0 6 −4.0143 0
0 0 0.4 −0.0393
Sedemikian hingga,
 det(𝜆𝐼 − (𝐴 + 𝐵𝐹))
𝜆 − 0.1𝑓1 + 0.1497 −𝑓2 + 0.0455 −0.1𝑓3 + 0.3118 −0.1𝑓4 + 0.1797
−0.1𝑓1 𝜆 − 0.1𝑓2 + 6.0143 −0.1𝑓3 − 0.2663 −0.1𝑓4 − 0.1343
=| |
0 −6 𝜆 + 4.0143 0
0 0 −0.4 𝜆 + 0.0393
det(𝜆𝐼 − (𝐴 + 𝐵𝐹))
= 𝜆4 + (10.2176 − 0.1𝑓1 − 0.1𝑓2 − 0.6𝑓3 )𝜆3 + (24.4467
− 1.0022𝑓1 − 0.4203𝑓2 ) 𝜆2
+ (3.9978 − 2.0884𝑓1 − 0.0765𝑓2 − 0.1134𝑓3 − 0.24𝑓4 )𝜆 + 0.0844
− 0.0024𝑓2 − 0.0035𝑓3 − 0.0359𝑓4 − 0.0052𝑓1

Jika diinginkan nilai eigen 𝜆 nya adalah -1, -2, -3, dan -4, maka didapat:

𝑝(𝜆) = (𝜆 + 1)(𝜆 + 2)(𝜆 + 3)(𝜆 + 4) = 𝜆4 + 10 𝜆3 + 35 𝜆2 + 50𝜆 + 24

Karena det(𝜆𝐼 − (𝐴 + 𝐵𝐹)) = 𝑝(𝜆), dengan menyamakan koefisiennya didapat:

10.2176 − 0.1𝑓1 − 0.1𝑓2 − 0.6𝑓3 = 10 ⟹ 0.1𝑓1 + 0.1𝑓2 + 0.6𝑓3 = 0.2176 (∗)

24.4467 − 1.0022𝑓1 − 0.4203𝑓2 = 35 ⟹ −1.0022𝑓1 − 0.4203𝑓2 = 10.5533 (∗∗)

3.9978 − 2.0884𝑓1 − 0.0765𝑓2 − 0.1134𝑓3 − 0.24𝑓4 = 50

⟹ −2.0884𝑓1 − 0.0765𝑓2 − 0.1134𝑓3 − 0.24𝑓4 = 46.0022 (∗∗∗)

0.0844 − 0.0024𝑓2 − 0.0035𝑓3 − 0.0359𝑓4 − 0.0052𝑓1 = 24

⟹ −0.0052𝑓1 − 0.0024𝑓2 − 0.0035𝑓3 − 0.0359𝑓4
= 23.9156 (∗∗∗∗)

Dari persamaan (∗), (∗∗), (∗∗∗), (∗∗∗∗), didapat:

𝑓1 = 59.5715
𝑓2 = −167.1566
𝑓3 = 18.2935
𝑓4 = −665.4101

Maka matriks F adalah 𝐹 = [59.5715 −167.1566 18.2935 −665.4101]

Sehingga sistemnya menjadi:
𝑥̇ (𝑡) = (𝐴 + 𝐵𝐹)𝑥(𝑡)

5.8075 −16.7612 1.5176 −66.7207

5.9572 −22.7300 2.0957 −66.4067
𝑥̇ (𝑡) = [ ] 𝑥(𝑡)
0 6 −4.0143 0
0 0 0.4 −0.0393

Untuk sistem endemik penyakit, diketahui bahwa

 −0.0712 −0.0455 −0.7042 −0.3760 0.1
−0.0785 −6.0143 0.6587 0.3305 0.1
𝑥̇ (𝑡) = [ ] 𝑥(𝑡) + [ ] 𝑢(𝑡)
0 6 −4.0143 0 0
0 0 0.4 −0.0393 0
Maka misalkan 𝐹 = [𝑓1 𝑓2 𝑓3 𝑓4 ], sehingga
−0.0712 −0.0455 −0.7042 −0.3760 0.1
−0.0785 −6.0143 0.6587 0.3305 0.1
𝐴 + 𝐵𝐹 = [ ] + [ ] [𝑓1 𝑓2 𝑓3 𝑓4 ]
0 6 −4.0143 0 0
0 0 0.4 −0.0393 0
0.1𝑓1 − 0.0712 0.1𝑓2 − 0.0455 0.1𝑓3 − 0.7042 0.1𝑓4 − 0.3760
0.1𝑓1 − 0.0785 0.1𝑓2 − 6.0143 0.1𝑓3 + 0.6587 0.1𝑓4 + 0.3305
=[ ]
0 6 −4.0143 0
0 0 0.4 −0.0393
Akan dicari persamaan karakteristik dari 𝐴 + 𝐵𝐹, yaitu dengan cara
det(𝜆𝐼 − (𝐴 + 𝐵𝐹))
𝜆 − 0.1𝑓1 + 0.0712 −0.1𝑓2 + 0.0455 −0.1𝑓3 + 0.7042 −0.1𝑓4 + 0.3760
−0.1𝑓1 + 0.0785 𝜆 − 0.1𝑓2 + 6.0143 −0.1𝑓3 − 0.6587 −0.1𝑓4 − 0.3305
=| |
0 −6 𝜆 + 4.0143 0
0 0 −0.4 𝜆 + 0.0393
= 𝜆4 + (−0.1𝑓1 − 0.1𝑓2 + 10.1391)𝜆3 + (−1.0022𝑓1 − 0.4046𝑓2 −
0.6𝑓3 + 21.2984)𝜆2 + (−1.6175𝑓1 − 0.0128𝑓2 − 0.0192𝑓3 −
0.24𝑓4 + 1.1198)𝜆 + (0.1075𝑓1 + 0.00011517𝑓2 + 0.00017213𝑓3
+0.0018𝑓4 − 0.0844)

Jika diinginkan nilai eigen 𝜆 nya adalah -1, -2, -3, dan -4, maka didapat:

𝑝(𝜆) = (𝜆 + 1)(𝜆 + 2)(𝜆 + 3)(𝜆 + 4) = 𝜆4 + 10 𝜆3 + 35 𝜆2 + 50𝜆 + 24

Karena det(𝜆𝐼 − (𝐴 + 𝐵𝐹)) = 𝑝(𝜆), dengan menyamakan koefisiennya didapat:

−0.1𝑓1 − 0.1𝑓2 + 10.1391 = 10

−0.1𝑓1 − 0.1𝑓2 = −0.1391 (𝑖)
−1.0022𝑓1 − 0.4046𝑓2 − 0.6𝑓3 + 21.2984 = 35
−1.0022𝑓1 − 0.4046𝑓2 − 0.6𝑓3 = 13.7016 (𝑖𝑖)
−1.6175𝑓1 − 0.0128𝑓2 − 0.0192𝑓3 − 0.24𝑓4 + 1.1198 = 50
−1.6175𝑓1 − 0.0128𝑓2 − 0.0192𝑓3 − 0.24𝑓4 = 48.8802 (𝑖𝑖𝑖)
0.1075𝑓1 + 0.00011517𝑓2 + 0.00017213𝑓3 + 0.0018𝑓4 − 0.0844 = 24
0.1075𝑓1 + 0.00011517𝑓2 + 0.00017213𝑓3 + 0.0018𝑓4 = 24.0844 (𝑖𝑣)

Dari persamaan (𝑖), (𝑖𝑖), (𝑖𝑖𝑖), dan (𝑖𝑣) diperoleh penyelesaiannya yaitu
 𝑓1 = 256.5
𝑓2 = −255.1
𝑓3 = −279.3
𝑓4 = −1896.5
Maka matriks F adalah 𝐹 = [256.5 −255.1 −279.3 −1896.5]
Sehingga sistemnya menjadi:
𝑥̇ (𝑡) = (𝐴 + 𝐵𝐹)𝑥(𝑡)

25.5788 −25.5555 −28.6342 −190.0260

25.5715 −31.5243 −27.2713 −189.3195
𝑥̇ (𝑡) = [ ] 𝑥(𝑡)
0 6 −4.0143 0
0 0 0.4 −0.0393

3.7 Observer
Akan dikonstruksikan suatu pengamat (observer) yang memiliki bentuk
𝑥̂̇(𝑡) = 𝐴𝑥̂(𝑡) + 𝐵𝑢(𝑡) + 𝐾(𝑦(𝑡) − 𝑦̂(𝑡))
dimana 𝑦̂(𝑡) = 𝐶𝑥̂(𝑡), dengan kata lain akan dicari matriks K, sedemikian hingga agar error
dari observernya mendekati atau kovergen ke nol, maka matriks (𝐴 − 𝐾𝐶) harus stabil
asimtotik (𝜆 < 0).
Dari perhitungan sebelumnya telah diperoleh model linear sistem bebas penyakit, yaitu
 0.1497  0.0455  0.3118  0.1797
 0  6.0143 0.2663 0.1343 

X  X
 0 6  4.0143 0 
 
 0 0 0.4  0.0393

𝑦(𝑡) = [1 0 0 0]𝑋
Misalkan 𝐾 = [𝑘1 𝑘2 𝑘3 𝑘4 ]𝑇 , maka diperoleh bahwa
−0.1497 −0.0455 −0.3118 −0.1797 𝑘1
0 −6.0143 0.2663 0.1343 𝑘
𝐴 − 𝐾𝐶 = [ ] − [ 2 ] [1 0 0 0]
0 6 −4.0143 0 𝑘3
0 0 0.4 −0.0393 𝑘4
−0.1497 − 𝑘1 −0.0455 −0.3118 −0.1797
−𝑘2 −6.0143 0.2663 0.1343
=[ ]
−𝑘3 6 −4.0143 0
−𝑘4 0 0.4 −0.0393
Sedemikian hingga,
 𝜆 + 𝑘1 + 0.1497 0.0455 0.3118 0.1797
𝑘2 𝜆 + 6.0143 −0.2663 −0.1343
det(𝜆𝐼 − (𝐴 − 𝐾𝐶) = | |
𝑘3 −6 𝜆 + 4.0143 0
𝑘4 0 −0.4 𝜆 + 0.0393
𝑑𝑒𝑡(𝜆𝐼 − (𝐴 − 𝐾𝐶) = 𝜆4 + (𝑘1 + 10.217607402597402)𝜆3 + (10.067857142857143𝑘1
−0.045454545454545𝑘2 − 0.311759545454545𝑘3
−0.179749902597403𝑘4 + 24.446774359925790)𝜆2
+(22.939210816326533𝑘1 − 2.054810519480520𝑘2
−1.971263410714286𝑘3 − 1.808739071892393𝑘4
+3.998320041121123)𝜆 + (0.563396655247813𝑘1
−0.512054312152133𝑘2 − 0.509005332606679𝑘3
+4.328214816074741𝑘4 + 0.084363161493601)
Misalkan diinginkan pole-pole dari (𝐴 − 𝐾𝐶) adalah −0.5, −1, −1 + 𝑖, dan −1 − 𝑖, hal ini
berarti bahwa
𝑑𝑒𝑡(𝜆𝐼 − (𝐴 − 𝐾𝐶)) = (𝜆 + 0.5)(𝜆 + 1)(𝜆 + 1 − 𝑖)(𝜆 + 1 + 𝑖)
= 𝜆4 + 3.5𝜆3 + 5.5𝜆2 + 4𝜆 + 1
Dengan demikian diperoleh bahwa,
𝑘1 + 10.217597402597402 = 3.5 ⟹ 𝑘1 = −6.717597402597402
10.067857142857143𝑘1 − 0.045454545454545𝑘2 − 0.311759545454545𝑘3
− 0.179749902597403𝑘4 + 24.446774359925790 = 5.5
−0.045454545454545𝑘2 − 0.311759545454545𝑘3 − 0.179749902597403𝑘4
= 48.685036632653059 (∗)
22.939210816326533𝑘1 − 2.054810519480520𝑘2 − 1.971263410714286𝑘3
− 1.808739071892393𝑘4 + 3.998320041121123 = 4
− 2.054810519480520𝑘2 − 1.971263410714286𝑘3 − 1.808739071892393𝑘4
= 154.0980629562682 (∗∗)
0.563396655247813𝑘1 − 0.512054312152133𝑘2 − 0.509005332606679𝑘3
+ 4.328214816074741𝑘4 + 0.084363161493601 = 1
−0.512054312152133𝑘2 − 0.509005332606679𝑘3 + 4.328214816074741𝑘4
= 4.700308746431172 (∗∗∗)
Dari persamaan (∗), (∗∗), dan (∗∗∗), dengan bantuan MATLAB diperoleh penyelesaian
𝑘2 = 0.899008546742930
𝑘3 = −1.648503686627838
𝑘4 = −0.076648893580230
Jadi diperoleh bentuk observer dari model linear sistem bebas penyakit adalah
−0.1497 −0.0455 −0.3118 −0.1797 0.1
0 −6.0143 0.2663 0.1343 0.1
𝑥̂̇(𝑡) = [ ] 𝑥̂(𝑡) + [ ] 𝑢(𝑡)
0 6 −4.0143 0 0
0 0 0.4 −0.0393 0
−6.717597402597402
0.899008546742930
+[ ] (𝑦(𝑡) − [1 0 0 0]𝑥̂(𝑡))
−1.648503686627838
−0.076648893580230

Sedangkan untuk model linear sistem endemik penyakit, diperoleh bahwa sistemnya
berbentuk
−0.0712 −0.0455 −0.7042 −0.3760
𝑋̇ = [ −0.0785 −6.0143 0.6587 0.3305 ] 𝑋
0 6 −4.0143 0
0 0 0.4 −0.0393
𝑦(𝑡) = [1 0 0 0]𝑋
Akan dikonstruksikan matriks 𝐾, misalkan 𝐾 = [𝑘1 𝑘2 𝑘3 𝑘4 ]𝑇 , maka diperoleh
bahwa
−0.0712 − 𝑘1 −0.0455 −0.7042 −0.3760
−0.0785 − 𝑘2 −6.0143 0.6587 0.3305
𝐴 − 𝐾𝐶 = [ ]
−𝑘3 6 −4.0143 0
−𝑘4 0 0.4 −0.0393
Sehingga,
𝜆 + 0.0712 + 𝑘1 0.0455 0.7042 0.3760
0.0785 + 𝑘2 𝜆 + 6.0143 −0.6587 −0.3305
det(𝜆𝐼 − (𝐴 − 𝐾𝐶)) = | |
𝑘3 −6 𝜆 + 4.0143 0
𝑘4 0 −0.4 𝜆 + 0.0393
det(𝜆𝐼 − (𝐴 − 𝐾𝐶)) = 𝜆4 + (𝑘1 + 10.139118333937578)𝜆3 + (10.067857142857143𝑘1
−0.045454545454545𝑘2 − 0.704175126224167𝑘3
−0.375957692982213𝑘4 + 21.298597592942841)𝜆2
+(20.584717331708806𝑘1 − 4.409304004098246𝑘2
−4.443099346594617𝑘3 − 3.785341448535219𝑘4
+1.120853403431979)𝜆 + (1.368124920123941 × 10−16 𝑘1
−1.075450967399947𝑘2 − 1.078011564941375𝑘3
−9.047561847860013𝑘4 − 0.084400390310855)
Misalkan diinginkan pole-pole dari (𝐴 − 𝐾𝐶) adalah −0.5, −1, −1 + 𝑖, dan −1 − 𝑖, hal ini
berarti bahwa
𝑑𝑒𝑡(𝜆𝐼 − (𝐴 − 𝐾𝐶) = (𝜆 + 0.5)(𝜆 + 1)(𝜆 + 1 − 𝑖)(𝜆 + 1 + 𝑖)
= 𝜆4 + 3.5𝜆3 + 5.5𝜆2 + 4𝜆 + 1
 Sehingga diperoleh bahwa,
𝑘1 + 10.139118333937578 = 3.5 ⟹ 𝑘1 = −6.639118333937578
10.067857142857143𝑘1 − 0.045454545454545𝑘2 − 0.704175126224167𝑘3
− 0.375957692982213𝑘4 + 21.298597592942841 = 5.5
−0.045454545454545𝑘2 − 0.704175126224167𝑘3 − 0.375957692982213𝑘4
= 51.043097347664414 (𝑖)
20.584717331708806𝑘1 − 4.409304004098246𝑘2 − 4.443099346594617𝑘3
− 3.785341448535219𝑘4 + 1.120853403431979 = 4
−4.409304004098246𝑘2 − 4.443099346594617𝑘3 − 3.785341448535219𝑘4
= 139.5435208324386 (𝑖𝑖)
1.368124920123941 × 10−16 𝑘1 − 1.075450967399947𝑘2 − 1.078011564941375𝑘3
− 9.047561847860013𝑘4 − 0.084400390310855 = 1
−1.075450967399947𝑘2 − 1.078011564941375𝑘3 − 9.047561847860013𝑘4
= 1.084400390310856 (𝑖𝑖𝑖)
Dari persamaan (𝑖), (𝑖𝑖), dan (𝑖𝑖𝑖), maka dapat diperoleh penyelesaiannya yaitu
𝑘2 = 42.902459316677110
𝑘3 = −77.392236025978804
𝑘4 = 4.001722708861413
Dengan demikian diperoleh bentuk observer dari model linear sistem endemik penyakit
adalah
−0.0712 −0.0455 −0.7042 −0.3760 0.1
−0.0785 −6.0143 0.6587 0.3305 𝑥̂(𝑡) + 0.1 𝑢(𝑡)
𝑥̂̇(𝑡) = [ ] [ ]
0 6 −4.0143 0 0
0 0 0.4 −0.0393 0
−6.639118333937578
42.902459316677110
+[ ] (𝑦(𝑡) − [1 0 0 0]𝑥̂(𝑡))
−77.392236025978804
4.001722708861413

3.8 Simulasi
Simulasi dilakukan dengan nilai awal parameter yang diberikan pada table berikut:
Parameter 𝛽 𝜎 𝑟1 𝑟2 𝜇 𝜀 𝑝 𝛼 𝜔 𝑞 𝛿

Nilai 0.85 6 4 0.025 1/70 0.8 0.1 0.5 0.1 0.1 1/22

Sedangkan untuk nilai populasi awalnya adalah

 Populasi Nilai awal

Suspectible (𝑥̇ ) 0.493

Exposed (𝑒̇ ) 0.035

Infection Acute (𝑦̇ ) 0.035

Carrier (𝑒̇ ) 0.007

Sehingga dengan menggunakan bantuan software MATLAB (source code terlampir),

didapatkan hasil simulasinya sebagai berikut

Grafik Hasil Simulasi Populasi Suspectible

Grafik Hasil Simulasi Populasi Exposed

 Grafik Hasil Simulasi Populasi Infetion Acute

Grafik Hasil Simulasi Populasi Carrier

Apabila dilakukan simulasi dengan memberikan nilai pada proporsi individu suspectible
yang tervaksinasi (𝑝) yang berbeda, dalam hal ini dilakukan dengan nilai 𝑝 = 0.01, 0.1, 0.5,
dan 0.75, diperoleh hasil simulasinya sebagai berikut

Terlihat dari hasil simulasi individu tingkat Suspectible dengan melihat grafik diatas, jelas
bahwa apabila proporsi individu suspectible yang tervaksinasi kecil atau vaksinasi hanya
dilakukan dalam skala kecil. Maka jumlah individu Suspectible akan meningkat, hal ini
berakibat pula pada tingkat individu Exposed dan Infection Acute yang disajikan dalam
grafik berikut

Jelas bahwa jumlah individu exposed dan infection acute lebih tinggi apabila vaksinasi
dilakukan dalam skala kecil.

Sedangkan apabila nilai 𝑝 diperbesar maka jumlah individu yang sehat atau sembuh
(Recover) akan meningkat mendekati angka 1 dimana menujukkan dalam populasi tersebut
hamper semuanya sehat atau menimbulkan kondisi bebas penyakit.
 BAB IV
PENUTUP
4.1 Kesimpulan
Dari Model SEICR penyebaran HBV mempunyai dua titik kesetimbangan yaitu titik
kesetimbangan bebas penyakit dan titik kesetimbangan endemik virus, titik kesetimbangan
bebas penyakit stabil jika rasio reproduksi dasar 𝑅0 < 1 dan titik kesetimbangan endemik
virus stabil pada saat rasio reproduksi dasar 𝑅0 > 1. Didapatkan pula dengan adanya suatu
variabel pengontrol yang digunakan sebagai matriks input pada sistem, maka didapati
bahwa sistem yang ada merupakan sistem terkontrol. Sistem yang tersedia juga merupakan
sistem yang teramati. Feedback Controller dan juga Observer dari sistem juga dapat
dibentuk sehingga menyebabkan sistemnya teramati dan stabil. Selain itu, dari simulasi
infeksi HBV di Salatiga, diperoleh hasil bahwa untuk mengendalikan infeksi HBV di
Salatiga langkah pertama adalah meningkatkan laju pengobatan. Hal ini dilakukan dengan
memberikan vaksin kepada setiap individu yang rentan terinfeksi dan memberikan
pengobatan kepada penderita Hepatitis B. Keefektivitasan pengobatan sangat berpengaruh
terhadap penyebaran HBV kepada individu yang rentan terinfeksi.
4.2 Saran
Setelah membahas aplikasi pemodelan penyebaran penyakit Hepatitis B dengan
menggunakan model SEICR, penggunaan variabel kontrol berupa terapi untuk
mengahalangi penambahan tingkat penyebaran virus Hepatitis B, yang dalam modelnya
difungsikan sebagai matriks input perlu dikaji ulang. Sebaiknya variabel kontrol tersebut
dapat di ikut sertakan juga kedalam persamaan diferensial atau model SEICR nya.
 DAFTAR PUSTAKA

[1] Ahmad, N., & Kusnanto, H. (2017). Prevalensi infeksi virus Hepatitis B pada bayi dan anak
yang dilahirkan ibu dengan HBsAg positif. Journal of Community Medicine and Public
Health, 33, 515–520.
[2] Anggriani, N., Supriatna, A., Subartini, B., & Wulantini, R. (2015). Kontrol Optimum pada
Model Epidemik SIR dengan Pengaruh Vaksinasi dan Faktor Imigrasi. Jurnal Matematika
Integratif, 11(2), 111–118.
[3] Blass, A. (1978).Bulettin of The American Mathematical Society. American Mathematical
Society, 84, 41–47. Dontwi, I. K., Frempong, N. K., Bentil, D. E., & Adetunde, I. 2010.
Mathematical modeling of Hepatitis C Virus transmission among injecting drug users and
the impact of vaccination. 1–6.
[4] Driessche, P. Van Den, & Watmough, J. (2002).Reproduction numbers and subthreshold
endemic equilibria for compartmental models of disease transmission. 180, 29–48.
[5] Finizio, N., & Ladas, G. (1971).An Introduction to Differential Equations. In The
Mathematical Gazette (Vol. 55). https://doi.org/10.2307/3613341 Giesecke, J. (2017).
Modern Infectious Diseaase Epidemiology. Third Edition. New York: CRC Press.
[6] Kamyad, A. V., Akbari, R., Heydari, A. A., & Heydari, A. (2014). Mathematical Modeling
of Transmission Dynamics and Optimal Control of Vaccination and Treatment for
Hepatitis B Virus. 2014, 80–90.
[7] Khan, T., Ullah, Z., Ali, N., & Zaman, G. (2019). Modeling and control of the hepatitis B
virus spreading using an epidemic model. Chaos, Solitons and Fractals. 124, 1–9.
https://doi.org/10.1016/j.chaos.2019.04.033
[8] Liang, P., Zu, J., & Zhuang, G. (2018). A Literature Review of Mathematical Models of
Hepatitis B Virus Transmission Applied to Immunization Strategies From 1994 to 2015.
28(76), 221–229.
[9] Martin, N. K., Vickerman, P., & Hickman, M. (2011). Mathematical modelling of hepatitis
C treatment for injecting drug users. Journal of Theoretical Biology, 274(1), 58–66.
https://doi.org/10.1016/j.jtbi.2010.12.041
[10] Riwayanti, Gita Purnamasari. 2019. “Analisis Kestabilan Model Matematika SEICR pada
Penyebaran Hepatitis B dengan Pemberian Vaksinasi”. Fakultas Sains dan Teknologi.
Universitas Islam Negeri Sunan Ampel. Surabaya.
[11] Soebiono. 2016. “Sistem Linear dan Optimal Kontrol”. Version 2.2.1. Surabaya : Jurusan
Matematika. Institut Teknologi Sepuluh Nopember.
Lampiran: Source code program simulasi Model Matematika Penyebaran Virus Hepatitis B.
Program dibagi menjadi 2 bagian, yang pertama berupa pembuatan method/function lalu yang
kedua berfungsi untuk memanggil function yang telah dibuat sebelumnya.
Program 1
function [Ddv_Div] = DEdef_scriptfile(~,D)
% Inisialisasi Parameter
beta = 0.85;
sigma= 6;
r1=4;
r2=0.025;
mu=1/70;
epsilon=0.8;
p=0.1;
alpha=0.5;
omega=0.1;
q=0.1;
delta=1/22;

x=D(1);
e=D(2);
y=D(3);
c=D(4);

Ddv_Div = [mu*omega*(1-epsilon*c)-beta*x*(y+alpha*c)-(1-
omega)*p*x+delta*(1-(x+e+y+c))-mu*x;
beta*x*(y+alpha*c)+mu*omega*epsilon*c-sigma*e-mu*e;
sigma*e-r1*y-mu*y;
q*r1*y-r2*c-mu*c];
end

Program 2
% Simulasi Model Matematika SEICR pada Penyebaran Virus Hepatitis B
domain=[0 100]; % Waktu
IC1=0.493; % Nilai awal Suspectible
IC2=0.035; % Nilai awal Exposed
IC3=0.035; % Nilai awal Infection Acute
IC4=0.007; % Nilai awal Carrier

IC=[IC1 IC2 IC3 IC4];

[IVsol, DVsol]=ode23('DEdef_scriptfile',domain,IC);

% Plot grafik Suspectible

figure(1)
plot(IVsol, DVsol(:,1),'k')
title('Suspectible');
ylabel('Populasi');
xlabel('waktu');
grid on;
% Plot grafik Exposed
figure(2)
plot(IVsol, DVsol(:,2),'r')
title('Exposed');
ylabel('Populasi');
xlabel('waktu');
grid on;

% Plot grafik Infection Acute

figure(3)
plot(IVsol, DVsol(:,3),'b')
title('Infection Acute');
ylabel('Populasi');
xlabel('waktu');
grid on;

% Plot grafik Carrier

figure(4)
plot(IVsol, DVsol(:,4),'g')
title('Carrier');
ylabel('Populasi');
xlabel('waktu');
grid on;