## Analisis Model Penipisan Sumber Daya Hutan Bakau akibat Pembukaan Tambak Ikan dengan Waktu Tunda
Analysis of Mangrove Forest Resource Depletion Models due to The Opening of Fish Pond Land with Time Delay
Ummi Usrotus S. Wakhidah *) , Kartika Nugraheni, Winarni Program Studi Matematika, Jurusan Matematika dan Teknologi Informasi, Institut Teknologi Kalimantan
* E-mail: ummyusrotussw@gmail.com

## ABSTRACT
Mangrove  forest  is  an  ecosystem  with  many  resources  and  high  biological  diversity  that  contain many species of animal, such as proboscis monkeys. One of the problems that keep happening on mangrove forest is the opening of fish pond. Opening of fish pond can damage mangrove resource by the degradation of the mangrove habitat. This study aims to analyze the depletion stability model of mangrove forest resources with time delay, and to show effect of time delay and the presence of change in the equilibrium point (Hopf bifurcation) on depletion stability model of mangrove forest resources due to fish pond openings. The depletion stability model of mangrove forest resources can be  modeled  as  a  system  of  nonlinear  differential  equations.  In  the  numerical  simulation  results, critical value of delay (𝜏 0 ) is 7.05 and the transverse conditions are not fulfilled. This mean there is no change in stability and the delay time parameter (𝜏) does not affect the stability of the system. System  will remain stable  when the  mangrove  forests, proboscis  monkeys, and  fish ponds are in equilibrium. In other words, mangrove forests, proboscis monkeys, and fish ponds can coexist even though time delays, and the analysis using the Hopf bifurcation cannot be carried out.
Keywords: Hopf bifurcation, mangrove forest, time delay.

## PENDAHULUAN
Indonesia  merupakan  gugusan  pulau yang memiliki  luas  hutan  bakau  terbesar  di  dunia (Onrizal,  2010).  Hutan  bakau  didominasi  oleh spesies  pohon  bakau  diantaranya  tanjung, bogem,  bakau  dan  api-api  yang  merupakan populasi  vegetasi  pesisir  tropis.  Jenis  spesies pohon  tersebut  dimanfaatkan  oleh  biota  laut untuk tumbuh dan berkembang di daerah pasang surut  pantai  berlumpur.  Dilihat  dari  segi ekosistem akuatik, hutan bakau dikenal sebagai tempat  mengasuh  dan  membesarkan  berbagai jenis  hewan  akuatik  seperti  kepiting,  udang, ikan dan kerang-kerangan (Majid et al ., 2016). Selain  berfungsi  sebagai  tempat  tumbuh kembang biota laut, hutan bakau juga memiliki peran penting untuk mencegah terjadinya intrusi air laut ke daratan, erosi, dan abrasi pantai yaitu dengan mengendapkan lumpur di akar pohon.
Saat  ini  keberadaan  hutan  bakau  di Indonesia  semakin  berkurang.  Jumlah  luasan hutan  bakau  dari  tahun  ke  tahun  semakin menurun mengindikasikan bahwa terjadi degradasi  hutan  bakau  yang  cukup  tinggi. Adanya  penyusutan  luasan  hutan  bakau  salah satunya  disebabkan  oleh  aktivitas  konversi hutan bakau menjadi tambak (Riwayati, 2014). Pengembangan budidaya tambak dengan
mengeksploitasi  bakau  di  Indonesia  sudah dimulai  sejak  tahun  1800,  dan  berdasarkan penelitian  Ilman et  al. (2016)  dalam  kurun waktu  dua  dekade  kedepan  penyebab  utama kerusakan  ekosistem  bakau  disebabkan  oleh aktivitas  konversi  hutan  bakau  menjadi  lahan tambak  jika  dibandingkan  dengan  aktivitas lainnya seperti pembangunan pesisir, penebangan dan bencana alam.
Kemudian  sebagian  besar  wilayah  hutan bakau  menjadi  habitat  salah  satu  primata  yaitu bekantan  (Kartono et  al. ,  2008).  Bekantan membutuhkan  pohon  di  hutan  bakau  sebagai sarana  pergerakan,  pohon  tidur,  dan  sumber pakan.  Bagian  tumbuhan  yang  digemari  oleh bekantan  untuk  di  konsumsi  adalah  bagian umbut,  pucuk,  bunga,  buah,  daun  muda  dan batang  muda (Syarif  & Maulana, 2019). Jenis- jenis  tumbuhan  yang  menjadi  konsumsi bekantan  di  hutan  bakau  meliputi  tumbuhan bakau dan tumbuhan liar. Salah satu jenis pohon bakau  yang  teridentifikasi  sebagai  sumber makanan  paling  umum  bagi  bekantan  ialah Sonneratia  caseolaris (Noorhidayah, et  al , 2007).  Berdasarkan  penelitian  Akbar et  al.
(2002)  di  Cagar  Alam  Pulau  Kaget  bekantan melakukan  penggundulan  sebagian  pucuk  dan memakan  ranting S.  caseolaris sehingga mengganggu  pertumbuhan  anakan  bakau.
Secara  tidak  langsung,  keberadaan  bekantan dapat menghambat pertumbuhan bakau.
Pengembangan  ilmu  di  bidang  matematika turut memberikan peran penting sebagai kajian untuk  mencegah  masalah  penipisan  sumber daya hutan bakau. Penelitian terhadap penipisan hutan bakau telah dilakukan oleh Nugraheni et al. (2019).  Pada  tahun  2019,  Nugraheni et  al. membuat model matematis yang diformulasikan untuk  memeriksa  masalah  penipisan  sumber daya hutan bakau karena populasi bekantan dan pembukaan lahan untuk tambak ikan di kawasan konservasi.  Pada  model  tersebut  penipisan sumber  daya  hutan  bakau  bergantung  pada waktu sekarang.
Berdasarkan  penelitian  yang  dikemukakan oleh  Nugraheni et  al. (2019)  dianggap  perlu dilakukan  penelitian  lebih  lanjut  mengenai penipisan sumber daya hutan bakau menggunakan  waktu  tunda.  Waktu  tunda ditambahkan  untuk  menyatakan  selang  waktu yang  dibutuhkan  untuk  mengkonversi  hutan bakau  menjadi  lahan  tambak  ikan  dari  tambak pertama  ke tambak  kedua, dan seterusnya. Hal ini didasarkan pada proses fenomena lingkungan  yang  terjadi  secara  real  dalam berbagai permasalahan nyata sering melibatkan waktu  tunda.  Pada  penelitian  ini  dianalisis kestabilan  model penipisan  sumber daya  hutan bakau  akibat  pembukaan  lahan  tambak  ikan dengan  waktu  tunda,  serta  pengaruh  waktu tunda  dan  keberadaan  bifurkasi  Hopf  terhadap kestabilan  model  tersebut.  Hasil  dari  analisis kemudian  disimulasikan  secara  numerik  untuk memperlihatkan dinamika model yang disajikan dalam bentuk grafik.

