Skema Numerik Persamaan Leslie Gower dengan Pemanenan CAUCHY – ISSN: 2086-0382 1 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) = [x1(n), x2(n), …, xk(n)]T dan �⃗�= [g1, g2, …, gk]T . Bentuk umum dari suatu sistem dinamik diskrit adalah �⃗� (n + 1) =�⃗� (�⃗� (n)) (1) Bentuk (1) disebut sebagai sistem dinamik diskrit nonlinear jika gi 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 = (aij), ∀i = 1, 2, …, k adalah matriks yang bernilai real, dan �⃗� (n) = [x1(n), x2(n), … , xk(n)]T , T menyatakan transpose dari vektor A. Jika diberikan nilai awal �⃗� (n) =�⃗�0 maka persamaan (2) memiliki solusi umum �⃗� (n) = An �⃗� (0) (3) Jika matriks A dapat didiagonalkan maka solusi umum persamaan (2) adalah �⃗� (n) = c1𝜆1 𝑛�⃗�1+ … + c3𝜆𝑘 𝑛�⃗�𝑘 (4) Skema Numerik Persamaan Leslie Gower dengan Pemanenan CAUCHY – ISSN: 2086-0382 35 Dengan ci 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 p2 – 4q < 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, p2 – 4q < 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. 𝑑𝑥 𝑑𝑡 =(r1 ─ a1y ─ b1x) x (5a) 𝑑𝑦 𝑑𝑡 = (𝑟2 − 𝑎2 𝑦 𝑥 )y (5b) dengan semua parameter bernilai positif. Kedua parameter, yaitu r1 dan r2 berturut-turut menunjukkan pertumbuhan intrinsik prey dan predator. Parameter a1 merupakan parameter interaksi pemangsaan predator terhadap prey. Parameter b1 adalah parameter antraksi antar prey, sedangkan b2 adalah parameter interaksi antar predator [9].. Model predator prey yang digunakan pada penelitian ini adalah sebagai berikut. 𝑑𝐻 𝑑𝑡 = (r1 ─ a1P ─ b1H) H ─ c1H (6a) 𝑑𝑃 𝑑𝑡 =(𝑟2 − 𝑎2 𝑃 𝐻 )P ─ c2P (6b) Dengan H dan P berturut-turut adalah kepadatan spesies prey dan predator pada waktu t. Pada model ini, diasumsikan 0 < ci < ri, i = 1, 2 untuk mengontrol kepadatan prey dan predator. Model ini mempunyai dua titik kesetimbangan, yaitu E1 dan E2 berturut turut E1 =( 𝑟1−𝑐1 𝑏1 ,0)dan titik kesetimbangan E2 = (H2, P2) dengan H2 = (𝑟1−𝑐1)𝑎2 𝑎1(𝑟2−𝑎2)+𝑎2𝑏1 dan P2 = (𝑟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. Hn+1 = Hn + (k1 – a1Pn – b1Hn) Hn (7a) Pn+1 = Pn +(𝑘2 − 𝑎2 𝑃𝑛 𝐻𝑛 )Pn (7b) Dengan k1 = (r1 – c1) dan k2 = (r2 – c2). 2. Titik Kesetimbangan Model Misalkan titik kesetimbangan dari (7a) dan (7b) adalah (H*, P*) maka haruslah berlaku H* = H* + (k1 – a1 P* – b1 H*) H* (8a) P* = P* +(𝑘2 − 𝑎2 𝑃∗ 𝐻∗ ) P* (8b) Trija Fayeldi 36 Volume 3 No. 4 Mei 2015 Berdasarkan (8a) dan (8b), diperoleh H* = 0 atau k1 – a1 P* – b1 H* = 0 (9a) dan P* = 0 atau (𝑘2 − 𝑎2 𝑃∗ 𝐻∗ )= 0. (9b) Berdasarkan (8b), jelas H* ≠ 0. Substitusikan P* = 0 ke (9a) k1 − a1 (0) – b1 H* = 0 k1 – b1 H* = 0 k1 = b1H* H* = 𝑟1−𝑐1 𝑏1 (10) Sehingga diperoleh E1 =( 𝑟1−𝑐1 𝑏1 ,0). Untuk k1 – a1 P* – b1 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 a1k2H* − a2k1 + a2b1H* = 0 H* (a1k2 + a2b1) = a2k1 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 E2 =( 𝑎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 + (k1 – a1P – b1H) 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 E1 diperoleh: JE1 =( 1 − 𝑘1 − 𝑎1𝑘1 𝑏1 0 1 + 𝑘2 ). Nilai-nilai eigen dari JE1 adalah λ1 = 1 – k1 dan λ2 = 1 + k2. Jelas |λ2| = |1 + k2| = 1 + k2 > 1. Lemma Jika 0 < k1 < 2 maka |λ1| < 1. Bukti: 0 < k1 < 2 −2 < −k1 < 0 −2 + 1 < −k1 + 1 < 0 + 1 −1 < 1 – k1 < 1 |1 – k1| < 1 (terbukti) Berdasarkan lemma tersebut, dapat dituliskan teorema berikut. Teorema 2 Jika k1, k2 > 0 maka 1. E1 saddle jika k1 < 2 2. E1 source jika k1 > 2 3. E1 nonhiperbolik jika k1 = 2 Adapun matriks Jacobi untuk E2 adalah JE2 =( 1 − 𝑏1𝐻 ∗ −𝑎1𝐻 ∗ 𝑘2 2 1 − 𝑘2 ). Misalkan, p dan q berturut-turut adalah trace dan determinan JE2 maka dapat ditulis p = 2 – (b1H* + k2) q = 1 – (b1H* + k2) + k1k2 Lemma 2 Jika k1, k2 > 0 maka 1 – p + q > 0 Bukti: Misalkan k1, k2 > 0. 1 – p + q = 1 – [2 – (b1H* + k2)] + 1 – (k2 + b1H*) + k1k2 = b1H* + k2 + – (b1H* + k2) + k1k2 = k1k2 > 0 (terbukti) Lemma 3 Jika H* > 𝑘2 𝑏1 (k1 – 1) maka q < 1. Bukti: Misalkan, H* > 𝑘2 𝑏1 (k1 – 1) Maka Skema Numerik Persamaan Leslie Gower dengan Pemanenan CAUCHY – ISSN: 2086-0382 37 b1H* > k2(k1 – 1) b1H* > k1k2 – k2 b1H* + k2 > k1k2 b1H* + k2 − k1k2 > 0 −( b1H* + k2) + k1k2 < 0 1 −( b1H* + k2) + k1k2 < 1 q < 1 (terbukti) Lemma 4 Jika H* < 1 2𝑏1 [4 + k2 (k1 – 2)] maka 1 + p + q >0. Bukti: Misalkan H* < 1 2𝑏1 [4 + k2 (k1 – 2)] H* < 1 𝑏1 [2 + 𝑘2 2 (k1 – 2)] b1H* + k2 < 2 + 𝑘1𝑘2 2 2 + 𝑘1𝑘2 2 > b1H* + k2 4 + k1k2 > 2(k2 + b1H*) 4 + k1k2 − 2(b1H* + k2) > 0 1 + [2 – (b1H* + k2)] + [1 – (b1H* + k2) + k1k2] > 0 1 + p + q > 0 (terbukti). Berdasarkan Lemma 2, Lemma 3, dan Lemma 4 dapat dibuat Teorema berikut. Teorema 3 Jika 𝑘2 𝑏1 (k1 – 1) < H* < 1 2𝑏1 [4 + k2(k1 – 2)] maka titik kesetimbangan E2 akan stabil. 4. Simulasi Numerik Berikut akan disajikan simulasi numerik dari sistem (7a) dan (7b) dengan c1 = 0,1; r1 = 0.3; c2= 0.1; r2 = 0.2; a1 = 0.1; b1 = 0.1; a2 = 0.1. Pada simulasi tersebut, diperoleh nilai eigen pada titik kesetimbangan E1 = (3, 0) adalah |λ1| = 0,7 dan |λ2| = 1,1 yang berarti titik kesetimbangan E1 saddle. Adapun nilai eigen yang diperoleh pada titik kesetimbangan E2 = (1.5, 1.5) adalah |λ1| = 0,8832 dan |λ2| = 0,8832 yang berarti titik kesetimbangan E2 stabil. Gambar 1 Plot sistem (7a) dan (7b) dengan c1 = 0,1; r1 = 0.3; c2= 0.1; r2 = 0.2; a1 = 0.1; b1 = 0.1; a2 = 0.1. Simulasi numerik berikutnya menggunakan parameter c1 = 1; r1 = 3.5; c2= 0.1; r2 = 0.2; a1 = 0.3; b1 = 0.2; dan a2 = 0.5. Pada simulasi tersebut, diperoleh nilai eigen pada titik kesetimbangan E1 = (12.5, 0) adalah |λ1| = 1,5 dan |λ2| = 1,1 yang berarti titik kesetimbangan E1 source. Adapun nilai eigen yang diperoleh pada titik kesetimbangan E2 = (9.62, 1.92) adalah |λ1| = 0,89 dan |λ2| = 0,87 yang berarti titik kesetimbangan E2 stabil. Gambar 2 Plot sistem (7a) dan (7b) dengan c1 = 1; r1 = 3.5; c2= 0.1; r2 = 0.2; a1 = 0.3; b1 = 0.2; dan a2 = 0.5. Trija Fayeldi 38 Volume 3 No. 4 Mei 2015 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. PENDAHULUAN TINJAUAN PUSTAKA 1. Sistem Dinamik Diskrit 2. Sistem Dinamik Diskrit Linear 3. Model PredatorPrey Leslie-Gower METODOLOGI PENELITIAN HASIL DAN PEMBAHASAN 1. Diskritisasi Model 2. Titik Kesetimbangan Model 3. Analisis Kestabilan Model 4. Simulasi Numerik KESIMPULAN Bibliography