## SKEMA NUMERIK PERSAMAAN LESLIE GOWER DENGAN PEMANENAN

## Trija Fayeldi
Jurusan Pendidikan Matematika Universitas Kanjuruhan Malang
Email: trija_fayeldi@yahoo.com

## ABSTRAK
Model  predator-prey  merupakan  salah  satu  model  yang  sangat  banyak  diteliti  dan  dimodifikasi Pada  penelitian  ini,  akan  dibangun  suatu  skema  numerik  untuk  menyelesaikan  suatu  model  predator- prey. Model yang dipilih pada penelitian ini adalah model kontinu Leslie-Gower yang telah dimodifikasi pada  pemanenan.  Model  kontinu  tersebut  akan  didiskritisasi  dengan  menggunakan  metode  Euler. Selanjutnya,  akan  diamati  perilaku  model  hasil  diskritisasi  tersebut  dan  pada  akhirnya  akan  diamati apakah  perubahan  parameter  dapat  mempengaruhi  kestabilan  model  dengan  menggunakan  perangkat lunak  Matlab.  Hasil  penelitian  menujukkan  bahwa  perilaku  model  hasil  diskritisasi  konsisten  dengan model kontinunya. Pada simulasi numerik dengan Matlab, ditunjukkan bahwa perubahan parameter akan mempengaruhi kestabilan model.
Kata kunci: predator-prey, Leslie-Gower, stabil, saddle, source, pemanenan.

## PENDAHULUAN
Leslie-Gower memperkenalkan model dengan  populasi  predator  tumbuh  mengikuti model  logistik  [1].  Pada  model  ini,  diasumsikan bahwa  bahwa  carrying  capacity  (kapasitas pendukung) predator sebanding dengan banyaknya  populasi  prey.  Model  ini  kemudian diteliti  dan  dikembangkan  oleh  beberapa peneliti.  Aguirre  mengembangkan  penelitian model  predator-prey  Leslie-Gower  dengan penambahan  efek  [2].  Hasil  penelitian  ini menunjukkan bahwa model Leslie-Gower dengan penambahan efek Allee dapat memiliki dua limit cycle. Aziz mengkaji model Leslie-Gower dengan mengasumsikan fungsi respon yang menyatakan besarnya  pemangsaan  predator  terhadap  prey mengikuti  fungsi  respon  Holling  tipe  II  [3]. Mereka meneliti kestabilan global dan keterbatasan  sistem  tersebut.  Model  yang dikembangakan  oleh  Aziz,  dkk.  kemudian dimodifikasi  dengan  penambahan  waktu  tunda [4]. Dari hasil penelitian, diperoleh bahwa waktu tunda  mempengaruhi  perilaku  dinamik  model tersebut.    Penelitian  ini  dikembangkan  dengan mengamati bifurkasi dari [5]. Chen memperkenalkan  perlindungan  terhadap  prey dan  menunjukkan  bahwa  perlindungan  tidak mempengaruhi  perilaku  model  secara  signifikan [6]. Selain itu, pengaruh efek impulsif pada model diteliti  oleh  Song  [7].  Zhu  melakukan  penelitian model  Leslie-Gower  dengan  fungsi  respon Holling  tipe  II  tanpa  adanya  kontrol  impulsif [8] dan  mengkaji  eksistensi  dan  kestabilan  global solusi periodik positif model dengan mengkonstruksi suatu fungsi Lyapunov..

## TINJAUAN PUSTAKA
1. Sistem Dinamik Diskrit
Bentuk  umum  dari  suatu  sistem  dinamik diskrit dapat dilihat pada Definisi berikut.

## Definisi 1
Misalkan 𝑥⃗ ( n )  =  [ x 1 ( n ), x 2 ( n ),  …, x k ( n )] T dan 𝐺⃗ = [ g 1 , g 2 , …, g k ] T . Bentuk  umum  dari  suatu  sistem  dinamik diskrit adalah
𝑥⃗ ( n + 1) = 𝐺⃗ ( 𝑥⃗ ( n )) (1)
Bentuk  (1)  disebut  sebagai  sistem  dinamik diskrit  nonlinear  jika g i memuat  perkalian antarvariabel tak bebas, ∀ i = 1, 2, …, k .
Definisi 2 Misalkan g adalah fungsi dari R ke R . Titik p *  dikatakan  sebagai  titik  kesetimbangan dari g jika memenuhi g ( p *) = p *.