## METODE

## Model Matematika
Model matematika dari permasalahan ini melibatkan 3  variabel,  yaitu  konsentrasi  sumber  daya  hutan bakau (𝑀(𝑡)) ,  populasi  bekantan (𝑁(𝑡)) ,  dan konsentrasi tambak ikan (𝑃(𝑡)) yang direkonstruksi dari  model  Nugraheni et  al .  (2019)  dengan menambahkan  waktu  tunda.  Penambahan  waktu tunda terdapat pada fungsi konsentrasi lahan tambak ikan  persamaan  konsentrasi  sumber  daya    hutan bakau (𝑀(𝑡)) dan  pada  fungsi  sumber  daya  hutan bakau  persamaan  konsentrasi  lahan  tambak  ikan (𝑃(𝑡)) .
Bakau  tumbuh  alami  tanpa  adanya  persaingan dinyatakan dengan 𝑠𝑀 . Laju pertumbuhan konsentrasi  sumber  daya  hutan  bakau  akan  terbatas oleh  luas  lahan  yang  tersedia  sehingga  kapasitas maksimum habitat hutan bakau dinyatakan dengan 𝐿 dan b merupakan respon konsentrasi bakau terhadap tekanan  lingkungan  sebagaimana  ditunjukkan  oleh
rasio 𝑀 𝐿 .  Hutan  bakau  dimanfaatkan  oleh  bekantan sebagai  tempat  keberlangsungan  hidup,  sehingga konsentrasi  hutan  bakau  akan  berkurang  seiring dengan  pertumbuhan  populasi  bekantan  yang dinyatakan  dengan 𝛼𝑀𝑁 .  Tingkat  penipisan  hutan bakau  akibat  pembukaan  tambak  ikan  dinyatakan dengan 𝑒𝑀 2 𝑃 .
Pertumbuhan  bekantan  secara  alami  dan  dapat dipengaruhi  oleh  eksternal  dinyatakan  dengan 𝑟𝑁 . Laju pertumbuhan konsentrasi bekantan akan terbatas oleh  luas  lahan  yang  tersedia  sehingga  kapasitas maksimum  populasi  bekantan  dinyatakan  dengan 𝐾 dan 𝑐 merupakan respon populasi bekantan terhadap tekanan  lingkungan  sebagaimana  ditunjukkan  oleh rasio 𝑁 𝐾 . Bertambahnya konsentrasi populasi bekantan disebabkan  oleh  tingkat  pertumbuhan  bekantan karena  adanya  interaksi  dengan  sumber  daya  hutan bakau yang dinyatakan dengan 𝛽𝛼𝑀𝑁 .
Berkurangnya konsentrasi populasi bekantan disebabkan  oleh  tingkat  kematian  alami  bekantan yang dinyatakan dengan 𝜇𝑁 .
Konsentrasi  pembukaan  lahan  tambak  ikan dipengaruhi  oleh  berkurangnya  sumber  daya  hutan bakau  akibat  proses  penebangan  hutan  yang dinyatakan  sebagai 𝛾𝑀 .  Konsentrasi  lahan  tambak ikan  dapat  berkurang  apabila  lahan  tambak mengalami  penurunan  kualitas  lahan  secara  alami yang dinyatakan dengan 𝑑𝑃 .
Berdasarkan  penjelasan  tersebut  dapat  disajikan dalam model matematika sebagai berikut:
𝑑𝑀 𝑑𝑡 = 𝑠𝑀 (1 − 𝑏 𝐿 𝑀) − 𝛼𝑀𝑁 − 𝑒𝑀 2 𝑃(𝑡 − 𝜏),
𝑑𝑁 𝑑𝑡 = 𝑟𝑁 (1 − 𝑐 𝐾 𝑁) + 𝛽𝛼𝑀𝑁 − 𝜇𝑁, (1) 𝑑𝑃 𝑑𝑡 = 𝛾𝑀(𝑡 − 𝜏) − 𝑑𝑃.
dengan    nilai    awal 𝑀(0) ≥ 0,   𝑁(0) ≥ 0, dan 𝑃(0) ≥ 0 .

## Waktu Tunda dan Bifurkasi Hopf
Model  Persamaan  (1)  merupakan  model  persamaan diferensial  nonlinier  dengan 𝜏 > 0 menyatakan selang  waktu  yang  dibutuhkan  untuk  mengkonversi hutan  bakau  menjadi  tambak  ikan  dari  tambak pertama  ke  tambak  kedua,  dan  seterusnya.  Secara matematis  selang  waktu  atau  penundaan  pada konversi  hutan  bakau  disebut  sebagai  waktu  tunda ( time delay ). Pada sistem dinamis, kestabilan sistem dapat  dipengaruhi  oleh  waktu  tunda  (Montolalu  & Mongi,  2017).  Titik  kesetimbangan  yang  pada awalnya  bersifat  stabil  lalu  diberikan  waktu  tunda dapat  berubah  menjadi  tidak  stabil,  begitupun sebaliknya. Perubahan kestabilan dari sistem apabila terjadi  di  titik  kesetimbangan  akibat  adanya perubahan parameter disebut bifurkasi lokal (Puspita, 2016). Bifurkasi lokal terdiri dari bifurkasi fold dan bifurkasi  Hopf.  Pada  penelitian  ini  yang  dibahas hanya bifurkasi Hopf. Bifurkasi Hopf terjadi jika titik kesetimbangan  mempunyai  sepasang  akar  imajiner murni, yaitu 𝜆(𝜏 ∗ ) = 𝑝(𝜏 ∗ ) + 𝑖𝑞(𝜏 ∗ ) dengan
𝑝(𝜏 ∗ ) = 0 , 𝑞(𝜏 ∗ ) ≠ 0 dan memenuhi syarat  kondisi transversal 𝑑 = 𝑑 𝜏 [𝑅𝑒 𝜆(𝜏 ∗ )] ≠ 0 (Tu, 1994).

