Limits: Journal of Mathematics and Its Applications E-ISSN: 2579-8936 P-ISSN: 1829-605X Vol. 17, No. 1, Juli 2020, 19-33 DOI: http://dx.doi.org/10.12962/limits.v17i1.6809
Bifurkasi Periode Ganda dan Neimark-Sacker pada Model Diskret Leslie-Gower dengan Fungsi Respon Ratio-Dependent
Reza Mokodompit 1 , Nurwan 2 , Emli Rahmi 3,* 1,2,3 Jurusan Matematika, Universitas Negeri Gorontalo, Gorontalo, Indonesia e-mail: emlirahmi@ung.ac.id
Diajukan:16 April 2020,  Diperbaiki: 3 Juni 2020, Diterima:17 Juni 2020
Abstrak
Dinamika  model  Leslie-Gower  dengan  fungsi  respon ratio-dependent yang  didiskretisasi menggunakan  skema  Euler  maju  adalah  fokus  utama  pada  artikel  ini.  Analisis  diawali  dengan mengidentifikasi eksistensi dari titik ekuilibrium dan kestabilan lokalnya. Diperoleh empat titik ekuilibrium yaitu  titik  kepunahan  kedua  populasi  dan  titik  kepunahan predator yang  selalu  tidak  stabil,  dan  titik kepunahan prey dan eksistensi kedua populasi yang stabil kondisional. Selanjutnya dipelajari eksistensi dari bifurkasi periode ganda dan Neimark-Sacker di sekitar titik eksistensi kedua populasi sebagai akibat perubahan parameter h (time-step) . Dari hasil analisis ditemukan bahwa bifurkasi periode ganda  terjadi setelah melewati ℎ ℎ atau ℎ ℎ dan bifurkasi Neimark-Sacker terjadi setelah melewati ℎ ℎ . Di akhir pembahasan, diberikan simulasi numerik yang mendukung hasil analisis sebelumnya.
Kata Kunci : model predator-prey , Leslie-Gower, proses diskretisasi, bifurkasi.

## Abstract
In this paper, the dynamics of discrete Leslie-Gower model with ratio-dependent functional response constructed by forward Euler scheme is studied. Firstly, we identify the existence of equilibrium points and their local stability. We obtain four equilibrium points i.e the extinction point of both population and the predator extinction point which are always unstable, the prey extinction point and co-existence point which are conditionally locally stable. Furthermore, the existence of period-doubling and Neimark-Sacker bifurcation are investigated driven by time-step size parameter (h). It is shown that the period-doubling bifurcation occurs when h passes through ℎ or ℎ and Neimark-Sacker bifurcation occurs when h passes through ℎ . We end our result by present some numerical simulations which support previous analytical results.
Keywords: predator-prey model, Leslie-Gower, discretization process, bifurcation.

## 1 Pendahuluan
Dinamika model predator dan prey telah menjadi salah satu topik penting dalam matematika ekologi  karena  keberadaannya  dalam  interaksi  antar  makhluk  di  alam  semesta  [1].  Pada  tahun 1926,  untuk  pertama  kalinya  Lotka  dan  Volterra  secara  terpisah memperkenalkan  suatu  model predator-prey sederhana yang dikenal dengan model Lotka-Volterra [2]. Model Lotka-Volterra mengasumsikan  bahwa predator dan prey tumbuh  secara  eksponensial.  Untuk  memperbaiki asumsi ini, tahun 1948 Leslie dan Gower melakukan modifikasi pada model Lotka-Volterra yang
Bifurkasi Periode Ganda dan Neimark-Sacker pada Model Diskret Leslie-Gower dengan Fungsi Respon Ratio-Dependent
dikenal dengan model Leslie-Gower [3]. Pada model Leslie-Gower kedua populasi diasumsikan tumbuh secara logistik dimana daya dukung predator sebanding dengan kepadatan prey . Namun mengingat  bahwa prey tidak  tersedia  dalam  jumlah  melimpah  di  alam  maka  Aziz-Alaoui  dan Okiye memodifikasi model Leslie-Gower dengan mengasumsikan adanya mangsa lain selain prey ketika kepadatan prey cukup rendah [4].
Dalam  interaksi  antara predator dan prey ,  proses  pemangsaan predator terhadap prey dinyatakan  oleh  fungsi  respon predator [2].  Holling  [5],  memperkenalkan  fungsi  respon  yang hanya bergantung pada ketersediaan prey saja yang disebut prey-dependent . Lebih lanjut, Arditi dan  Ginzburg  [6]  memperkenalkan  fungsi  respon  yang  bergantung  pada  rasio prey terhadap predator yang disebut ratio-dependent . Beberapa penelitian tentang model predator-prey dengan fungsi  respon ratio-dependent seperti  yang  dilakukan  di  [6]–[9]    menunjukkan  bahwa  fungsi respon ratio-dependent sangat sesuai dengan kondisi dimana predator membutuhkan waktu untuk mencari dan mencerna prey -nya.
Kajian  tentang  model predator-pre y  terus  berkembang  dengan  menambahkan  berbagai asumsi  yg  membuat  model  semakin  realistis  dan  kompleks  seperti  adanya  pemanenan  pada predator [10]–[12], adanya efek Allee pada predator [13], [14], infeksi penyakit pada prey [15], dan  sebagainya.  Hal  ini  menyebabkan  solusi  tidak  dapat  ditentukan  secara  analitik,  sehingga dibutuhkan  pendekatan  secara  numerik.    Salah  satu  pendekatan  numerik  yang  dapat  dilakukan untuk mencari solusi dari model kontinu adalah dengan melakukan diskretisasi terhadap model tersebut. Proses diskretisasi ini akan menghasilkan model diskret yang dapat mendekati solusi dari model  kontinu.  Hasil  kajian  di  [16]  dan  [17]  pada  model  diskret  Lotka-Volterra  menunjukkan bahwa dinamika model diskret lebih kaya jika dibandingkan dengan model kontinu.