## 2. Sistem Dinamik Diskrit Linear
Bentuk  umum  persamaan  beda  orde  satu dengan   variabel bebas adalah
𝑥⃗ ( n + 1) = A 𝑥⃗ ( n ) (2)
Dengan A = ( a ij ), ∀ i = 1, 2, …, k adalah matriks yang bernilai real, dan 𝑥⃗ ( n ) = [ x 1 ( n ), x 2 ( n ), … , x k ( n )] T , T menyatakan transpose dari vektor A . Jika diberikan nilai awal 𝑥⃗ ( n ) = 𝑥⃗ 0 maka persamaan (2) memiliki solusi umum
𝑥⃗ ( n ) = A n 𝑥⃗ (0) (3)
Jika  matriks A dapat  didiagonalkan  maka  solusi umum persamaan (2) adalah
𝑥⃗ ( n ) = c 1 𝜆 1 𝑛 𝑣⃗ 1 + … + c 3 𝜆 𝑘 𝑛 𝑣⃗ 𝑘 (4)
Dengan c i adalah  sembarang  konstanta  dan 𝑣 𝑖 ⃗⃗⃗⃗ adalah  vektor  eigen  yang  bersesuaian  dengan nilai eigen λ i .
Teorema 1 Misalkan F (λ) = λ 2 ─ p λ + q maka kondisi-kondisi berikut:
1. F (1) = 1 ─ p + q > 0 2. F ( ─ 1) = 1 + p + q > 0 3.  F (0) = q < 1 merupakan  syarat  perlu  dan  cukup  agar persamaan tersebut memiliki akar-akar karakteristik |λ i |, i =1, 2.
Lemma 1 Jika F (λ) = λ 2 ─ p λ + q dengan 𝜆 1 dan 𝜆 2 adalah dua akar dari F (λ) = 0 maka 1. |λ 1 |  <  1  dan  |λ 2 |  <  1  jika  dan  hanya  jika F (1) > 0 , F (−1) > 0, dan F (0) < 1.
2. |λ 1 | < 1 dan |λ 2 | > 1 atau |λ 1 | > 1 dan |λ 2 | < 1 jika dan hanya jika F (1) > 0 , F (−1) < 0. 3. |λ 1 |  >  1  dan  |λ 2 |  >  1  jika  dan  hanya  jika F (1) > 0 , F (−1) > 0, dan F (0) > 1. 4.  λ 1 dan λ 2 adalah kompleks dan   |λ 1 |  = |λ 2 | = 1 jika dan hanya jika F (1) > 0 , F (1) > 0, F (0) = 1, dan p 2 – 4 q < 0.
Akibat 1 1. Titik kesetimbangan bersifat stabil asimtotik ( sink ) F (1) > 0 , F (−1) > 0, dan F (0) < 1,
2. Titik  kesetimbangan  bersifat  tak  stabil pelana ( saddle ) jika F (1) > 0 , F (−1) < 0,
3. Titik  kesetimbangan  bersifat  tak  stabil ( source ) jika F (1) > 0 , F (−1) > 0, dan F (0) > 1,
4.  Titik  kesetimbangan  merupakan  titik kesetimbangan non-hyperbolic jika F (1) > 0, p 2 – 4 q < 0, dan q = 1.