## Persamaan Diferensial Waktu Tunda
Persamaan  diferensial  tundaan  dapat  ditulis  dalam bentuk 𝑥̇(𝑡) = 𝑓(𝑥(𝑡), 𝑥(𝑡 − 𝜏)), (2) dengan 𝑥̇ = turunan 𝑥 terhadap 𝑡 , 𝜏 = waktu tunda, 𝜏 ≤ 0 (Kuang, 1993). Roussel  (2019)  mengatakan  variabel  yang tertunda ditunjukkan oleh nilai subskrip keterlambatan (𝑥 𝜏 ≡ 𝑥(𝑡 − 𝜏)) .  Analisis  kestabilan linear  pada  persamaan  diferensial  waktu  tunda dengan  memisalkan 𝜀𝑥(𝑡) merupakan  perpindahan dari  titik  kesetimbangan  yang  diasumsikan  sangat kecil  sekali  nilainya  dalam  interval  terbuka [𝑡 0 − 𝜏 𝑚𝑎𝑥 , ∞) , sedemikian sehingga diperoleh
𝑥 = 𝑥 ∗ + 𝜀𝑥
(3)
dan
𝑥̇ = 𝜀𝑥̇ = 𝑓(𝑥 ∗ + 𝜀𝑥, 𝑥 ∗ + 𝜀𝑥 𝜏 1 , … , 𝑥 ∗ + 𝑥 𝜏 𝑛 , (4) masing-masing jumlah 𝜀𝑥, 𝜀𝑥 𝜏 1 , … , 𝑥 𝜏 𝑛 adalah kecil.
Persamaan  diferensial  titik  kesetimbangan  dapat dilinearisasi menggunakan metode deret Taylor biasa
𝜀𝑥̇ ≈ J 0 𝜀𝑥 + J 𝜏 1 𝜀𝑥 𝜏 1 + ⋯ + J 𝜏 𝑛 𝜀𝑥 𝜏 𝑛 . (5)
Persamaan  diferensial  waktu  tunda  linear  memiliki solusi eksponensial, dapat ditulis
𝜀𝑥(𝑡) = 𝐴𝜀 𝜆𝑡 . (6)
Kemudian  Persamaan (6) disubstitusikan ke Persamaan (5), diperoleh
𝜆𝐴 = (J 0 + 𝑒 −𝜆𝜏 1 J 𝜏 1 + ⋯ + 𝑒 −𝜆𝜏 𝑛 J 𝜏 𝑛 )𝐴. (7)
Oleh karena 𝐴 ≠ 0 , maka diperoleh
|J 0 + 𝑒 −𝜆𝜏 1 J 𝜏 1 + ⋯ + 𝑒 −𝜆𝜏 𝑛 J 𝜏 𝑛 − 𝜆𝐼| = 0, (8)
dengan
𝐽 0 , J 𝜏 1 , … , J 𝜏 𝑛 = matriks 𝑛 × 𝑛 , 𝐴 = parameter, 𝜆 = akar imajiner murni, 𝐼 = matriks identitas.
Persamaan  (8)  disebut  persamaan  karakteristik  titik kesetimbangan.

## Eksistensi Nilai Waktu Tunda Kritis
Forde & Nelson (2004) mengatakan titik kesetimbangan  akan  bersifat  stabil  ketika  akar  dari persamaan  karakteristik  memiliki  bagian  real  yang bernilai negatif pada saat waktu tunda tidak ada (𝜏 = 0) .  Namun,  saat 𝜏 ≠ 0 akar  karakteristik  akan mengalami  perubahan  kestabilan.  Terdapat  nilai kritis  tundaan  di  mana  perubahan  kestabilan  terjadi ketika Persamaan (9) memiliki akar imajiner murni.
𝐺(𝜆, 𝜏) ≡ 𝐺 1 (𝜆) + 𝐺 2 (𝜆)𝑒 −𝜆𝜏 = 0, (9)
dengan 𝐺 1 dan 𝐺 2 = polinomial dalam 𝜆 , 𝜏 = panjang waktu tunda diskrit. Misalkan  akar  Persamaan  (9)  adalah 𝜆 = 𝑖𝜔, 𝜔 > 0, 𝜔 ∈ ℝ , sehingga Persamaan (9) menjadi 𝐺 1 (𝑖𝜔) + 𝐺 2 (𝑖𝜔)𝑒 −(𝑖𝜔)𝜏 = 0, (10)
jika  bagian  polinomial  dipisahkan  menjadi  bagian real  dan  imajiner,  dan  suku  eksponensial  diubah dalam bentuk trigonometri diperoleh
𝑅 1 (𝜔) + 𝑖𝑄(𝜔) + (𝑅 2 (𝜔) + 𝑖𝑄 2 (𝜔))
(cos(𝜔𝜏) − 𝑖 sin(𝜔𝜏)) = 0, (11)
karena murni imajiner, 𝑅 1 dan 𝑅 2 adalah polinomial genap  dari 𝜔 . 𝑄 1 dan 𝑄 2 adalah  polinomial  ganjil. Agar Persamaan (11) terpenuhi, maka  kedua bagian real dan imajiner harus sama dengan nol, diperoleh
−𝑅 1 (𝜔) = 𝑅 2 (𝜔) cos(𝜔𝜏) + 𝑄 2 (𝜔)
sin(𝜔𝜏) = 0,
(12)
𝑄 1 (𝜔) = 𝑅 2 (𝜔) sin(𝜔𝜏) − 𝑄 2 (𝜔) cos(𝜔𝜏) = 0.
Persamaan  (12)  dikuadratkan  dan  dijumlahkan sehingga diperoleh
𝑅 1 (𝜔) 2 + 𝑄 1 (𝜔) 2 = 𝑅 2 (𝜔) 2 + 𝑄 2 (𝜔) 2 , (13)
dari  Persamaan  (13)  dapat  dilihat  dua  hal.  Pertama, bentuk  trigonometri  menghilang  dan  waktu  tunda 𝜏 telah  dieliminasi.  Kedua,  Persamaan  (13)  ini  sama dengan polinomial berderajat dua.
Misalkan  ditemukan 𝜔 ∗ memenuhi  Persamaan (13) serta Persamaan (12) dikali dengan 1 𝑐𝑜𝑠 𝜔 ∗ 𝜏 dan diubah ke dalam bentuk sudut 𝑠𝑒𝑐 sehingga diperoleh
𝑠𝑒𝑐 𝜔 ∗ 𝜏 = 𝑅 2 (𝜔 ∗ ) + 𝑄 2 (𝜔 ∗ ) 𝑡𝑎𝑛 𝜔 ∗ 𝜏 −𝑅 1 (𝜔 ∗ ) , (14a)
𝑠𝑒𝑐 𝜔 ∗ 𝜏 = 𝑅 2 (𝜔 ∗ ) 𝑡𝑎𝑛 𝜔 ∗ 𝜏 − 𝑄 2 (𝜔 ∗ ) 𝑄 1 (𝜔 ∗ ) . (14b)
Persamaan (14a) disubstitusikan ke Persamaan (14b) diperoleh
𝜏 𝑛 ∗ = 1 𝜔 ∗ 𝑎𝑟𝑐 tan { 𝑄 2 (𝜔 ∗ )𝑅 1 (𝜔 ∗ ) − 𝑅 2 (𝜔 ∗ )𝑄 1 (𝜔 ∗ ) 𝑄 1 (𝜔 ∗ )𝑄 2 (𝜔 ∗ ) + 𝑅 1 (𝜔 ∗ )𝑅 2 (𝜔 ∗ ) } + 2𝑛𝜋 𝜔 ∗ , 𝑛 = 0,1,2, …. (15)
Hasil pada Persamaan (15) menunjukkan terdapat tak  hingga 𝜏 𝑛 ∗ yang  memenuhi  Persamaan  (9).  Oleh karena itu, pilih 𝑛 = 0 agar efektif. Jadi nilai 𝜏 𝑛 ∗ = 𝜏 0 yang merupakan titik kritis.

## HASIL DAN PEMBAHASAN
Penentuan  Titik  Kesetimbangan  Model tanpa Waktu Tunda
Titik kesetimbangan dari model Persamaan (1) dilambangkan  dengan 𝐸 𝑖 .  Titik  kesetimbangan diperoleh pada saat kondisi berikut
𝑑𝑀 𝑑𝑡 = 𝑑𝑁 𝑑𝑡 = 𝑑𝑃 𝑑𝑡 = 0,
dan ditentukan saat 𝜏 = 0 atau saat model tanpa waktu tunda, sehingga Persamaan (1) menjadi 𝑑𝑀 𝑑𝑡 = 𝑠𝑀 (1 − 𝑏 𝐿 𝑀) − 𝛼𝑀𝑁 − 𝑒𝑀 2 𝑃(𝑡 − 𝜏) = 0,
𝑑𝑁 𝑑𝑡 = 𝑟𝑁 (1 − 𝑐 𝐾 𝑁) + 𝛽𝛼𝑀𝑁 − 𝜇𝑁                  = 0, 𝑑𝑃 𝑑𝑡 = 𝛾𝑀(𝑡 − 𝜏) − 𝑑𝑃                                           = 0, dengan  menggunakan  perhitungan  aljabar diperoleh enam titik kesetimbangan, yaitu  𝐸 1 (0,0,0),
68                                      Analisis Model Penipisan Sumber Daya …                     (Wakhidah, dkk)
 𝐸 2 (0, 𝐾 𝑟𝑐 (𝑟 − 𝜇), 0),
 𝐸 3,4 (𝑀 1,2 ∗ , 0, 𝛾 𝑑 𝑀 1,2 ∗ ),
 𝐸 5,6 (𝑀 3,4 ∗ , 𝐾 𝑟𝑐 (𝛽𝛼𝑀 3,4 ∗ − 𝜇 + 𝑟), 𝛾 𝑑 𝑀 3,4 ∗ ), dengan nilai
𝑀 1,2 ∗ = − 𝑑 2𝑒𝛾 ( 𝑏 𝐿 𝑠 ± √( 𝑏 𝐿 𝑠) 2 + 4𝑒 𝛾 𝑑 𝑠), 𝑀 3,4 ∗ = − 𝑑 2𝑒𝛾 ( 𝑏 𝐿 𝑠 + 𝛽𝛼 2 𝐾 𝑟𝑐 ± √( 𝑏 𝐿 𝑠 + 𝛽𝛼 2 𝐾 𝑟𝑐 ) 2 + 4𝑒 𝛾 𝑑 (𝑠 + (𝜇 − 𝑟)𝛼 𝐾 𝑟𝑐 ))
(Nugraheni et al ., 2019). Analisis Kestabilan Model dengan Waktu Tunda Analisis kestabilan dilakukan dengan pendekatan  linear  pada  persamaan  diferensial waktu tunda dengan memisalkan 𝜀𝑀 ̅(𝑡) , 𝜀𝑁 ̅(𝑡) , dan 𝜀𝑃̅(𝑡) merupakan  perpindahan  dari  titik kesetimbangan  yang  diasumsikan  sangat  kecil sekali nilainya, sedemikian sehingga diperoleh
𝑀(𝑡) = 𝑀 ∗ + 𝜀𝑀 ̅(𝑡) , 𝑁(𝑡) = 𝑁 ∗ + 𝜀𝑁 ̅(𝑡) , (16)
𝑃(𝑡) = 𝑃 ∗ + 𝜀𝑃̅(𝑡) .
Persamaan  (16)  disubstitusikan  ke  dalam Persamaan  (1)  diperoleh  persamaan  linear sebagai berikut.
[ 𝑀 ̅ ′ (𝑡) 𝑁 ̅ ′ (𝑡) 𝑃̅ ′ (𝑡) ] = [ −2 𝑏 𝐿 𝑆𝑀 ∗ + 𝑆 − 𝛼𝑁 ∗ − 2𝑒𝑀 ∗ 𝑃 ∗ −𝛼𝑀 ∗ 0
𝛽𝛼𝑁 ∗ −2 𝑐 𝐾 𝑟𝑁 ∗ + 𝑟 + 𝛽𝛼𝑀 ∗ − 𝜇 0
0
0 −𝑑 ] [ 𝑀 ̅ 𝑁 ̅ 𝑃̅ ] + [ 0 0 −𝑒𝑀 ∗ 𝑀 ∗ 0 0 0 𝛾 0 0 ] [ 𝑀 ̅ 𝜏 𝑁 ̅ 𝜏 𝑃̅ 𝜏 ] ,
Setelah dilakukan linearisasi, diperoleh persamaan
𝑑𝑈 ⃗⃗ (𝑡) 𝑑𝑡 = 𝐴 2 𝑈 ⃗⃗ (𝑡) + 𝐵 2 𝑈 ⃗⃗ (𝑡 − 𝜏) . (17)
dengan
𝑈 ̅(𝑡) = [ 𝑀(𝑡) 𝑁(𝑡) 𝑃(𝑡) ], 𝐴 2 = [ 𝑎 11 𝑎 12 0 𝑎 21 𝑎 22 0 0 0 𝑎 33 ], 𝐵 2 = [ 0 0 𝑏 13 0 0 0 𝑏 31 0 0 ] .
Entri-entri dari A 2 dan B 2 adalah
𝑎 11 = −2 𝑏 𝐿 𝑠𝑀
∗ + 𝑠 − 𝛼𝑁 ∗ − 2𝑒𝑀 ∗ 𝑃 ∗ , 𝑎 21 = 𝛽𝛼𝑁 ∗ , 𝑎 12 = −𝛼𝑀 ∗ , 𝑎 22 = −2 𝑐 𝐾 𝑟𝑁 ∗ + 𝑟 + 𝛽𝛼𝑀 ∗ − 𝜇 ,
𝑎 33 = −𝑑 , 𝑏 13 = −𝑒𝑀 ∗ 𝑀 ∗ , dan 𝑏 31 = 𝛾 .
Kemudian  dari  Persamaan  (17)  diperoleh penurunan  persamaan  karakteristik  model waktu  tunda  dengan  memisalkan 𝑈 ̅(𝑡) = 𝐶̅𝑒 𝜆𝑡 yaitu det(𝐴 2 + 𝐵 2 𝑒 −𝜆𝜏 − 𝜆𝐼) = 0 ,  sehingga akar  persamaan  karakteristik  tundaan  dari model Persamaan (1) ialah 𝐺(𝜆, 𝜏) = 𝜆 3 + 𝑎 1 𝜆 2 + 𝑎 2 𝜆 + 𝑎 3 +(𝑎 4 𝜆 + 𝑎 5 )(𝑒 −𝜏𝜆 ) 2 , 0 = 𝜆 3 + 𝑎 1 𝜆 2 + 𝑎 2 𝜆 + 𝑎 3 +(𝑎 4 𝜆 + 𝑎 5 )(𝑒 −𝜏𝜆 ) 2 , (18)
dengan 𝑎 1 = −𝑎 11 − 𝑎 22 − 𝑎 33 , 𝑎 2 = 𝑎 11 𝑎 22 + 𝑎 11 𝑎 33 + 𝑎 22 𝑎 33 − 𝑎 12 𝑎 21 , 𝑎 3 = −𝑎 11 𝑎 22 𝑎 33 + 𝑎 12 𝑎 21 𝑎 33 , 𝑎 4 = 𝑏 13 𝑏 31 , dan 𝑎 5 = −𝑎 22 𝑏 13 𝑏 31 .

## Penentuan Nilai Kritis Tundaan
Kriteria  perubahan  kestabilan  pada  Persamaan (18) diperoleh dengan memisalkan 𝜆 = 𝑖𝜔, 𝜔 > 0, 𝜔 ∈ ℝ , sehingga Persamaan (18) menjadi −𝑖𝜔 3 − 𝑎 1 𝜔 2 + 𝑎 2 𝑖𝜔 + 𝑎 3 +𝑎 4 𝑖𝜔 cos 2𝜔𝜏 + 𝑎 4 𝜔 sin 2𝜔𝜏 + 𝑎 5 cos 2𝜔𝜏 − 𝑎 5 𝑖 sin 2𝜔𝜏 = 0.
(19)
dengan 𝑒 −𝑖2𝜔𝜏 = cos 2𝜔𝜏 − 𝑖 sin 2𝜔𝜏 .  Agar Persamaan  (19)  bernilai  nol,  maka  bagian  real dan  bagian  imajiner  harus  bernilai  nol  yang dikerjakan secara terpisah, diperoleh 𝑎 1 𝜔 2 − 𝑎 3 = 𝑎 5 cos 2𝜔𝜏 + 𝑎 4 𝜔 sin 2𝜔𝜏 , (20)
−𝜔 3 + 𝑎 2 𝜔 = 𝑎 5 sin 2𝜔𝜏 − 𝑎 4 𝜔 cos 2𝜔𝜏 .
Mengeliminasi
𝜏
dilakukan dengan menguadratkan  Persamaan  (20),  sehingga diperoleh 𝑎 1 2 𝜔 4
− 2𝑎 1 𝑎 3 𝜔 2 + 𝑎 3 2 = 𝑎 5 cos 2 2𝜔𝜏 +2𝑎 4 𝑎 5 𝜔 (cos 2𝜔𝜏) (sin 2𝜔𝜏) + 𝑎 4 2 𝜔 2 sin 2 2𝜔𝜏 , 𝜔 6 − 2𝑎 2 𝜔 4 +𝑎 2 2 𝜔 2 = 𝑎 5 2 sin 2 2𝜔𝜏 − 2𝑎 4 𝑎 5 𝜔(cos 2𝜔𝜏)
(sin 2𝜔𝜏)
+𝑎 4 2 𝜔 2 cos 2 2𝜔𝜏 . Selanjutnya kedua persamaan tersebut dijumlahkan dan karena cos 2 2𝜔𝜏 + sin 2 2𝜔𝜏 = 1 ,  diperoleh  polinomial  berderajat enam
𝜔 6 + 𝑑 1 𝜔 4 + 𝑑 2 𝜔 2 + 𝑑 3 = 0 , (21)
dengan 𝑑 1 = 𝑎 1 2 − 2𝑎 2 , 𝑑 2 = −2𝑎 1 𝑎 3 + 𝑎 2 2 − 𝑎 4 2 , 𝑑 3 = 𝑎 3 2 − 𝑎 5 2 . Kemudian dimisalkan 𝑣 = 𝜔 2
untuk memudahkan  proses  perhitungan,  sehingga Persamaan  (21)  menjadi ℎ(𝑣) = 𝑣 3 + 𝑑 1 𝑣 2 + 𝑑 2 𝑣 + 𝑑 3 . Nilai akar dari Pesramaan (21) dapat diketahui  dengan  menggunakan  lemma  yang diberikan oleh Ruan dan Wei sebagai berikut:
Lemma 1 Misalkan ∆= 𝑑 1 2 − 3𝑑 2 , a)
Jika 𝑑 3 < 0 ,  maka  Persamaan ℎ(𝑣)
paling  sedikit  memiliki  satu  akar  yang bernilai positif. b) Jika 𝑑
3 ≥ 0 dan ∆ < 0 , maka Persamaan ℎ(𝑣) tidak  memiliki  akar  yang  bernilai
positif.
c) Jika 𝑑
3 ≥ 0 dan ∆ ≥ 0 , maka Persamaan ℎ(𝑣) memiliki  akar-akar  yang  bernilai positif  jika  dan  hanya  jika 𝑣 1 = 1 3 (−𝑑 1 + √∆) > 0 dan ℎ(𝑣 1 ) ≤ 0 .
Persamaan  (21)  memiliki  akar  imajiner murni 𝑖𝜔 , sehingga dengan mensubstitusikan 𝜔 ke Persamaan (21) diperoleh nilai tundaan kritis.
𝜏 𝑘 ∗ = 1 2𝜔 0 𝑎𝑟𝑐 tan { −(𝑎 1 𝑎 4 + 𝑎 5 )𝜔 0 3 + (𝑎 3 𝑎 4 − 𝑎 2 𝑎 5 ) 𝑎 4 𝜔 0 4 + (𝑎 2 𝑎 4 − 𝑎 1 𝑎 5 )𝜔 0 2 + 𝑎 3 𝑎 5 } + 𝑘𝜋 𝜔 0 , 𝑘 = 0,1,2,3, …
Dipilih 𝜏 0 merupakan  nilai  terkecil  dari  tak  terhingga 𝜏 𝑘 ∗ sebagai  nilai  kritis  tundaan  yang dinyatakan dalam bentuk
𝜏 0 = 1 2𝜔 0 𝑎𝑟𝑐 tan { −(𝑎 1 𝑎 4 + 𝑎 5 )𝜔 0 3 + (𝑎 3 𝑎 4 − 𝑎 2 𝑎 5 ) 𝑎 4 𝜔 0 4 + (𝑎 2 𝑎 4 − 𝑎 1 𝑎 5 )𝜔 0 2 + 𝑎 3 𝑎 5 }. (22)
Kondisi Transversal Pada  saat  kritis  tundaan 𝜏 0 ,  persamaan karakteristik  (21)  memiliki  sepasang  akar imajiner,  maka  memungkinkan  terjadinya bifurkasi  atau  perubahan  struktur  orbit.  Syarat terjadinya  bifurkasi  Hopf  adalah  memenuhi kondisi transversal dinyatakan sebagai berikut:
𝑑 𝑑𝜏 𝑅𝑒(𝜆(𝜏))| 𝜏=𝜏0,𝜆=𝑖𝜔0 ≠ 0.
Menurut  Forde  dan  Nelson  (2004)  kondisi tersebut terpenuhi jika dan hanya jika
𝑅 1 (𝑣 ∗ )𝑅 1 ′ 𝑣 ∗ + 𝑄 1 (𝑣 ∗ )𝑄 1 ′ (𝑣 ∗ ) ≠ 𝑅 2 (𝑣 ∗ )𝑅 2 ′ (𝑣 ∗ ) + 𝑄 2 (𝑣 ∗ )𝑄 2 ′ (𝑣 ∗ ). (23)
Misalkan 𝜆(𝜏) = 𝛼(𝜏) ± 𝑖𝜔(𝜏) sebagai    akar Persamaan  (18)   yang      memenuhi 𝛼(𝜏 0 ) = 0 dan 𝑣 ∗ = 𝜔(𝜏 0 ) = 𝜔 0 ,  dari  Persamaan  (20) diperoleh
𝑅 1 (𝜔) = −𝑎 1 𝜔 2 + 𝑎 3 , (24)
𝑅 2 (𝜔) = 𝑎 5 , 𝑄 1 (𝜔) = 𝜔 3 + 𝑎 2 𝜔 , 𝑄 2 (𝜔) = 𝑎 4 𝜔 .
Kemudian Persamaan (24) diturunkan terhadap 𝜔 , sehingga diperoleh
𝑅 1 ′ (𝜔) = −2𝑎 1 𝜔 , (25)
𝑅 2 ′ (𝜔) = 0 ,
𝑄 1 ′(𝜔) = 3𝜔 2 + 𝑎 2 ,
𝑄 2 ′(𝜔) = 𝑎 4 .
Selanjutnya Persamaan (25) disubstitusikan ke Persamaan (23), diperoleh
𝑅 1 (𝜔)𝑅 1 ′ (𝜔) +𝑄 1 (𝜔)𝑄 1 ′ (𝜔) ≠ 𝑅 2 (𝜔)𝑅 2 ′ (𝜔) +𝑄 2 (𝜔)𝑄 2 ′(𝜔) 3𝜔 5 +(2𝑎 1 2 + 4𝑎 2 )𝜔 3 +(𝑎 2 2 − 2𝑎 1 𝑎 3 )𝜔 ≠ 𝑎 4 2 𝜔 . (26)
Jika  Persamaan  (26)  terpenuhi  maka Persamaan Karakteristik (18) memenuhi kondisi  transversal.  Terpenuhinya  kondisi transveral mengakibatkan terjadinya  perubahan kestabilan    titik    kesetimbangan,    yaitu    stabil asimtotik  saat 𝜏 < 𝜏 0 dan tidak stabil saat 𝜏 > 𝜏 0 , atau sebaliknya dengan nilai 𝑅𝑒(𝜆) ≠ 0 .