## 2 Metode Penelitian
Adapun tahapan-tahapan dan prosedur analisis yang dilakukan dalam penelitian ini adalah sebagai berikut.
2.1  Menentukan Model
Aziz-Alaoui dan Okiye [4] mengkaji keterbatasan solusi dan kestabilan global model Leslie- Gower dengan fungsi respon Holling tipe II sebagai berikut:
𝑑𝑥 𝑑𝑡 𝑟 𝑏 𝑥 𝜇𝑦 𝑥 𝛽 𝑥 , 𝑑𝑦 𝑑𝑡 𝜉 𝑘𝑦 𝑥 𝛿 𝑦 , (1)
dengan 𝑥 0 0 dan 𝑦 0 0, 𝑥 dan 𝑦 masing-masing menyatakan kepadatan populasi prey dan predator pada  saat 𝑡; 𝑟,  𝑏 , 𝜇, 𝛽, 𝜉, 𝑘 ,  dan 𝛿 adalah  parameter-parameter  yang  bernilai positif. Parameter tersebut didefinisikan sebagai berikut: 𝑟 dan 𝜉 masing-masing menyatakan laju pertumbuhan intrinsik prey dan predator , 𝑏 adalah  parameter yang  mengukur  kompetisi  yang terjadi antar prey , 𝜇 adalah nilai maksimum yang dapat dicapai oleh laju reduksi per kapita 𝑥 ketika dimangsa 𝑦 , 𝑘 adalah nilai maksimum yang dapat dicapai oleh laju reduksi per kapita 𝑦 sebagai akibat kompetisi antar predator , dan 𝛽, 𝛿 berturut-turut menyatakan  besarnya perlindungan yang diberikan oleh lingkungan terhadap 𝑥 dan 𝑦 .
Fungsi respon predator yang digunakan pada persamaan (1) adalah fungsi respon Holling tipe II yang prey-dependent , yaitu 𝑝 𝑥 . Dengan mengasumsikan bahwa rasio kepadatan populasi prey dan predator berpengaruh  ketika  pemangsaan  terjadi  seperti  yang  diajukan  oleh Arditi dan Ginzburg [6] maka bentuk fungsi respon ratio-dependent yang mengikuti Holling tipe II adalah 𝑝 .
Pilih 𝑏 , 𝑘 , 𝛼 𝜂𝛿 dan dengan menggunakan fungsi respon ratio-dependent , maka persamaan (1)  menjadi:
𝑑𝑥 𝑑𝑡 𝑟𝑥 1 𝑥 𝐾 𝜇𝑥𝑦 𝑥 𝛽𝑦 , 𝑑𝑦 𝑑𝑡 𝜉𝑦 1 𝑦 𝜂𝑥 𝛼 . (2)
Model pada persamaan (2) adalah model Leslie-Gower dengan fungsi respon ratio-dependent . 2.2  Non-Dimensionalisasi
Sebelum didiskretisasi, terlebih dahulu dilakukan non-dimensionalisasi terhadap model pada persamaan  (2).  Non-dimensionalisasi  model  akan  memangkas  beberapa  parameter  sehingga memudahkan  proses  analisis.  Non-dimensionalisasi  dapat  dilakukan  karena  tidak  mengubah kualitatif model [18], [19].
Dengan  memilih 𝑥, 𝑦, 𝑡 → 𝑥, 𝑦, 𝑟𝑡 ,  diperoleh  model  non-dimensional  dari  model pada persamaan (2) sebagai berikut.
𝑑𝑥 𝑑𝑡 𝑥 1 𝑥 𝑚𝑥𝑦 𝑥 𝑦 , 𝑑𝑦 𝑑𝑡 𝑠𝑦 1 𝑏𝑦 𝑥 𝑎 , (3)
22
dengan 𝑚 , 𝑠 , 𝑏 , 𝑎 , dan 𝛼, 𝛽, 𝜂, 𝜇, 𝜉, 𝑟, 𝐾 0 .  Model diskret dari model pada persamaan (3) adalah model yang akan dikaji dalam artikel ini.
2.3  Analisis Model Diskret
Analisis model pada persamaan (3) menggunakan definisi dan dua lemma berikut.
Definisi 1. [20] Suatu titik 𝑥 yang berada dalam domain 𝑓 dikatakan sebagai titik ekulibrium dari suatu persamaan beda 𝑥 𝑓 𝑥 jika 𝑓 𝑥 𝑥 .
Lemma 2. [21] Misalkan 𝜆 dan 𝜆 adalah nilai-nilai eigen dari matriks Jacobi 𝐷𝑓 𝑥 sehingga:
(i) Jika |𝜆 | 1 dan |𝜆 | 1 , maka titik ekuilibrium dari 𝐷𝑓 𝑥 stabil asimtotik lokal (sink).
(ii) Jika |𝜆 | 1 dan |𝜆 | 1 , maka titik ekuilibrium dari 𝐷𝑓 𝑥 tidak stabil (source).
(iii) Jika |𝜆 | 1 dan |𝜆 | 1 (atau, |𝜆 | 1 dan |𝜆 | 1 ), maka titik ekuilibrium dari 𝐷𝑓 𝑥 tidak stabil (saddle).
(iv) Jika |𝜆 | 1 atau |𝜆 | 1 , maka titik kesetimbangan dari 𝐷𝑓 𝑥 non-hiperbolik.
Lemma 3. [22] Jika diasumsikan 𝐹 1 0 dengan 𝜆 , 𝜆 adalah akar-akar karakteristik dari 𝐹 𝜆 𝜆 𝑇𝑟𝑎𝑐𝑒 𝜆 𝐷𝑒𝑡 0 , maka: (i) |𝜆 | 1 dan |𝜆 | 1 jika dan hanya jika 𝐹 1 0 dan Det 1 .
(ii) |𝜆 | 1 dan |𝜆 | 1 , atau |𝜆 | 1 dan |𝜆 | 1 , jika dan hanya jika 𝐹 1 0 . (iii) |𝜆 | 1 dan |𝜆 | 1 jika dan hanya jika 𝐹 1 0 dan Det 1 . (iv) 𝜆 1 dan 𝜆 1 jika dan hanya jika 𝐹 1 0 dan 𝑇𝑟𝑎𝑐𝑒 0,2 .
(v) 𝜆 dan 𝜆 adalah bilangan kompleks dan |𝜆 | |𝜆 | jika dan hanya jika 𝑇𝑟𝑎𝑐𝑒 4𝐷𝑒𝑡
0 dan Det 1.
2.4  Simulasi Numerik
Simulasi numerik dalam artikel ini menggunakan aplikasi Anaconda versi 2019.3.0.0 dan Python  versi  3.7.3150.0.  Anaconda  adalah  distribusi open-source dari  bahasa  pemrograman Python dan R yang digunakan untuk komputasi ilmiah.