## 3. Model PredatorPrey Leslie-Gower
Pada  tahun  1948,  Leslie  dan  Gower memperkenalkan  model  dengan  populasi prey tumbuh  mengikuti  model  logistik  sehingga pertumbuhan predator juga  terbatas  disebabkan oleh predator memangsa prey . Selain itu, carrying capacity  predator sebanding  dengan  banyaknya populasi prey ,  sehingga  model  Leslie-Gower dinyatakan sebagai sistem persamaan diferensial nonlinear berikut.
𝑑𝑥 𝑑𝑡 = ( r 1 ─ a 1 y ─ b 1 x ) x (5a)
𝑑𝑦 𝑑𝑡 = (𝑟 2 − 𝑎 2 𝑦 𝑥 ) y (5b)
dengan  semua  parameter  bernilai  positif.  Kedua parameter,  yaitu r 1 dan r 2 berturut-turut
menunjukkan  pertumbuhan  intrinsik prey dan predator .  Parameter a 1 merupakan  parameter interaksi  pemangsaan predator terhadap prey . Parameter b 1 adalah  parameter  antraksi  antar prey ,  sedangkan b 2 adalah  parameter  interaksi antar predator [9] . . Model  predator  prey  yang digunakan  pada  penelitian  ini  adalah  sebagai berikut.
𝑑𝐻 𝑑𝑡 = ( r 1 ─ a 1 P ─ b 1 H ) H ─ c 1 H (6a) 𝑑𝑃 𝑑𝑡 = (𝑟 2 − 𝑎 2 𝑃 𝐻 ) P ─ c 2 P (6b)
Dengan H dan P berturut-turut adalah kepadatan spesies  prey  dan  predator  pada  waktu t .  Pada model ini, diasumsikan  0  < c i < r i , i =  1, 2  untuk mengontrol kepadatan prey dan predator. Model ini mempunyai dua titik kesetimbangan, yaitu E 1 dan E 2 berturut  turut E 1 = ( 𝑟 1 −𝑐 1 𝑏 1 , 0) dan  titik kesetimbangan E 2 =  ( H 2 , P 2 )    dengan H 2 = (𝑟 1 −𝑐 1 )𝑎 2 𝑎 1 (𝑟 2 −𝑎 2 )+𝑎 2 𝑏 1 dan P 2 = (𝑟 1 −𝑐 1 )(𝑟 2 −𝑐 2 ) 𝑎 1 (𝑟 2 −𝑎 2 )+𝑎 2 𝑏 1 . Menurut Zhu  [8],  titik  kesetimbangan  ( H *, P *)  bersifat stabil global.

## METODOLOGI PENELITIAN
Metodologi  penelitian  yang  dilakukan pada  penelitian  ini  adalah  studi  literatur.  Pada penelitian  ini,  persamaan  model  Leslie-Gower kontinu  (6a)  dan  (6b)  didiskritisasi  dengan menggunakan metode Euler. Kemudian, dilakukan  pencarian  titik  kesetimbangan  dan analisis  kestabilan  pada  titik  kesetimbangan tersebut  secara  lokal.  Akhirnya,  dilakukan simulasi  numerik  untuk  beberapa  kasus parameter model.

## HASIL DAN PEMBAHASAN

## 1. Diskritisasi Model
Pandang  kembali  (6a)  dan  (6b).  Dengan menggunakan  beda  maju,  persamaan  (6a)  dan (6b) dapat didiskritisasi menjadi seperti berikut.
H n +1 = H n + ( k 1 – a 1 P n – b 1 H n ) H n (7a) P n +1 = P n + (𝑘 2 − 𝑎 2 𝑃 𝑛 𝐻 𝑛 ) P n (7b)
Dengan k 1 = ( r 1 – c 1 ) dan k 2 = ( r 2 – c 2 ).