## Simulasi Numerik
Pada  tahapan  ini  diberikan  beberapa  simulasi numerik yang bersesuaian dengan hasil analisis yang dipaparkan sebelumnya. Simulasi numerik    dilakukan    dengan  mensubstitusikan nilai    parameter    asumsi 𝑏 = 0,3 ; 𝐿 = 90 ; 𝛼 = 0,11 ; 𝑒 = 0,1 ; 𝐾 = 50 ; 𝜇 = 0,44 ; 𝛾 = 0,36 dan menganut nilai  parameter  milik Nugraheni et  al .  (2019)    yaitu 𝑠 = 0,8 ; 𝑟 = 0,3 ; 𝑐 = 0,1 ; 𝛽 = 0,4 ; 𝑑 = 0,5 . Nilai parameter yang  dipilih  memenuhi  salah  satu  syarat
terjadinya  bifurkasi  Hopf,  yaitu  memiliki  akar imajiner murni  pada  Lemma  1  poin  (c)  untuk titik    kesetimbangan 𝐸 6 ,    dengan    nilai 𝑑 1 = 0,7822; 𝑑 2 = −0,0009; 𝑑 3 = 0,000018; ∆= 0,614; dan ℎ(𝑣 1 ) = 0; sehingga  diperoleh 𝜔 0 = 0,0246 . Namun, saat parameter disubstitusikan  ke  Persamaan  (26)  nilai parameter  tersebut  tidak  memenuhi  kondisi transversal
3𝜔 0 5 +(2𝑎 1 2 + 4𝑎 2 )𝜔 0 3 +(𝑎 2 2 − 2𝑎 1 𝑎
3 )𝜔 0 = 0.0033 ≠ 𝑎 4 2 𝜔 = 0.0033 .
Selanjutnya, nilai parameter disubstitusi ke Persamaan  (22)  diperoleh  nilai  kritis  tundaan terkecil 𝜏 0 = 7,05 . Untuk mengamati bifurkasi Hopf akibat variasi nilai parameter 𝜏 , dilakukan pengamatan  dengan  menggunakan  nilai  awal (3,189;  0,5375;  2,29) .  𝜏 < 𝜏 0 Simulasi model dengan waktu tunda 𝜏 < 𝜏 0 ditunjukkan dengan memilih nilai 𝜏 = 2 .
(a) (b)
(c) (d)
Gambar  1.  Solusi  sistem  untuk 𝜏 = 2 ,  (a) pertumbuhan  bakau,  bekantan,  dan tambak  ikan,  (b)  pertumbuhan bakau,  (c)  pertumbuhan  bekantan,
(d) pertumbuhan tambak ikan
Pada  Gambar  1  (a)  titik  berwarna  merah merupakan  nilai  awal  yang  bergerak  ke  arah dalam,  yaitu  menuju  kestabilan  pada  titik kesetimbangan 𝐸 6 (3,1892;  0,5382;  2,2962) dan orbit periodiknya stabil asimtotik.
Kemudian  pada  bidang  solusi  (b),  (c),  dan  (d) tingkat  pertumbuhan  konsentrasi  hutan  bakau, bekantan, dan tambak ikan berosilasi diawal dan semakin  lama  simpangannya  semakin  kecil hingga  mencapai  kondisi  setimbang.  Dengan demikian, sistem dengan nilai parameter 𝜏 < 𝜏 0 bersifat stabil asimtotik pada titik kestabilan 𝐸 6 .
 𝜏 = 𝜏 0 Simulasi model ditunjukkan dengan memilih nilai 𝜏 = 7,05 . Nilai tersebut merupakan nilai yang sama dengan 𝜏 0 = 7,05 .