## 3 Hasil dan Pembahasan
Model non-dimensional pada persamaan (3) didiskretisasi menggunakan skema Euler maju [20]. Model diskret dari model pada persamaan (3) adalah sebagai berikut:
𝑥 𝑥 ℎ 𝑥 1 𝑥 𝑚𝑥 𝑦 𝑥 𝑦 , 𝑦 𝑦 ℎ 𝑠𝑦 1 𝑏𝑦 𝑥 𝑎 ,
(4)
dengan ℎ 0 adalah  parameter  yang  menyatakan  besar  langkah  dalam  numerik  ( time-step ). Berdasarkan  Definisi  1,  titik  ekuilibrium  dari  model  pada  persamaan  (4)  diperoleh  dengan menyelesaikan secara simultan kedua persamaan berikut:
𝑥 𝑥 ℎ 𝑥 1 𝑥 𝑚𝑥𝑦 𝑥 𝑦 , 𝑦 𝑦 ℎ 𝑠𝑦 1 𝑏𝑦 𝑥 𝑎 ,
(5)
sehingga diperoleh empat titik ekuilibrium sebagai berikut:
𝐸 0,0 ,   𝐸 1,0 , 𝐸 0, 𝑎 𝑏 , 𝐸 𝑥 ∗ , 𝑦 ∗ (6)
dengan
𝑥 ∗ 𝑏 1 𝑎 𝑚 𝐷 ∗ 2 𝑏 1 ,
𝑦 ∗ 𝑥 ∗ 𝑎 𝑏 ,
𝐷 ∗ 𝑚 𝑚 𝑚 𝑚 , 𝑚 , 𝑏 1 𝑎 1 2𝑏 𝑎𝑏 𝑎 1 𝑏 1 .