## 2.  Titik Kesetimbangan Model
Misalkan  titik  kesetimbangan  dari  (7a) dan  (7b) adalah ( H *, P *) maka haruslah berlaku H * = H * + ( k 1 – a 1 P * – b 1 H *) H * (8a)
P * = P * + (𝑘 2 − 𝑎 2 𝑃 ∗ 𝐻 ∗ ) P * (8b)
Berdasarkan (8a) dan (8b), diperoleh H * = 0 atau k 1 – a 1 P * – b 1 H * = 0 (9a)
dan P * = 0 atau
(𝑘 2 − 𝑎 2 𝑃 ∗ 𝐻 ∗ ) = 0. (9b)
Berdasarkan (8b), jelas H * ≠ 0. Substitusikan P * = 0 ke (9a) k 1 − a 1 (0) – b 1 H * = 0 k 1 – b 1 H * = 0 k 1 = b 1 H* H * = 𝑟 1 −𝑐 1 𝑏 1 (10)
Sehingga diperoleh E 1 = ( 𝑟 1 −𝑐 1 𝑏 1 , 0) .
Untuk k 1 – a 1 P * – b 1 H * = 0 diperoleh
P * = 𝑘 1 −𝑏 1 𝐻 ∗ 𝑎 1 . (11)
Substitusikan (11) ke (9b), diperoleh
(𝑘 2 − 𝑎 2 𝑃 ∗ 𝐻 ∗ ) = 0 (𝑘 2 − 𝑎 2 (𝑘 1 −𝑏 1 𝐻 ∗ ) 𝑎 1 𝐻 ∗ ) = 0 𝑎 1 𝑘 2 𝐻 ∗ −𝑎 2 (𝑘 1 −𝑏 1 𝐻 ∗ ) 𝑎 1 𝐻 ∗ = 0 a 1 k 2 H * − a 2 k 1 + a 2 b 1 H * = 0 H * ( a 1 k 2 + a 2 b 1 ) = a 2 k 1 H * = 𝑎 2 𝑘 1 𝑎 1 𝑘 2 +𝑎 2 𝑏 1 (12)
Substitusikan (12) ke (11), diperoleh
P * = 𝑘 1 −𝑏 1 𝐻 ∗ 𝑎 1 P * = 𝑘 1 −𝑏 1 ( 𝑎2𝑘1 𝑎1𝑘2+𝑎2𝑏1 )
𝑎 1 . P * = 𝑘 1 (𝑎 1 𝑘 2 +𝑎 2 𝑏 1 )−𝑎 2 𝑏 1 𝑘 1 𝑎 1 (𝑎 1 𝑘 2 +𝑎 2 𝑏 1 ) P * = 𝑘 1 (𝑎 1 𝑘 2 +𝑎 2 𝑏 1 −𝑎 2 𝑏 1 ) 𝑎 1 (𝑎 1 𝑘 2 +𝑎 2 𝑏 1 )
P * = 𝑎 1 𝑘 1 𝑘 2 𝑎 1 (𝑎 1 𝑘 2 +𝑎 2 𝑏 1 ) P * = 𝑘 1 𝑘 2 𝑎 1 𝑘 2 +𝑎 2 𝑏 1 .
Sehingga diperleh E 2 = ( 𝑎 2 𝑘 1 𝑎 1 𝑘 2 +𝑎 2 𝑏 1 , 𝑘 1𝑘2 𝑎 1 𝑘 2 +𝑎 2 𝑏 1 ) .

## 3.  Analisis Kestabilan Model
Analisis  kestabilan  model  dilakukan  dengan mengamati perilaku model secara lokal di sekitar titik kesetimbangannya. Misalkan,
F ( H , P ) = H + ( k 1 – a 1 P – b 1 H ) H (13a) G ( H , P ) = P + (𝑘 2 − 𝑎 2 𝑃 𝐻 ) P (13b)
Matriks  Jacobi  dari  sistem  (13a)  dan  (13b) adalah
J = ( 𝜕𝐹 𝜕𝐻 𝜕𝐹 𝜕𝑃 𝜕𝐺 𝜕𝐻 𝜕𝐺 𝜕𝑃 )
= ( 1 + (𝑘 1 − 𝑎 1 𝑃 − 2𝑏 1 𝐻) −𝑎 1 𝑃 𝑎 2 𝑃 2
𝐻 2 1 + (𝑘 2 − 2𝑎 2 𝑃 𝐻 ) ) (14)
.
Untuk E 1 diperoleh:
J E 1 = ( 1 − 𝑘 1 − 𝑎 1 𝑘 1 𝑏 1 0 1 + 𝑘 2 ) .
Nilai-nilai eigen dari  JE1 adalah λ 1 = 1 – k 1 dan λ 2 = 1 + k 2 . Jelas |λ 2 | = |1 + k 2 | = 1 + k 2 > 1.
Lemma Jika 0 < k 1 < 2 maka |λ 1 | < 1.
Bukti: 0 < k 1 < 2 −2 < − k 1 < 0 −2 + 1 < − k 1 + 1 < 0 + 1 −1 < 1 – k 1 < 1 |1 – k 1 | < 1 (terbukti)
Berdasarkan  lemma  tersebut,  dapat  dituliskan teorema berikut.