(a) (b) (c) (d)
Gambar 2. Solusi sistem untuk 𝜏 = 7,05 , (a) pertumbuhan bakau, bekantan, dan tambak ikan, (b) pertumbuhan bakau, (c) pertumbuhan bekantan, (d) pertumbuhan tambak ikan
Pada  Gambar  2  (a),  titik  berwarna  merah merupakan  nilai  awal  yang  bergerak  ke  arah dalam,  yaitu  menuju  kestabilan  pada  titik kesetimbangan 𝐸 6 (3,1892;  0,5382;  2,2962) dan orbit periodiknya stabil asimtotik.
Kemudian  pada  bidang  solusi  (b),  (c),  dan  (d) tingkat  pertumbuhan  konsentrasi  hutan  bakau, bekantan, dan tambak ikan berosilasi diawal dan
semakin  lama  simpangannya  semakin  kecil hingga  mencapai  kondisi  setimbang.  Dengan demikian, sistem dengan nilai parameter 𝜏 = 𝜏 0 bersifat stabil asimtotik pada titik kestabilan 𝐸 6 .  𝜏 > 𝜏 0
Simulasi  model  dengan  waktu  tunda 𝜏 > 𝜏 0 ditunjukkan dengan memilih nilai 𝜏 = 9 .
(a) (b) (c) (d)
Gambar  3.  Solusi  sistem  untuk 𝜏 = 9 ,  (a)  pertumbuhan  bakau,  bekantan,  dan  tambak  ikan,  (b) pertumbuhan bakau, (c) pertumbuhan bekantan, (d) pertumbuhan tambak ikan
Pada  Gambar  3  (a),  titik  berwarna  merah merupakan  nilai  awal  yang  bergerak  ke  arah dalam,  yaitu  menuju  kestabilan  pada  titik kesetimbangan 𝐸 6 (3,1892;  0,5382;  2,2962) dan orbit periodiknya stabil asimtotik.
Kemudian  pada  bidang  solusi  (b),  (c),  dan  (d) tingkat  pertumbuhan  konsentrasi  hutan  bakau, bekantan, dan tambak ikan berosilasi diawal dan semakin  lama  simpangannya  semakin  kecil hingga  mencapai  kondisi  setimbang.  Dengan demikian, sistem dengan nilai parameter 𝜏 > 𝜏 0 bersifat stabil asimtotik pada titik kestabilan 𝐸 6 . Dalam  penelitian  ini,  berdasarkan  hasil analisis  kondisi  transversal  tidak  terpenuhi, sehingga  tidak  terjadi  perubahan  kestabilan. Hasil simulasi menunjukkan waktu tunda tidak berpengaruh  pada  sistem  di  sekitar  titik kesetimbangan 𝐸 6 . Hal tersebut mengindikasikan  bahwa  berbagai  macam perubahan nilai waktu tunda yang disubstitusikan  ke  dalam  model  tidak  akan mempengaruhi  kestabilan  sistem.  Analisis menggunakan  teori  bifurkasi  Hopf  tidak  bisa dilakukan karena tidak terjadi bifurkasi.

