PENYELESAIAN NUMERIK MODEL PREDATOR- PREY DENGAN SKEMA BEDA HINGGA TAK- STANDAR
Rina Reorita 1 , Renny 2
1,2 Jurusan Matematika, Universitas Jenderal Soedirman
Email : 1 rina.reorita@unsoed.ac.id
Abstract. Interaction between predator and prey can be represented as a system of non- linear differential equation which is difficult to be solved analytically. In this research, a predator-prey model with an addition of harvesting factor is discretized into a system of difference  equation  using  non-standard  finite  difference  scheme.  The  analysis  result shows  that  the  developed  scheme  has  qualitative  property  which  is  consistent  to  the continuous system.
Keywords: predator prey model, non-standard finite difference scheme, harvesting
Abstrak. Interaksi  yang  terjadi  antara  predator  dan  prey  dapat  direpresentasikan  ke dalam sistem persamaan diferensial non linier yang sulit untuk diselesaikan secara eksak. Pada  penelitian  ini,  model  predator-prey  yang  ditambahkan  dengan  faktor  pemanenan didiskritisasi ke dalam bentuk sistem persamaan beda dengan menggunakan skema beda hingga  tak-standar.  Hasil  analisis  menunjukkan  bahwa  skema  yang  dikembangkan memiliki sifat kualitatif yang konsisten dengan sistem kontinu.
Kata kunci: model predator-prey, skema beda hingga tak-standar, pemamenan

## I. PENDAHULUAN
Populasi dalam suatu ekosistem tidak dapat berdiri sendiri, tetapi saling berinteraksi. Bentuk  interaksi  yang  terjadi  di  dalam  ekosistem  antara  lain  berupa  predasi  dan  persaingan memangsa ( competitors prey ). Predasi merupakan proses pemangsaan terhadap satu spesies, yaitu satu spesies memangsa spesies lain. Dalam predasi, individu  yang terlibat mempunyai fungsi  sebagai  mangsa  ( prey )  dan  pemangsa  (predator).  Seiring  dengan  interaksi  tersebut terdapat  rangkaian  peristiwa  memakan  dan  dimakan  yang  menjadikan  ekosistem  tetap seimbang.
Dalam  bidang  pemodelan  matematika,  model  predator -prey masih  menjadi  hal  yang sangat menarik untuk dikaji. Dalam [1], model predator -prey yang pertama dicetusakan adalah adalah  model  Lotka  Volterra  pada  sekitar  tahun  1920.  Model  ini  menggambarkan  interaksi antara  satu  spesies predator dan  satu  spesies prey .  Model  ini  berbentuk  sistem  persamaan diferensial  non  linier  sehingga  sulit  untuk  diselesaikan  secara  eksak.  Dalam  [2],  model predator-prey dimodifikasi dengan menambahkan faktor pemanenan. Sifat-sifat dinamik dari model tersebut dianalisis secara rinci.
Salah satu alternatif lain untuk menyelesaikan sistem persamaan diferensial non linier adalah  dengan  metode  pendekatan  numerik.  Kelebihan  dari  metode  numerik  adalah penyelesaian yang dihasilkan dapat dibuat sedekat mungkin dengan nilai sesungguhnya. Selain itu,  tampilan  perhitungan  juga  dapat  disimulasikan.  Dengan  demikian,  pada  penelitian  ini
http://dx.doi.org/10.14710/jfma.v2i1.24
JOURNAL OF FUNDAMENTAL MATHEMATICS AND APPLICATIONS (JFMA) VOL. 2 NO. 1 (Jun 2019) Available online at www.jfma.math.fsm.undip.ac.id
model  predator -prey akan  diselesaikan  secara  numerik  dengan  menggunakan  skema  beda hingga  tak-standar.  Skema  beda  hingga  tak-standar  merupakan  metode  numerik  yang menstransformasikan sistem persamaan diferensial ke dalam bentuk persamaan beda dengan mengaproksimasi turunan pertama dengan generalisasi beda maju dan mendekati suku lainnya dengan pendekatan non-local [3].
Langkah-langkah yang dilakukan dalam penelitian ini antara lain menentukan asumsi- asumsi  model  yang  sesuai  dan  mengidentifikasi  variabel-variabel  dan  parameter  yang digunakan dalam model, menurunkan model, menentukan penyelesaian model secara kualitatif dengan  metode  skema  beda  hingga  tak-standar.  Analisis  dilakukan  untuk  mengetahui kestabilan titik kesetimbangan.