## Teorema 2
Jika k 1 , k 2 > 0 maka 1. E 1 saddle jika k 1 < 2 2. E 1 source jika k 1 > 2 3. E 1 nonhiperbolik jika k 1 = 2
Adapun matriks Jacobi untuk E 2 adalah
J E 2 = ( 1 − 𝑏 1 𝐻 ∗ −𝑎 1 𝐻 ∗ 𝑘 2 2 1 − 𝑘 2 ) .
Misalkan, p dan q berturut-turut adalah trace dan determinan J E 2 maka dapat ditulis p = 2 – ( b 1 H * + k 2 ) q = 1 – ( b 1 H * + k 2 ) + k 1 k 2 Lemma 2 Jika k 1 , k 2 > 0 maka 1 – p + q > 0
Bukti: Misalkan k 1 , k 2 > 0. 1 – p + q = 1 – [2 – ( b 1 H * + k 2 )] + 1 – ( k 2 + b 1 H *) + k 1 k 2 = b 1 H * + k 2 + – ( b 1 H * + k 2 ) + k 1 k 2 = k 1 k 2 > 0 (terbukti)
Lemma 3 Jika H * > 𝑘 2 𝑏 1 ( k 1 – 1) maka q < 1.
Bukti: Misalkan, H * > 𝑘 2 𝑏 1 ( k 1 – 1) Maka
b 1 H * > k 2 ( k 1 – 1) b 1 H * > k 1 k 2 – k 2 b 1 H * + k 2 > k 1 k 2 b 1 H * + k 2 − k 1 k 2 > 0 −( b 1 H * + k 2 ) + k 1 k 2 < 0 1 −( b 1 H * + k 2 ) + k 1 k 2 < 1 q < 1 (terbukti)
Lemma 4 Jika H * < 1 2𝑏 1 [4 + k 2 ( k 1 – 2)] maka 1 + p + q >0.
Bukti: Misalkan H * < 1 2𝑏 1 [4 + k 2 ( k 1 – 2)]
H * < 1 𝑏 1 [2 + 𝑘 2 2 ( k 1 – 2)] b 1 H * + k 2 < 2 + 𝑘 1 𝑘 2 2 2 + 𝑘 1 𝑘 2 2 > b 1 H * + k 2 4 + k 1 k 2 > 2( k 2 + b 1 H *) 4 + k 1 k 2 − 2( b 1 H * + k 2 ) > 0 1 + [2 – ( b 1 H * + k 2 )] + [1 – ( b 1 H * + k 2 ) + k 1 k 2 ] > 0 1 + p + q > 0 (terbukti).
Berdasarkan  Lemma  2,  Lemma  3,  dan  Lemma  4 dapat dibuat Teorema berikut.
Teorema 3 Jika 𝑘 2 𝑏 1 ( k 1 – 1) < H * < 1 2𝑏 1 [ 4 + k 2 ( k 1 – 2)] maka titik kesetimbangan E 2 akan stabil.