## 3.1  Eksistensi Titik Ekuilibrium
Berdasarkan Definisi 1, diturunkan Lemma 4 sebagai berikut. Lemma  4. ℝ 𝑥, 𝑦 |𝑥 0, 𝑦 0, 𝑥, 𝑦 ∈ ℝ . Jika 𝐸 𝑥, 𝑦 , 𝑖 0,1,2,3 , adalah solusi dari persamaan (5) , maka 𝐸 adalah titik ekuilibrium dari persamaan (4) jika dan hanya jika 𝐸 ∈ ℝ .
Bukti .  Berdasarkan  kondisi  biologis,  jumlah  populasi  selalu  tak  negatif,  maka  solusi  dari persamaan  (4) harus tak negatif. ■
Dari Lemma 4, diturunkan Teorema 5 berikut, yang menjelaskan eksistensi titik ekuilibrium dari model pada persamaan (4).
Teorema 5. Eksistensi titik ekuilibrium
(i) 𝐸 , 𝐸 , 𝑑𝑎𝑛 𝐸 adalah titik ekuilibrium.
(ii) 𝐸 adalah titik ekuilibrium jika 𝐷 ∗
0 dan (ii.i) 0 𝑎 𝑚 𝑏 1 𝐷 ∗ , atau (ii.ii) 0 𝑎 𝑚 𝐷 ∗ 𝑏 1 .
𝐷 ∗ 0 memberikan dua titik interior, yaitu 𝑥 𝑥 ∗ dan 𝑥 𝑥 ∗ , sedangkan 𝐷 ∗ 0 hanya memberikan satu titik interior yaitu 𝑥 𝑥 𝑥 ∗ .
24
Bukti .
(i) 𝐸 , 𝐸 , 𝐸 ∈ ℝ , maka 𝐸 , 𝐸 , 𝐸 adalah titik ekuilibrium persamaan (4).
(ii) Jika 𝐷 ∗ 0 , maka 𝐷 ∗ 0 , sehingga syarat kepositifan titik ekuilibrium bergantung pada (ii.i) 𝑏 1 𝑎 𝑚 𝐷 ∗ 0 , atau (ii.ii) 𝑏 1 𝑎 𝑚 𝐷 ∗ 0 .
Diperoleh
𝑥 ∗ 𝑏 1 𝑎 𝑚 𝐷 ∗ 2 𝑏 1 0 𝑥 ∗ 𝑏 1 𝑎 𝑚 𝐷 ∗ 2 𝑏 1 0 (7)
Karena 𝑥 ∗ , 𝑥 ∗ 0 , maka 𝑦 ∗ , 𝑦 ∗ 0 , sehingga 𝐸 ∈ ℝ . Akibatnya, 𝐸 adalah titik ekuilibrium persamaan  (4).  Jika 𝐷 ∗ 0 ,  maka 𝐷 ∗ 0 ,  sehingga  syarat  kepositifan  titik  ekuilibrium bergantung pada 𝑏 1 𝑎 𝑚 0 , sehingga diperoleh
𝑥 ∗ 𝑥 ∗ 𝑥 ∗ 𝑏 1 𝑎 𝑚 2 𝑏 1 0. Karena 𝑥 ∗ 0 , maka 𝑦 ∗ 0 . 𝐴 kibatnya 𝐸 adalah titik ekuilibrium persamaan (4). ■ 3.2  Kestabilan Titik Ekuilibrium
Dinamika model pada persamaan (4) di sekitar 𝐸 , 𝑖 0,1,2,3 dipelajari dengan melakukan pelinearan terhadap model pada persamaan (4). Matriks Jacobi dari model pada persamaan (4) di sekitar titik ekuilibrium diberikan oleh:
𝐽 𝑥, 𝑦 ⎣ ⎢ ⎢ ⎢ ⎡1 ℎ 1 2𝑥
ℎ𝑚𝑦 𝑥 𝑦 ℎ𝑚𝑥 𝑥 𝑦 ℎ𝑏𝑠𝑦 𝑥 𝑎 1 ℎ𝑠 2ℎ𝑏𝑠𝑦 𝑥 𝑎 ⎦ ⎥ ⎥ ⎥ ⎤ .
Berdasarkan  Lemma  2  dan  Lemma  3,  diturunkan  beberapa  teorema  berikut  yang menjelaskan kestabilan model pada persamaan (4).
Teorema 6. 𝐸 selalu source.
Bukti . Dari 𝑑𝑒𝑡 𝐽 𝐸 𝜆𝐼 0 , diperoleh 𝜆 1 ℎ ; 𝜆 1 ℎ𝑠. Akibatnya, |𝜆 | 1 dan
|𝜆 | 1 . Sehingga berdasarkan Lemma 2(ii), 𝐸 tidak stabil ( source) .
■
Teorema 7. 𝐸 memiliki tiga kondisi sebagai berikut.
(i) 𝐸 saddle jika 0 ℎ 2 .
(ii) 𝐸 non-hiperbolik jika ℎ 2 .
(iii) 𝐸 source jika ℎ 2 .
Bukti . Dari 𝐷𝑒𝑡 𝐽 𝐸 𝜆𝐼 0 , diperoleh nilai eigen 𝜆 1 ℎ dan 𝜆 1 ℎ𝑠. Perhatikan bahwa  nilai ℎ dan 𝑠 selalu  positif  sehingga |𝜆 | 1. Jika 0 ℎ 2 maka |𝜆 | 1, sehingga menurut  Lemma  2(iii)  titik  ekuilibrium 𝐸 tidak  stabil (saddle). Jika ℎ 2 maka |𝜆 | 1, sehingga menurut Lemma 2(ii) titik ekuilibrium 𝐸 tidak stabil (source). Selanjutnya jika ℎ 2
maka |𝜆 | 1 , sehingga menurut Lemma 2(iv) titik ekuilibrium 𝐸 non-hiperbolik. ■ Teorema 8. 𝐸 memiliki empat kondisi sebagai berikut.
(i) 𝐸 sink jika 1 𝑚 1 dan 0 ℎ𝑠 2 .
(ii) 𝐸 saddle jika
(ii.i) 1 𝑚 1 dan ℎ𝑠 2 ; atau
(ii.ii) 𝑚 1 atau 0 𝑚 1 , dan 0 ℎ𝑠 2 . (iii) 𝐸 source jika 𝑚 1 atau 0 𝑚 1 ; dan ℎ𝑠 2 .
(iv) 𝐸 non-hiperbolik jika 𝑚 1 atau 𝑚 1 ; atau ℎ𝑠 2 . Bukti. Dari 𝐷𝑒𝑡 𝐽 𝐸 𝜆𝐼 0 ,  diperoleh  nilai  eigen 𝜆 1 ℎ ℎ𝑚 dan 𝜆 1 ℎ𝑠.
Perhatikan bahwa ℎ, 𝑠, dan 𝑚 selalu positif, maka kondisi-kondisi berikut akan berlaku:  Jika 1 𝑚 1 , maka dengan operasi aljabar diperoleh 1 1 ℎ 1 𝑚 1 atau |λ | 1 dan
jika 𝑚 1
atau 0 𝑚 1 maka diperoleh |λ | 1 . Jika 0 ℎ𝑠 2, maka diperoleh |λ |
1 dan  jika ℎ𝑠 2 ,  maka |λ | 1 .    Selanjutnya  jika 𝑚 1 atau 𝑚 1 atau ℎ𝑠
2 maka |λ | 1 atau |λ | 1 . Dengan menggunakan Lemma 2(i)-(iv), Teorema 8 terbukti.      ■
Teorema 9. 𝐽𝑖𝑘𝑎 𝑥 ∗ √ maka 𝐸 memiliki tiga kondisi kestabilan yang bergantung pada parameter ℎ sebagai berikut:
(i) 𝐸 sink jika memenuhi salah satu kondisi berikut:
(i.i) Λ 0 dan 0 ℎ ℎ ; atau (i.ii) Λ 0 𝑑𝑎𝑛 0 ℎ ℎ ,
(ii) 𝐸 saddle jika Λ 0 dan ℎ ℎ ℎ ,
26
(iii) 𝐸 source jika memenuhi salah satu kondisi berikut:
(iii.i) Λ 0 dan ℎ ℎ ; 𝑎𝑡𝑎𝑢
(iii.ii) Λ 0 dan ℎ ℎ ,
(iv) 𝐸 non-hiperbolik jika memenuhi salah satu kondisi berikut:
(iv.i) Λ 0 dan ℎ ℎ  𝑎𝑡𝑎𝑢 ℎ ℎ
(iv.ii) Λ 0 dan ℎ ℎ , dengan ℎ √ ,    ℎ ,     ℎ √ , dan Λ 𝑥 ∗ 𝑠 ∗ ∗ ∗ ,   𝛬 𝑠𝑥 ∗ ∗ ∗
,    𝛬 𝛬 4𝛬 .
Bukti. Matriks Jacobi 𝐽 𝑥, 𝑦 dihitung di titik eksistensi kedua populasi, 𝐸 𝑥 ∗ , 𝑦 ∗ ) diberikan oleh:
𝐽 𝐸 ⎣ ⎢ ⎢ ⎢ ⎡1 ℎ 𝑚𝑏𝑥 ∗ 𝑥 ∗ 𝑎 𝑏 1 𝑥 ∗ 𝑎 𝑥 ∗ ℎ𝑚𝑏 𝑥 ∗ 𝑏 1 𝑥 ∗ 𝑎 ℎ𝑠 𝑏 1 ℎ𝑠 ⎦ ⎥ ⎥ ⎥ ⎤ .
Persamaan  karakteristiknya  dapat  dituliskan: 𝜆 𝑇𝑟𝑎𝑐𝑒 𝐽 𝐸 𝜆 𝐷𝑒𝑡 𝐽 𝐸 0, dengan 𝑇𝑟𝑎𝑐𝑒 𝐽 𝐸 2 ℎΛ dan 𝐷𝑒𝑡 𝐽 𝐸 1 ℎΛ ℎ Λ . Misalkan 𝐹 𝜆 𝜆 𝑇𝑟𝑎𝑐𝑒 𝐽 𝐸 𝜆 𝐷𝑒𝑡 𝐽 𝐸 , maka 𝐹 1 Λ ℎ , 𝐹 1 Λ ℎ 2Λ ℎ 4, dan 𝐷𝑒𝑡 𝐽 𝐸 Λ ℎ Λ ℎ 1. Dari  persamaan  karakteristik  diperoleh  nilai  eigen 𝜆 , 1 Λ √Λ. Nilai dari 𝐹 1 dapat dinyatakan sebagai berikut:
(i) 𝐹 1 0 jika Λ >0 dan 0 ℎ ℎ atau ℎ ℎ .
(ii) 𝐹 1 0 dan 𝑇𝑟𝑎𝑐𝑒 𝐽 𝐸 0,2 jika Λ >0 dan ℎ ℎ atau ℎ ℎ .
(iii) 𝐹 1 0 jika Λ >0 dan ℎ ℎ ℎ . Dari 𝐷𝑒𝑡 𝐽 𝐸 1 diperoleh ℎ , dengan ℎ ℎ ℎ .
Pada kasus Λ 0, mengakibatkan 𝜆 , ∈ ℂ. Adapun nilai 𝐷𝑒𝑡 𝐽 𝐸 dinyatakan sebagai berikut:
(i) 𝐷𝑒𝑡 𝐽 𝐸
1 dan 𝜆 , 1 jika 0 ℎ ℎ . (ii) 𝐷𝑒𝑡 𝐽 𝐸 1 dan 𝜆 , 1 jika ℎ ℎ . (iii) 𝐷𝑒𝑡 𝐽 𝐸 1 dan 𝜆 ,
1 jika ℎ ℎ .
Dengan menggunakan Lemma 2 dan 3, maka Teorema 9 terbukti.
■
Titik ekuilibrium 𝐸 dan 𝐸 merupakan titik yang tidak stabil secara struktural, sedangkan titik ekuilibrium 𝐸 dan 𝐸 merupakan titik yang stabil secara kondisional. Hal ini berarti bahwa kedua populasi tidak akan punah secara bersamaan, atau populasi predator tidak akan mungkin punah.
3.3  Eksistensi Bifurkasi Periode Ganda dan Neimark-Sacker
Selanjutnya akan ditunjukkan fenomena perubahan dinamika di sekitar titik ekuilibrium eksistensi kedua  populasi  yaitu  bifurkasi  periode  ganda  dan  bifurkasi  Neimark-Sacker.  Bifurkasi  periode ganda  terjadi  jika  perubahan  suatu  nilai  parameter  mengakibatkan  perubahan  kestabilan  titik ekuilibrium disertai munculnya solusi periodik dengan periode-2. Selanjutnya, bifurkasi Neimark- Sacker  terjadi  jika  perubahan  suatu  nilai  parameter  mengakibatkan  perubahan  kestabilan  titik ekuilibrium disertai munculnya suatu limit-cycle yang mengisolasi titik ekuilibrium tersebut. Pada artikel  ini  ditunjukkan  bahwa  kedua  bifurkasi  ini  terjadi  pada  titik  eksistensi  kedua  populasi dengan menggerakkan nilai dari parameter h atau time-step nya.
Teorema 10. Titik ekuilibrium 𝐸 𝑥 ∗ , 𝑦 ∗ kehilangan kestabilannya:
(i) via bifurkasi periode ganda ketika ℎ melewati ℎ atau ℎ ,
(ii) via bifurkasi Neimark-Sacker ketika ℎ melewati ℎ .
Bukti . Kajian lengkap tentang teori bifurkasi dapat dilihat pada Elaydi [20]. Berdasarkan Teorema 9, 𝐸 non-hiperbolik pada saat Λ 0 dan ℎ ℎ  atau ℎ ℎ ; dan Λ 0 dan ℎ ℎ . Jika Λ 0
dan ℎ ℎ  atau ℎ ℎ maka  salah  satu  nilai  eigen  bernilai  -1.  Akibatnya  solusi  periode-2 muncul ketika h melewati ℎ ℎ  atau ℎ ℎ sehingga terjadi bifurkasi periode ganda. Pada saat
Λ 0 dan ℎ ℎ ,  terdapat  sepasang  nilai  eigen  kompleks  dengan  modulusnya  bernilai  1. Akibatnya limit-cycle muncul ketika h melewati ℎ . Dengan demikian terjadi bifurkasi Neimark- Sacker.
■
3.4  Simulasi Numerik
Pada tahapan ini diberikan beberapa simulasi numerik yang bersesuaian dengan hasil analisis yang  dipaparkan  sebelumnya.  Simulasi  numerik  diawali  dengan  menunjukkan  kestabilan  lokal untuk 𝐸 dan 𝐸 , dan dilanjutkan dengan memberikan simulasi eksistensi dari bifurkasi periode ganda dan Neimark-Sacker. Parameter yang digunakan dipilih berdasarkan kondisi yang diberikan pada hasil analisis sebelumnya dikarenakan tidak adanya data rill yang bersesuaian dengan model. Oleh karena itu, diberikan nilai parameter sebagai berikut.
𝑚 1.1,   𝑠 0.9, 𝑏 0.5, 𝑎 0.1, ℎ 0.1 (8)
Dengan menggunakan nilai parameter (8) diperoleh titik ekuilibrium 𝐸 0,2 yang stabil asimtotik  lokal,  lihat  Gambar  1.  Hal  ini  menunjukkan  bahwa  populasi prey akan  punah  dan eksistensi  dari predator tetap  dipertahankan  jika  nilai  awal  dari  kedua  populasi  cukup  dekat
dengan  titik  ekuilibrium  tersebut.  Kondisi  ini  bersesuaian  dengan  asumsi  sebelumnya  bahwa predator memiliki  alternatif  makanan  lain  selain prey .  Ketika 𝑚 diturunkan  menjadi  0.3,  titik ekuilibrium 𝐸 berubah  menjadi  tidak  stabil  yang  disertai  dengan  munculnya  titik  ekuilibrium tunggal 𝐸 yang  stabil  asimtotik  lokal  seperti  yang  ditunjukkan  oleh  Gambar  2.  Hal  ini menunjukkan bahwa jika laju pemangsaan prey oleh predator diturunkan maka populasi prey dan predator akan terjaga eksistensinya.
(a)
(b)
(c)
Gambar 1. Simulasi numerik model pada persamaan (4) dengan nilai parameter pada persamaan (8) (a) Potret fase disekitar 𝐸 (b) dan (c) Simulasi deret waktu 𝑥 dan 𝑦 terhadap 𝑛 . (a) (b) (c) Gambar 2. Simulasi numerik model pada persamaan (4) dengan nilai parameter pada persamaan (8) dan 𝑚 0.3 (a) Potret fase disekitar 𝐸 (b) dan (c) Simulasi deret waktu 𝑥 dan 𝑦 terhadap 𝑛
Untuk menunjukkan eksistensi dari bifurkasi periode ganda, dipilih nilai parameter seperti pada persamaan (8) dengan mengganti nilai 𝑚 0.3 , dan 𝑠 0.33 . Selanjutnya dilakukan iterasi sebanyak 5000 kali dan di plot 200 titik terakhir. Proses ini dilakukan untuk 2.9 ℎ 3.4 dengan tujuan untuk melihat konvergensi dari solusi dengan variasi nilai ℎ , sehingga diperoleh diagram bifurkasi seperti pada Gambar 3.
Gambar 3. Diagram bifurkasi periode ganda pada model pada persamaan (4) dengan nilai parameter pada persamaan (8), 𝑚 0.3 , dan 𝑠 0.33 (a) 𝒉 𝟑. 𝟏
(b) 𝒉 𝟑. 𝟑
(c) 𝒉 𝟑. 𝟑𝟏 (d) 𝒉 𝟑. 𝟑𝟐 Gambar 4. Solusi periodik pada model pada persamaan (4) dengan nilai parameter pada
persamaan (8), 𝑚 0.3 , dan 𝑠 0.33
Gambar 5. Diagram bifurkasi Neimark-Sacker pada model pada persamaan (4) dengan nilai
parameter pada persamaan (8) dan 𝑚 0.3
(a) 𝒉 𝟐. 𝟑 (b) 𝒉 𝟐. 𝟑𝟏 Gambar 6. Potret fase model pada persamaan (4) dengan nilai parameter pada persamaan
(8) dan 𝑚 0.3
Perhatikan bahwa solusi yang pada awalnya konvergen ke titik ekuilibrium 𝐸 , kehilangan kestabilannya via bifurkasi periode ganda dan menuju ke solusi periode-2. Dari Gambar 3 juga diperlihatkan  munculnya  solusi  periodik  dengan  berbagai  periode  dari  solusi.  Untuk menunjukkannya, dipilih beberapa nilai ℎ . Pada saat ℎ 3.1 solusi konvergen ke solusi periode- 2, lihat Gambar 4(a). Nilai ℎ dinaikkan menjadi 3.3 sehingga muncul solusi periode-4 yang stabil seperti  pada  Gambar  4(b).  Ketika ℎ menjadi  3.31,  solusi  konvergen  ke  solusi  periode-8  yang ditunjukkan oleh Gambar 4(c), dan simulasi untuk solusi periodik diakhiri dengan memilih ℎ 3.32 dimana  muncul  solusi  periode-16  seperti  pada  Gambar  4(d).  Secara  biologis,  hal  ini
menunjukkan  bahwa  meskipun  titik  ekuilibrium 𝐸 tidak  stabil,  namun  eksistensi  dari  kedua populasi terjaga secara periodik.
Simulasi  numerik  ini  diakhiri  dengan  menunjukkan  eksistensi  dari  bifurkasi  Neimark- Sacker. Dengan memilih nilai parameter pada persamaan (8) dan 𝑚 0.3 , serta menggerakkan nilai ℎ diperoleh diagram bifurkasi seperti Gambar 5. Solusi yang pada awalnya konvergen ke titik ekuilibrium 𝐸 kehilangan kestabilannya via bifurkasi Neimark-Sacker yang ditandai munculnya limit-cycle yang mengisolasi titik ekuilibriumnya. Pada saat ℎ 2.3 , solusi masih konvergen ke
𝐸 yang ditunjukkan oleh Gambar 6(a), dan ketika ℎ 2.31 , titik ekuilibrium 𝐸 menjadi tidak stabil  disertai  munculnya limit-cycle yang  mengisolasi  titik  ekuilibrium  tersebut,  lihat  Gambar 6(b).  Meskipun  titik  ekuilibrium 𝐸 kehilangan  kestabilannya,  namun  solusi  bergerak  secara periodik  di  sekitar limit-cycle .  Secara  biologis,  hal  ini  berarti  bahwa  dengan  kondisi  awal kepadatan kedua populasi yang cukup dekat dengan titik ekuilibrium 𝐸 , eksistensi kedua populasi tetap terjaga.