## KESIMPULAN
Perubahan kestabilan tidak terjadi di sistem oleh karena nilai kritis tundaan (𝜏 0 ) adalah 7,05 serta kondisi transversal tidak terpenuhi.
Parameter waktu tunda (𝜏) tidak berpengaruh terhadap kestabilan sistem, sistem akan  tetap  bersifat  stabil  pada  saat  konsentrasi dalam  keadaan  setimbang.  Dengan  kata  lain, jika  waktu  tunda  diasumsikan  terjadi  pada pembukaan lahan tambak ikan, ketiga konsentrasi  bisa  hidup  secara  berdampingan. Dengan  demikian,  tidak  terjadi  bifurkasi  Hopf pada model yang digunakan dalam penelitian ini meskipun diperhitungkan adanya  selang  waktu yang dibutuhkan untuk membuka lahan tambak ikan,  sehingga  analisis  menggunakan  bifurkasi Hopf tidak bisa dilakukan.

## DAFTAR PUSTAKA
Akbar  A,  Priyatno  E,  Basaing  HA,  Manaon  & Rizani  F.  2002. Pengaruh  Pemagaran Terhadap  Pertumbuhan  Awal  Permudaan Rambai  (Sonneratia  caesolaris  ( L )  Engl.) di  Cagar  Alam  Pulau  Kaget .  Buletin Teknologi Reboisasi 9: 42-53. Yogyakarta: Pusat Penelitian dan Pengembangan
Bioteknologi  dan  Pemuliaan  Tanaman Hutan.
Forde  J  &  Nelson  P.  2004.  Applications  of
Sturm Sequences to Bifurcation Analysis of
Delay  Differential  Equation Models. Journal  of  Mathematical  Analysis  and Applications . 300 (2): 273-284.
Ilman M, Dargusch P, Dart P & Onrizal. 2016.
A Historical Analysis of The Drivers Ofloss and Degradation of Indonesia’s Mangroves. Land  Use  Policy . 54 (1):448- 459
Kartono  AP,  Ginting  A  &  Santoso  N.  2008.
Karakteristik  Habitat  dan  Wilayah  Jelajah Bekantan  di  Hutan  Mangrove  Desa  Nipah Panjang Kecamatan Batu Ampar Kabupaten Kubu Raya Provinsi Kalimantn Barat. Media Konservasi . 13 (3): 1-6. Kuang  Y.  1993. Delay  Differential  Equations with Applications in Population Dynamics .
Boston: Academic Press.
Majid I, Muhdar MHIA, Rohman F & Syamsuri
I.  2016.  Konservasi  Hutan  Mangrove  di
Pesisir  Pantai  kota  Ternate  Terintegrasi dengan Kurikulum Sekolah. Jurnal
Bioedukasi . 4 (2): 488-496.
Montolalu  CEJC  &  Mongi  CE.  2017.  Analisis Bifurkasi  Hopf  Dan  Kestabilan  Model
Dinamik  Transmisi  Virus  Dengue  Dengan Waktu Tunda Dalam Pengaruhnya Terhadap  Kemunculan  Penyakit  Demam
Berdarah  Dengue.  d’Cartesian: Jurnal
Matematika dan Aplikasi . 6 (2): 86-95. Noorhidayah,  Sidiyasa  K  &  Ma’ruf  A.  2007.
Struktur  dan  Komposisi  Vegetasi  Habitat Bekantan (Nasalis Larvatus Wurmb.) Pada Hutan  Mangrove  di  Bagian  Hilir  Sungai Wain Kalimantan Timur. Jurnal Penelitian Hutan dan Konservasi Alam . 4 (2):107-116. Nugraheni  K,  Burhan  MJI,  Pancahayani  S  & Azka  M.  2019.  Stability  Analysis  of Mangrove  Forest  Resource  Depletion Models Due to The Opening of Fish Pond Land. Journal  of  Physics:  Conference Series . 1277 (1): 1-8.
Onrizal.  2010.  Perubahan  Tutupan  Hutan
Mangrove di Pantai Timur Sumatera Utara Periode 1977-2006. Jurnal Biologi Indonesia . 6 (2): 163-172.
Puspita  M.  2016.  Bentuk  Normal  Bifurkasi
Hopf  Pada  Sistem  Umum  Dua  Dimensi.
Jurnal Matematika UNAND . 5 (3): 15-23. Riwayati.  2014.  Manfaat  dan  Fungsi  Hutan
Mangrove Bagi Kehidupan. Jurnal Keluarga Sehat Sejahtera . 12 (2): 17-23. Roussel  MR.  2019. Nonlinear  Dynamics  A Hands-On  Introductory  Survey .  Amerika: Morgan & Claypool. Ruan S & Wei J. 2001. On The Zeros of A Third
Degree  Exponential  Polynomial  With Applications to A Delayed Model For The
Control  of  Testosterone  Secretion. IMA Journal  of  Mathematics  Applied  in
Medicine and Biology . 18 (1): 41-52. Syarif A & Maulana F. 2019. Keanekaragaman Jenis Tumbuhan Sebagai Alternatif Konsumsi  Bekantan  ( Nasalis  Larvatus )  di Desa  Lawahan  Kabupaten  Tapin. Jurnal Pendidikan Hayati . 5 (4): 180-186.
74                                      Analisis Model Penipisan Sumber Daya …                     (Wakhidah, dkk)