## 4.  Simulasi Numerik
Berikut  akan  disajikan  simulasi  numerik dari sistem (7a) dan (7b) dengan c 1 = 0,1; r 1 = 0.3; c 2 = 0.1; r 2 =  0.2; a 1 =  0.1; b 1 =  0.1; a 2 =  0.1.  Pada simulasi tersebut, diperoleh nilai eigen pada titik kesetimbangan E 1 =  (3,  0)  adalah  |λ 1 |  =  0,7  dan |λ 2 |  =  1,1  yang  berarti  titik  kesetimbangan E 1 saddle.  Adapun  nilai  eigen  yang  diperoleh  pada titik  kesetimbangan E 2 =  (1.5,  1.5)  adalah  |λ 1 |  = 0,8832  dan  |λ 2 |  =  0,8832  yang  berarti  titik kesetimbangan E 2 stabil.
Gambar  1  Plot  sistem  (7a)  dan  (7b)  dengan c 1 = 0,1; r 1 = 0.3; c 2 = 0.1; r 2 = 0.2; a 1 = 0.1; b 1 = 0.1; a 2 = 0.1.
Simulasi  numerik  berikutnya  menggunakan parameter c 1 = 1; r 1 =  3.5; c 2 = 0.1; r 2 =  0.2; a 1 = 0.3; b 1 = 0.2; dan a 2 = 0.5. Pada simulasi tersebut, diperoleh nilai eigen pada titik kesetimbangan E 1 =  (12.5,  0)  adalah  |λ 1 |  =  1,5  dan  |λ 2 |  =  1,1  yang berarti  titik  kesetimbangan E 1 source.  Adapun nilai eigen yang  diperoleh pada titik kesetimbangan E 2 =  (9.62,  1.92)  adalah  |λ 1 |  = 0,89  dan  |λ 2 |  =  0,87  yang  berarti  titik kesetimbangan E 2 stabil.
Gambar  2  Plot  sistem  (7a)  dan  (7b)  dengan c 1 = 1; r 1 = 3.5; c 2 = 0.1; r 2 = 0.2; a 1 = 0.3; b 1 = 0.2; dan a 2 = 0.5.

## KESIMPULAN
Berdasarkan  uraian  pada  hasil  dan pembahasan,  maka  dapat  diambil  kesimpulan bahwa perilaku model hasil diskritisasi konsisten dengan  model  kontinunya.  Selain  itu,  pada simulasi  numerik  dengan  Matlab,  ditunjukkan bahwa perubahan parameter
akan mempengaruhi kestabilan model.

## Bibliography
[1] P. H. Leslie, "Some Furhter Notes on The Use of Matrices in Population
Mathematics," Biometrika, vol. 35, pp. 213- 245, 1948.
[2]  E. P. Aguirre and E. Saez, "Two Limits Cycles in a Leslie-Gower Predator-Prey Model with Additive Allee Effect," Nonlinear Analysis:
Real World Applications, vol. 10, no. 3, pp.
1401-1416, 2009.
[3]  M. A. Aziz-Alaoui and M. D. Okiye,
"Boundedness and Global Stability for a
Predator-Prey Model with Modified Leslie-
Gower and Holling-Type II," Applied Math,
vol. 16, pp. 1069-1075, 2003.
[4]  A. F. Nindjin, M. Aziz-Alaoui and M. Cadivel,
"Analysis of a Predator-Prey Model with
Modified Leslie-Gower and Holling-type II
Schemes with Time Delay," Nonlinear
Analysis: Real World Applications, vol. 7, pp. 1104-1118, 2006.
[5]  R. Yafia and F. E. Adnani, "Limit Cycle and Numerical Simuations for Small and Large
Delays in a Predator-Prey with Modified
Leslie-Gower and Holling-type II Schemes,"
Nonlinear Analysis: Real World
Applications, vol. 9, no. 5, pp. 2055-2067, 2008.
[6]  L. Chen and L. Chen, "Permanence of a Discrete Periodic Volterra Model with
Mutual Interference," Discrete Dynamics in Nature and Society, 2009.
[7]  X. Song and Y. Li, "Dynamic Behaviors of The Periodic Predator-Prey Model with Modified Leslie-Gower Holling-type II
Schemes and Impulsive Effect," Nonlinear Analysis: Real World Applications, vol. 9, no. 1, pp. 64-79, 2008.
[8]  Y. Zhu and K. Wang, "Existence and Global Attractivity of Positive Periodic Solutions for a Predator-Prey Model with Modified Leslie-Gower Holling-type II Schemes,"
Journal of Mathematical Analysis and
Applications, vol. 384, no. 2, pp. 400-408,
2011.
[9]  H. F. Huo, X. Wang and C. C. Chavez,
"Dynamics of a Stage-Structured Leslie- Gower Predator-Prey odel," The Mathematical and Theoretical Biology
Institute, 2011.