## II. FORMULASI MODEL
Pembuatan  model  matematika  yang  bersumber  pada  kondisi  riil  diharapkan  dapat digunakan  untuk  memprediksi  perilaku  model  untuk  jangka  waktu  yang  lama  ( long  term behaviour ).   Dengan  adanya model matematika  yang merupakan  gambaran fenomena dunia nyata  ( real  world ),  dapat  diprediksi  perilaku  model  untuk  waktu  yang  akan  datang,  apakah unsur-unsur  yang  terdapat  dalam  model  akan  semakin  besar  dan  bertambah  tanpa  batas, berosilasi (bergerak naik turun) atau mendekati suatu nilai tertentu [4].
Misalkan ) ( t N menyatakan  jumlah prey pada  saat t dan ) ( t P menyatakan  jumlah predator pada saat t. Diasumsikan jika tidak ada predator, prey tumbuh secara logistik dengan carrying capacity K , sehingga model predator-prey dapat dituliskan sebagai berikut:
mP dP caNP dt dy aNP K N rN dt dN − − = −       − = 1 (1)
dengan r , K , a , c , dan d konstanta positif. Parameter r menyatakan laju pertumbuhan intrinsik prey, a menyatakan laju predasi, c menyatakan efisiensi predasi yang membagi laju kelahiran predator maksimum per kapita terhadap laju konsumsi maksimum per kapita, d menyatakan laju kematian predator, dan m menyatakan laju pemanenan [2].
Sistem persamaan (1) memiliki 3 titik kesetimbangan, yaitu TE 1 (0, 0), ( ) 0 , 2 K TE , dan
            + − + Kca m d a r ca m d
TE 1 , 3 . Agar komponen pada TE 3 positif, maka haruslah ( d + m ) < Kca .
Untuk  selanjutnya,  sistem  persamaan  (1)  akan  didiskritisasi  dengan  menggunakan skema beda hingga tak standar dengan pendekatan non lokal seperti dalam [5]. Variabel 0  t didiskritisasi pada titik-titik t n = nh dengan n = 0,1, 2, .... Pada penelitian ini, dipilih h h = ) (  dengan h >  0  menunjukkan  ukuran  langkah.  Akan  tetapi,  pada  bagian  Simulasi  akan dibandingkan  dengan r e h rh 1 ) ( − =  ,  sesuai  dengan  yang  dilakukan  dalam  [6].  Proses diskritisasi model untuk laju populasi prey terhadap waktu adalah sebagai berikut.
JOURNAL OF FUNDAMENTAL MATHEMATICS AND APPLICATIONS (JFMA) VOL. 2 NO. 1 (Jun 2019) Available online at www.jfma.math.fsm.undip.ac.id
. 1 1 ) ( 1 1 1 1 1 1 n n n n n n n n n n n n P aN K N rN h N N P aN K N rN h N N + + + + + + −       − = − −       − = −  Sementara itu, proses diskritisasi model untuk laju populasi predator adalah sebagai berikut . ) ( 1 1 1 1 1 1 + + + + + + − − = − − − = − n n n n n n n n n n n n mP dP P caN h P P mP dP P caN h P P  Dengan demikian, diperoleh bentuk sistem diskrit dari sistem persamaan (1) adalah ( ) ( ) ( ) . 1 1 1 1 1 1 1 mh P dh hcaN P haP N K hr N hr N n n n n n n n + + − =       + + + = + + + (2)

## III. ANALISIS TITIK KESETIMBANGAN
Untuk menguji kestabilan ketiga titik kesetimbangan, definisikan fungsi berikut ( )
      + + + = haP N K hr N hr P N f
1 1 ) , ( 1 (3a)
( ) ( ) . 1 1 ) , ( 2 mh P dh hcaN P N f + + − = (3b)
Matriks Jacobi dari sistem (2) yang dilinierisasi di sekitar titik kesetimbangan ( N *, P *) adalah sebagai berikut
       

         =
P P N f N P N f P P N f N P N f
J P N *) *, ( *) *, ( *) *, ( *) *, ( 2 2 1 1 *) *, ( (4)
dengan
( ) 2 1 * * 1 * 1 ) 1 ( *) *, (       + + + + =   haP N K hr haP hr N P N f 2 1 * * 1 * ) 1 ( *) *, (
      + + + − =   haP N K hr haN hr P P N f ( ) mh hcaP N P N f + =   1 * *) *, ( 2
http://dx.doi.org/10.14710/jfma.v2i1.24
JOURNAL OF FUNDAMENTAL MATHEMATICS AND APPLICATIONS (JFMA) VOL. 2 NO. 1 (Jun 2019) Available online at www.jfma.math.fsm.undip.ac.id
( ) mh dh hcaN P P N f + + − =   1 1 * *) *, ( 2 .
Titik kesetimbangan ( N *, P *) akan stabil jika semua nilai eigen ( i  , i =1, 2) matriks Jacobi pada persamaan (4) memenuhi
1  i  .
Matriks Jacobi (4) untuk TE 1 (0, 0) dapat dituliskan sebagai berikut
        + + − + = mh dh hr J TE 1 1 0 0 1 1 . (5)
Nilai eigen dari matriks Jacobi (5) adalah
1 1 + = hr  dan mh dh + − = 1 1 2  . Karena 1 1   , maka TE 1 tidak stabil.
Selanjutnya, matriks Jacobi (4) untuk TE 2 ( K , 0) dapat dituliskan sebagai berikut
          + + − + − + = mh dh hcaK hr haK hr J TE 1 1 0 1 1 1 2 . (6)
Nilai eigen dari matriks Jacobi (6) adalah
1 1 1 + = hr  dan mh dh hcaK + + − = 1 1 2  . Karena (1 +
hr ) >1 maka 1 0 1    . Dengan demikian, TE2 akan stabil jika memenuhi kedua syarat berikut:
(i) h ( caK – d ) > -2 + hr (ii) caK – d < r .
Sementara itu, untuk TE 3 , matriks Jacobi dapat ditulis sebagai berikut
                + + +       + − + + +       + − + =
mh dm mh caK m d hcr hr m d h hr caK m d hr J TE 1 1 1 1 1 ) ( 1 1 1 3 . (7)
Persamaan karakteristik matriks Jacobi (7) adalah
            +       + − +       + + +             +       + − + + + + − hr caK m d hr mh dm hr caK m d hr mh dm 1 1 1
1 1 1 1 1 1 1 2   0 ) 1 )( 1 ( 1 ) ( 2 = + +       + − + − mh hr caK m d m d cr h . (8)
Dari persamaan karakteristik (8), misalkan,
JOURNAL OF FUNDAMENTAL MATHEMATICS AND APPLICATIONS (JFMA) VOL. 2 NO. 1 (Jun 2019) Available online at www.jfma.math.fsm.undip.ac.id
            +       + − + + + + = hr caK m d hr mh dm A 1 1 1 1 1
(9)
) 1 )( 1 ( 1 ) ( 1 1 ) 1 ( 2 mh hr caK m d m d cr h caK m d hr dm B + +       + − + −             + − + + = . (10)
Untuk menentukan nilai eigen dari persamaan (8), dapat digunakan Lemma berikut [7, 8].
Lemma 1 . Akar-akar persamaan kuadrat 0 2 = + − B A   memenuhi 2 , 1 , 1 =  i i  jika dan hanya jika syarat-syarat berikut terpenuhi
(i) 1 + A + B > 0 (ii) 1 – A + B > 0 (iii) B < 1 Dari syarat kepositifan titik kesetimbangan TE 3 , ( d + m ) < caK , sehingga nilai A > 0. Sementara itu, TE 3 akan  stabil  jika  memenuhi  ketiga  syarat  pada  Lemma  1  yaitu:  (a)
) 1 )( 1 ( 1 ) ( 1 1 ) 1 ( 0 2 mh hr caK m d m d cr h caK m d hr dm + +        + − + −             + − + + 
(b)       + + +        + + + + caK dm hc m d hm m d caK m d h m ) ( ) ( 1 ) ( .
Gambar 1. Penyelesaian ( N n , P n ) untuk berbagai nilai awal
TE 2 (0,25, 0) TE 3 (0,2, 0,012)

## JOURNAL OF FUNDAMENTAL MATHEMATICS
AND APPLICATIONS (JFMA) VOL. 2 NO. 1 (Jun 2019) Available online at www.jfma.math.fsm.undip.ac.id
Stabil di TE 2 Stabil di TE 3 (a) Stabil di TE 2 Stabil di TE 3 (b) Stabil di TE 2 Stabil di TE 3 (c) Gambar 2 . Perbandingan simulasi (a) skema beda hingga tak-standar dengan , ) ( h h =  (b)
skema beda hingga tak-standar dengan r e h rh 1 ) ( − =  , dan (c) metode Runge-kutta orde 4

## IV. SIMULASI
Simulasi  dilakukan  untuk  melihat  kekonvergenan  sistem  pada  berbagai  nilai  awal. Pada simulasi ini, diambil nilai-nilai parameter r = 0,03, a = 0,5, c = 0,4, d = 0,01, m = 0,03, h = 0,2, K = 0,25. Dengan nilai-nilai parameter tersebut, untuk P 0 ≠ 0, penyelesaian konvergen ke TE 3 . Sementara itu, untuk P 0 = 0, penyelesaian konvergen ke TE 2 .
http://dx.doi.org/10.14710/jfma.v2i1.24
JOURNAL OF FUNDAMENTAL MATHEMATICS AND APPLICATIONS (JFMA) VOL. 2 NO. 1 (Jun 2019) Available online at www.jfma.math.fsm.undip.ac.id
Selanjutnya, simulasi dilakukan untuk melihat perbandingan skema beda hingga tak standar dengan metode Runge-kutta orde 4. Pada bagian ini, ukuran langkah ditentukan h = 0,2.  Hasil  simulasi  ditunjukkan  pada  Gambar  2.  Pada  Gambar  2,  terlihat  untuk  skema  beda hingga tak-standar baik  dengan , ) ( h h =  maupun dengan r e h rh 1 ) ( − =  tidak menunjukkan perbedaan  yang  signifikan.  Lain  halnya  dengan  metode  Runge  Kutta  orde  4,  terlihat  untuk iterasi n yang lebih kecil, grafik lebih cepat menuju TE 2 (0,25, 0) dan TE 3 (0,2, 0,012).

## V. KESIMPULAN
Skema beda hingga tak standar memberikan alternatif solusi untuk sistem persamaan diferensial  non  linier.  Sistem  tersebut  didiskritisasi  menjadi  sebuah  sistem  persamaan  beda. Model predator-prey dengan pemanenan yang telah dibentuk menjadi sistem persamaan beda dengan  menggunakan  skema  beda  hingga  tak-standar  memberikan  perilaku  dinamik  yang konsisten  dengan  sistem  kontinunya.  Pemilihan  fungsi h h = ) (  dan r e h rh 1 ) ( − =  tidak memberikan perbedaan yang signifikan terhadap perilaku sistem di sekitar titik kesetimbangan. Jika  dibandingkan  dengan  metode  Runge-kutta,  metode  Runge-kutta  orde  4  membutuhkan jumlah iterasi yang lebih sedikit agar sistem konvergen. Akan tetapi, untuk tingkat keakuratan tentu saja harus dibuktikan dengan membandingkan kedua metode terhadap data di lapangan.

## DAFTAR PUSTAKA
[1]   Edwards, C. H. dan Penney, D. E. 2000 : Elementary Differential Equations with Boundary Value Problems Fourth Edition . Prentice-Hall. New Jersey.
[2]  Li, B., Liu, S., Cui, J., dan Li, J. 2016. A Simple Predator-Prey Population Model with Rich Dynamics. Applied Sciences. 6(115). 150 – 168.
[3]  Suryanto, A. 2012. Skema Beda Hingga Tak-Standar untuk Model Epidemi dengan Laju Penularan  Tersaturasi  yang  Dimodifikasi ,  Konferensi  Nasional  Matematika  XVI, Jatinangor.
[4]  Allen, L.J.S. 2006. An Introduction to Mathematical Biology . Pearson Education. United States.
[5]  Mickens,  R.  E.  2000. Application  of  Nonstandard  Finite  Difference  Schemes ,  World Scientific. Publishing Co Pte. Ltd.
[6]  Mickens,  R.E.  2007.  Calculation  of  Denominator  Functions  for  Nonstandard  Finite Difference Schemes for Differential Equations Satisfying a Positivity Condition. Numer. Methods Partial Differential Equations . 23(3).672 – 691.
[7]  Brauer, F. And Castillo-Chavez, C. 2001. Mathematical Models in Population Biology and Epidemiology . Springer. New York.
[8]  Elaydi, S. 2005. An Introduction to Difference Equations. 3rd ed. Springer. New York.
JOURNAL OF FUNDAMENTAL MATHEMATICS AND APPLICATIONS (JFMA) VOL. 2 NO. 1 (Jun 2019) Available online at www.jfma.math.fsm.undip.ac.id