## 4 Simpulan
Model pada persamaan (4) merepresentasikan interaksi antara dua spesies dengan hubungan mangsa-pemangsa dimana prey dan predator masing-masing tumbuh secara logistik dengan daya dukung prey bergantung kepada lingkungannya dan daya dukung predator adalah kepadatan prey dan mangsa alternatifnya. Dari hasil analisis ditunjukkan bahwa kepunahan dari predator tidak akan  pernah  tercapai,  sedangkan  kepunahan  dari prey bisa  tercapai  dengan  kondisi  biologis tertentu.  Fenomena  menarik  juga  ditunjukkan  secara  biologis  di  sekitar  titik  kesetimbangan dimana eksistensi dari kedua populasi terjaga meskipun titik interior kehilangan kestabilannya via bifurkasi  periode-ganda  dan  Neimark-Sacker.  Kepadatan  populasi  bergerak  secara  periodik  di sekitar titik interior.
Meskipun dinamika dari model ini telah dibahas cukup lengkap, namun ada satu hal menarik yang tidak dieksplorasi pada artikel ini. Dari simulasi yang ditunjukkan pada Gambar 3, dapat dilihat  bahwa  terdapat  kemungkinan  munculnya  solusi  periode-3,  dimana  menurut  Elaydi  [20] merupakan salah satu indikator adanya strange-attractor yang memicu terjadinya chaos .

## 5 Ucapan Terima Kasih
Terimakasih kepada Bapak Hasan S. Panigoro atas masukan dan sarannya dalam penelitian dan penulisan artikel ini.
32

## 6 Daftar Pustaka
[1] A. A. Berryman, “The Orgins and Evolution of Predator-Prey Theory,” Ecology , vol. 73, no. 5, pp. 1530–1535, Oct. 1992, doi: 10.2307/1940005.
[2] J. D. Murray, Mathematical Biology: An Introduction , 3rd ed. New York, NY: Springer New York, 2002.
[3] P.  H.  Leslie,  “Some  Further  Notes  on  the  Use  of  Matrices in  Population  Mathematics,” Biometrika , vol. 35, no. 3–4, pp. 213–245, 1948, doi: 10.1093/biomet/35.3-4.213.
[4] M. A. Aziz-Alaoui and M. Daher Okiye, “Boundedness and global stability for a predator- prey model with modified Leslie-Gower and Holling-type II schemes,” Appl. Math. Lett. , vol. 16, no. 7, pp. 1069–1075, 2003, doi: 10.1016/S0893-9659(03)90096-6.
[5] C. S. Holling, “Some Characteristics of Simple Types of Predation and Parasitism,” Can. Entomol. , vol. 91, no. 07, pp. 385–398, 1959, doi: 10.4039/Ent91385-7.
[6] R. Arditi and L. R. Ginzburg, “Coupling in predator-prey dynamics: Ratio-Dependence,” J. Theor. Biol. , vol. 139, no. 3, pp. 311–326, 1989, doi: 10.1016/S0022-5193(89)80211-5.
[7] I. Hanski, “The functional response of predators: Worries about scale,” Trends Ecol. Evol. , vol. 6, no. 5, pp. 141–142, May 1991, doi: 10.1016/0169-5347(91)90052-Y.
[8] M. J. Bishop, B. P. Kelaher, M. P. L. Smith, P. H. York, and D. J. Booth, “Ratio-dependent response  of  a  temperate  Australian  estuarine  system  to  sustained  nitrogen  loading,” Oecologia , vol. 149, no. 4, pp. 701–708, Oct. 2006, doi: 10.1007/s00442-006-0481-5.
[9] A. Suryanto, I. Darti, H. S. Panigoro, and A. Kilicman, “A Fractional-Order Predator–Prey Model with Ratio-Dependent Functional Response and Linear Harvesting,” Mathematics , vol. 7, no. 11, p. 1100, Nov. 2019, doi: 10.3390/math7111100.
[10]  N.  Hasan,  R.  Resmawan,  and  E.  Rahmi,  “Analisis  Kestabilan  Model  Eko-Epidemiologi dengan Pemanenan Konstan pada Predator,” J. Mat. Stat. dan Komputasi , vol. 16, no. 2, p. 121, Dec. 2019, doi: 10.20956/jmsk.v16i2.7317.
[11]  H.  S.  Panigoro,  “Analisis  Dinamik  Sistem  Predator-Prey  Model  Leslie-Gower  dengan Pemanenan Secara Konstan terhadap Predator,” EULER , vol. 2, no. 1, pp. 1–12, 2014.
[12]  P. D. Ernawati and I. Darti, “Stability analysis of the Euler discretization for the harvesting Leslie-Gower predator-prey model,” Int. J. Pure Apllied Math. , vol. 105, no. 2, pp. 213– 221, Nov. 2015, doi: 10.12732/ijpam.v105i2.8.
[13]  H. S. Panigoro, E. Rahmi, N. Achmad, and S. L. Mahmud, “The Influence of Additive Allee Effect and Periodic Harvesting to the Dynamics of Leslie-Gower Predator-Prey Model,” Jambura J. Math. , vol. 2, no. 2, pp. 87–96, Mar. 2020, doi: 10.34312/jjom.v2i2.4566.
[14]  H. S. Panigoro and D. Savitri, “Bifurkasi Hopf pada Model Lotka-Volterra Orde-Fraksional dengan Efek Allee Aditif pada Predator,” Jambura J. Biomath. , vol. 1, no. 1, pp. 16–24, Jun. 2020, doi: 10.34312/jjbm.v1i1.6908.
[15]  H.  S.  Panigoro,  A.  Suryanto,  W.  M.  Kusumawinahyu,  and  I.  Darti,  “Dynamics  of  a Fractional-Order  Predator-Prey  Model  with  Infectious  Diseases  in  Prey,” Commun. Biomath. Sci. , vol. 2, no. 2, pp. 105–117, 2019, doi: 10.5614/cbms.2019.2.2.4.
[16]  M. Renisagayaraj, A. G. M. Selvam, and M. Meganathan, “Dynamics in a Discrete Prey- Predator System,” Int. J. Eng. Res. Dev. , vol. 6, no. 5, pp. 1–5, 2013.
[17]  M. Zhao, Z. Xuan, and C. Li, “Dynamics of a discrete-time predator-prey system,” Adv. Differ. Equations , vol. 2016, no. 1, p. 191, Dec. 2016, doi: 10.1186/s13662-016-0903-6.
[18]  L.  Edelstein-Keshet,  “Mathematical  models  in  biology,”  in Classics in applied mathematics , 2005.
[19]  H. S. Panigoro and E. Rahmi, “Modifikasi sistem predator-prey: dinamika model Leslie- Gower dengan daya dukung yang tumbuh logistik,” in SEMIRATA MIPAnet , 2017, pp. 94– 103.
[20]  S. N. Elaydi, Discrete Chaos: With Applications in Science and Engineering , 2nd ed., vol. 2007. San Antonio, Texas: Chapman & Hall/ CRC, 2007.
[21]  M. A. M. Abdelaziz, A. I. Ismail, F. A. Abdullah, and M. H. Mohd, “Bifurcations and chaos in a discrete SI epidemic model with fractional order,” Adv. Differ. Equations , vol. 2018, no. 1, p. 44, Dec. 2018, doi: 10.1186/s13662-018-1481-6.
[22]  L. G. Yuan and Q. G. Yang, “Bifurcation, invariant curve and hybrid control in a discrete- time predator–prey system,” Appl. Math. Model. , vol. 39, no. 8, pp. 2345–2362, Apr. 2015, doi: 10.1016/j.apm.2014.10.040.