Abstrak —Dalam sistem budidaya ikan, interaksi yang terjadi melibatkan manusia, organisme, makanan, dan kondisi lingkungan yang mendukung kehidupan didalamnya. Peran manusia juga penting dalam memperbanyak, menumbuhkan, serta meningkatkan kualitas organisme akuatik sehingga diperoleh keuntungan yang optimal. Pemodelan matematika menggunakan model prey-predator banyak digunakan dalam mempelajari ekologi akuatik. Model matematika pada prey- predator diperkenalkan pertama kali secara terpisah oleh A.J Lotka dan Vito Volterra pada tahun 1920. Model prey-predator ini telah banyak dikembangkan. Pada penelitian ini, menggabungkan model prey-predator dengan pemanenan proporsional dan respons fungsional Holling tipe II dengan persamaan diferensial bioekonomi. Tujuan penelitian ini adalah menganalisis kestabilan lokal dan global dari model prey- predator yang kemudian diinterpretasikan pada sistem budidaya ikan bandeng, sehingga dari sistem yang stabil ini diperoleh keuntungan ekonomi yang optimal, juga menyelidiki eksistensi bifurkasi dengan rate pemanenan pada prey sebagai parameter bifurkasi. Kestabilan lokal dianalisis dengan mencari nilai eigen dan kestabilan global dianalisis menggunakan fungsi Lyapunov. Serta menyelidiki kemungkinan parameter yang membuat sistem berubah kestabilannya. Dari hasil analisis yang telah dilakukan, dapat ditunjukkan bahwa titik setimbang 𝑬𝑬 𝟏𝟏 ( 𝒙𝒙 ∗ , 𝒚𝒚 ∗ ) stabil asimtotis yang bermakna bahwa populasi prey dan predator tidak pernah habis. Sistem ini stabil juga karena laju pertumbuhan pada prey lebih besar daripada koefisien predasi oleh predator. Keberadaan predator dalam tambak budidaya bandeng dapat menghambat pertumbuhan ikan bandeng karena mengurangi persediaan makanan dan mempersempit tempat hidup dari ikan bandeng. Dari sistem yang stabil ini kita dapat memperoleh keuntungan ekonomi yang optimal. Akan tetapi, pada sistem ini tidak terjadi bifurkasi dibuktikan dengan stabilitas sistem yang tidak berubah meskipun kita mengganti nilai dari parameter rate pemanenan pada prey 𝒆𝒆 𝟏𝟏 sebanyak 3 kali. Dari simulasi model juga diperoleh keuntungan ekonomi yang optimal ketika populasi prey dan predator sama-sama dipanen.
Kata Kunci —Sistem Budidaya Ikan Bandeng, Model Prey- Predator, Kestabilan Lyapunov, Bifurkasi .

## I. PENDAHULUAN
RESIK  merupakan salah satu kabupaten di Indonesia yang berada di wilayah pesisir. Hampir sepertiga dari wilayah Kabupaten Gresik adalah pesisir pantai dan memiliki sekitar 46% total luas tambak di Provinsi Jawa Timur. Pada data Kecamatan Manyar dalam Angka tahun  2021, Kecamatan Manyar memiliki potensi tambak paling luas diatara kecamatan lainnya. Jenis ikan budidaya tambak yang paling banyak di produksi adalah ikan bandeng. Gresik merupakan salah satu pusat produksi ikan bandeng di Indonesia. Namun, sebuah riset oleh Pusat Penelitian Sosial Ekonomi Kelautan dan Perikanan Kementerian Kelautan dan Perikanan menunjukkan bahwa sistem budidaya ikan
bandeng di Kabupaten Gresik justru tidak memerhatikan aspek keberlanjutan secara ekologi, ekonomi, sosial, dan teknologi. Oleh karena itu, penerapan sistem budidaya ikan Bandeng yang tepat diperlukan dalam pengelolaan sumber daya perikanan berkelanjutan, menjaga dari eksploitasi yang berlebihan serta menjaga kelestarian ekosistem [1].
Masalah atau fenomena ekologi yang terjadi ini dapat diselesaikan menggunakan penerapan ilmu matematika yang disebut pemodelan matematika. Pemodelan Matematika merupakan suatu proses merepresentasikan dan menjelaskan permasalahan pada dunia nyata ke dalam pernyataan matematis. Pembentukan model matematika disesuaikan dengan fenomena yang ada kemudian ditransformasikan ke dalam bentuk matematika berupa persamaan diferensial [2]. Pemodelan matematika menggunakan model prey-predator banyak digunakan dalam mempelajari ekologi akuatik. Model matematika pada prey-predator diperkenalkan pertama kali secara terpisah oleh A.J Lotka dan Vito Volterra pada tahun 1920. Model prey-predator adalah model yang digunakan untuk menggambarkan dinamika dari dua spesies yang berinteraksi, satu sebagai predator dan satunya sebagai prey. Dalam perkembangannya, banyak peneliti yang mengembangkan model ini. Dari hasil beberapa penelitian memperlihatkan bahwa model prey-predator dapat menggambarkan dinamika populasi ikan, menjelaskan pengaruh berbagai faktor pada populasi, serta memprediksi perilaku dinamis populasi [3]. Pada 2017, Meng Li dkk mengembangkan model prey-predator dengan menambahkan persamaan diferensial bioekonomi [3]. 2 tahun setelahnya, Yongzhen Pei dkk meneliti tentang model sistem pengelolaan perikanan dengan kebijakan panen selektif yang juga mengembangkan model prey-predator dengan pemanenan selektif yang menggabungkan dua penundaan usia dan ukuran untuk kedua spesies di daerah yang dilindungi dan tidak dilindungi [4].
Pada tugas akhir ini, akan dilakukan penelitian lebih lanjut

## Analisis Kestabilan Lyapunov dan Eksistensi Bifurkasi dari Model Sistem Budidaya Ikan Bandeng di Kecamatan Manyar, Kabupaten Gresik
Sulis Rizkiatul Fitri, Hariyanto, dan Nur Asiyah Departemen Matematika, Institut Teknologi Sepuluh Nopember (ITS) e-mail : hariyanto@matematika.its.ac.id
G
Gambar 1. Model matematika prey-predator dengan pemanenan prey nonlinear dan respons fungsional holling tipe II.
yaitu menyesuaikan model yang sudah ada sebelumnya dengan fenomena yang terjadi pada sistem budidaya ikan bandeng di Kecamatan Manyar, Kabupaten Gresik selanjutnya menganalisis kestabilan lokal menggunakan nilai eigen dan menganalisis kestabilan global menggunakan fungsi lyapunov dan menyelidiki eksistensi bifurkasi hopf dari model sistem budidaya ikan bandeng dengan mempertimbangkan faktor ekonomi dengan meninjau kembali persamaan diferensial bioekonomi. Kemudian dilakukan simulasi kestabilan dan bifurkasi menggunakan software MATLAB R2021a.

## A. Model Prey-Predator
Model matematika prey-predator Lotka Volterra diperkenalkan pertama kali secara terpisah oleh A.J Lotka pada tahun 1926. Model prey-predator adalah model yang digunakan untuk menggambarkan dinamika dari dua spesies yang berinteraksi, satu sebagai predator dan satunya sebagai prey.
Model prey-predator dalam penelitian ini diasumsikan bahwa prey-nya adalah ikan bandeng, dan predator-nya adalah ikan jenis lain yang lain yang mengganggu pertumbuhan ikan bandeng di dalam ekosistem tambak.
Meng Li dkk mengembangkan model prey-predator dengan menambahkan persamaan diferensial bioekonomi dengan pemanenan prey nonlinear dan respons fungsional Holling tipe II seperti pada Gambar 1. Maka diperoleh bentuk model matematika berikut [3]:
⎩ ⎪ ⎨ ⎪ ⎧ 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 = 𝑟𝑟𝑟𝑟 � 1 − 𝑑𝑑 𝑘𝑘 � − 𝑎𝑎 1 𝑑𝑑𝑥𝑥 𝑛𝑛+𝑑𝑑 − 𝑞𝑞𝑞𝑞𝑑𝑑 𝑚𝑚 1 𝑞𝑞+𝑚𝑚 2 𝑑𝑑 𝑑𝑑𝑥𝑥 𝑑𝑑𝑑𝑑 = 𝑠𝑠𝑠𝑠 �− 1 + 𝑎𝑎 2 𝑑𝑑 𝑛𝑛+𝑑𝑑 � , 0 = 𝑞𝑞𝑞𝑞𝑑𝑑 𝑚𝑚 1 𝑞𝑞+𝑚𝑚 2 𝑑𝑑 ( 𝑝𝑝𝑟𝑟 − 𝑐𝑐 ) − 𝑚𝑚 . (1)
dengan : 𝑟𝑟 = 𝑟𝑟 ( 𝑡𝑡 ) : total populasi prey pada saat t, 𝑠𝑠 = 𝑠𝑠 ( 𝑡𝑡 ) : total populasi predator pada saat t, 𝑟𝑟 : laju pertumbuhan prey,
𝑠𝑠 : laju kematian predator tanpa adanya makanan, 𝑘𝑘 : daya dukung prey,
𝑛𝑛 :  konstanta saturasi prey dan predator maksimal secara terpisah 𝑎𝑎 1 = 𝑎𝑎 2 𝑠𝑠 : koefisien predasi,
𝑞𝑞 : koefisien daya tangkap, 𝐸𝐸 : upaya pemanenan pada prey, 𝑚𝑚 1 , 𝑚𝑚 2 : konstanta positif yang sesuai, 𝑝𝑝 : harga per unit prey yang dipanen,
𝑐𝑐 : biaya panen,
𝑚𝑚 : keuntungan ekonomi. semua parameter diasumsikan positif karena pertimbangan biologis.

## B. Pelinieran
Untuk menganalisis kestabilan sistem non linear pada persamaan (4) digunakan analisis kestabilan lokal di sekitar titik setimbangnya. Untuk menganalisis kestabilan lokal tersebut, digunakan Deret Taylor untuk mencari hampiran di sekitar titik setimbang. Misal 𝑬𝑬 ∗ = ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) adalah titik setimbang dari sistem 𝒙𝒙̇ = 𝒇𝒇 𝟏𝟏 ( 𝒙𝒙 , 𝒚𝒚 ) dan 𝒚𝒚̇ = 𝒇𝒇 𝟐𝟐 ( 𝒙𝒙 , 𝒚𝒚 ) . Dengan memisalkan ∆𝑟𝑟 = 𝑟𝑟 − 𝑟𝑟 ∗ , ∆𝑠𝑠 = 𝑠𝑠 − 𝑠𝑠 ∗ diperoleh:
Gambar 3. Diagram alir penelitian.
Gambar  4. Model prey-predator dengan pemanenan prey-predator proposional dan respons fungsional holling tipe II.
Tabel 1.
Nilai pada parameter Parameter Nilai 𝑟𝑟 2,09 𝑠𝑠 0,5 𝑘𝑘 1,5 𝑛𝑛 0,9 𝑎𝑎 0,5 𝑏𝑏 2 𝑞𝑞 0,003 𝑒𝑒 1 1 𝑒𝑒 2 10 𝑝𝑝 1 3000 𝑝𝑝 2 700 𝑐𝑐 14600000
Gambar 2. Grafik bidang solusi populasi prey dan predator.
�𝑟𝑟̇𝑠𝑠̇� = �
𝜕𝜕𝑓𝑓 1 ( 𝑑𝑑 ∗ , 𝑥𝑥 ∗ ) 𝜕𝜕𝑑𝑑 𝜕𝜕𝑓𝑓 1 ( 𝑑𝑑 ∗ , 𝑥𝑥 ∗ ) 𝜕𝜕𝑥𝑥 𝜕𝜕𝑓𝑓 2 ( 𝑑𝑑 ∗ , 𝑥𝑥 ∗ ) 𝜕𝜕𝑑𝑑 𝜕𝜕𝑓𝑓 2 ( 𝑑𝑑 ∗ , 𝑥𝑥 ∗ ) 𝜕𝜕𝑥𝑥 � �∆𝑟𝑟 ∆𝑠𝑠� (2)
Matriks (2) bisa ditulis sebagai
�𝑟𝑟̇𝑠𝑠̇� = 𝐽𝐽�𝑓𝑓 ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) � �∆𝑟𝑟 ∆𝑠𝑠� dengan 𝐽𝐽�𝑓𝑓 ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) � adalah matriks Jacobian dari fungsi 𝑓𝑓 pada titik kesetimbangan 𝑟𝑟 ∗ , 𝑠𝑠 ∗ .

## C. Analisis Kestabilan
Analisis kestabilan sistem dapat diketahui berdasarkan pada kestabilan lokal dan kestabilan global. Kestabilan lokal ini dapat dianalisis berdasarkan nilai eigen yang didapat dari matriks Jacobian pada saat pelinieran di sekitar titik setimbang.
Didefinisikan suatu titik kesetimbangan pada sistem x ̇ = f( 𝐱𝐱 ) dikatakan stabil apabila nilai eigen bagian realnya negatif atau kompleks conjugate dengan bagian realnya negatif, dan tidak stabil apabila nilai eigen bagian realnya positif atau kompleks conjugate bagian realnya positif . Sedangkan untuk kestabilan global bisa dianalisis dengan fungsi lyapunov.

## D. Fungsi Lyapunov
Fungsi Lyapunov adalah salah satu fungsi yang dikonstruksi untuk memeriksa kestabilan global dari suatu sistem nonlinier.
Definisi fungsi Lyapunov adalah suatu fungsi yang memenuhi tiga pernyataan berikut ini [5]:
Diberikan fungsi V: D ∈ R n > R dan x ∗ ∈ D titik kesetimbangan sistem persamaan diferensial nonlinier. Fungsi V disebut fungsi Lyapunov jika memenuhi pernyataan berikut:
1. Fungsi V kontinu dan mempunyai turunan parsial pertama yang kontinu pada D 2. Fungsi V(x) > 0 untuk x ∈ D dengan x ≠ x ∗ dan V ( x ∗ ) = 0 dengan x = x ∗ (dengan titik kesetimbangan x ∗ merupakan titik minimum global).
3. Fungsi V ̇ (x) ≤ 0 untuk setiap x ∈ D.

## E. Bifurkasi Hopf
Bifurkasi Hopf adalah berubahnya jenis kestabilan suatu titik kesetimbangan karena munculnya sepasang nilai eigen yang bernilai imajiner murni.
Bifurkasi yang sesuai dengan keberadaan λ 1 , 2 = ± ω i, ω > 0 , dengan ω adalah bagian imaginer dari nilai eigen yang terkait. Maka bifurkasi yang terjadi disebut bifurkasi Hopf (atau Andronov Hopf) [6].
Dalam konsep bifurkasi dikenal suatu titik yang disebut titik bifurkasi. Secara umum sistem dinamik dapat ditulis sebagai berikut [7]:
𝐯𝐯̇ = F( 𝐯𝐯 , 𝛍𝛍 )
dengan:
𝐯𝐯 = 𝐯𝐯 (x, y) dan 𝛍𝛍 = 𝛍𝛍 (a 0 , b 0 , r 0 , … . ) .
Jika sebuah pasangan (v 0 , μ 0 ) memenuhi kriteria sebagai berikut:
1. F( 𝐯𝐯 𝟎𝟎 , μ 0 ) = 0 ,
2. J( 𝐯𝐯 , μ ) memiliki 2 nilai eigen kompleks λ 1 , 2 = a( 𝐯𝐯 , μ ) ± ib(v, μ ) di sekitar ( 𝐯𝐯 𝟎𝟎 , μ 0 ) ,
3. a( 𝐯𝐯 , μ ) = 0 dan b( 𝐯𝐯 , μ ) ≠ 0 ,
4. Nilai eigen ketiga λ 3 ≠ 0 , Maka titik ( 𝐯𝐯 𝟎𝟎 , μ 0 ) disebut titik bifurkasi Hopf [8].

## II. METODOLOGI PENELITIAN
Tahap-tahap yang dilakukan dalam penelitian ini adalah melakukan studi literatur, kemudian menyusun asumsi dan konstruksi model yang berkaitan dengan asumsi-asumsi terkait model sistem budidaya ikan bandeng, dilanjutkan dengan menganalisis kestabilan lokal  di sekitar titik kesetimbangan menggunakan nilai eigen, lalu menganalisis kestabilan global menggunakan fungsi lyapunov. Selanjutnya dilakukan eksistensi Bifurkasi Hopf, dilanjutkan dengan melakukan simulasi model dan analisis hasil simulasi. Setelah itu, melakukan penarikan kesimpulan dan penulisan tugas akhir. Diagram alir pengerjaan tugas akhir ini ditampilkan pada Gambar 2.

## III. ANALISIS DAN PEMBAHASAN
A. Model Prey-Predator dengan Penambahan Persamaan Diferensial Bioekonomi dengan Pemanenan Prey-Predator Proporsional dan Fungsi Holling tipe II
Dalam konsep biologi, suatu sistem bisa terbentuk karena adanya hubungan timbal balik antara makhluk hidup dengan sekitarnya. Hubungan timbal balik antar makhluk hidup dalam memenuhi makanannya terjadi pada rantai makanan yang terdiri dari prey dan predator. Model matematika prey- predator Lotka Volterra diperkenalkan pertama kali secara terpisah oleh A.J Lotka pada tahun 1926. Model prey- predator adalah model yang digunakan untuk menggambarkan dinamika dari dua spesies yang berinteraksi, satu sebagai predator dan satunya sebagai prey. Pada model prey-predator dalam penelitian ini, diasumsikan bahwa prey - nya adalah ikan bandeng, dan predator -nya adalah spesies lain yang menghambat pertumbuhan ikan bandeng baik berupa ikan jenis lain, tumbuhan, atau makhluk hidup lain yang mengganggu pertumbuhan ikan bandeng di dalam ekosistem tambak. Dalam kehidupan nyata, fenomena yang terjadi pada ilmu biologi berkaitan dengan ilmu ekonomi atau yang biasa disebut bioekonomi. Bioekonomi ini bertujuan untuk mengoptimalkan keuntungan pada usaha budidaya ikan. Model prey-predator yang dipadukan dengan persamaan diferensial bioekonomi ini diharapkan dapat mengoptimalkan keuntungan pada usaha budidaya ikan dengan tetap memperhitungkan aspek keberlanjutan sumberdaya ikan. Dalam tambak budidaya ikan bandeng, bukan hanya terdapat ikan bandeng tetapi ada ikan lain yang memperebutkan makanan yang sama dan kadang terdapat ikan predator dari bandeng. Namun, predator tidak selalu bisa memakan preynya secara langsung. Predator harus mencari prey untuk dimangsa terlebih dahulu. Fenomena ini termasuk tanda dari fungsi respon Holling tipe II. Pada respon Holling tipe II, tingkat konsumsi melambat karena predator menghabiskan waktunya untuk mencari prey.
Populasi prey 𝑟𝑟 mengalami pertumbuhan secara alami dengan laju pertumbuhan sebesar 𝑟𝑟 tetapi populasi prey berkurang sebesar 𝑟𝑟𝑑𝑑 𝑘𝑘 di setiap bertumbuhnya prey karena
keterbatasan daya dukung lingkungannya. Interaksi antara prey dan predator mengakibatkan populasi prey berkurang akibat predasi, sebanding dengan pertumbuhan populasi prey sebesar 𝑎𝑎 𝑛𝑛+𝑑𝑑 . Populasi prey juga berkurang akibat adanya pemanenan secara proporsional sebesar 𝑞𝑞𝑒𝑒 1 𝑟𝑟 . Adanya interaksi antara prey dan predator menyebabkan dampak negatif bagi predator karena terjadinya persaingan antar sesama predator dalam memperebutkan makanan, sehingga meningkatkan laju kematian sebesar 𝑠𝑠 . Interaksi antara prey dan predator juga menyebabkan laju pertumbuhan predator meningkat sebesar 𝑏𝑏 . Populasi predator juga berkurang akibat adanya pemanenan dan penangkapan predator secara proporsional sebesar 𝑞𝑞𝑒𝑒 2 𝑠𝑠 . Misalkan 𝑟𝑟 = 𝑟𝑟 ( 𝑡𝑡 ) adalah total populasi prey pada saat 𝑡𝑡 dan 𝑠𝑠 = 𝑠𝑠 ( 𝑡𝑡 ) mewakili total populasi predator pada saat 𝑡𝑡 , maka model prey-predator dengan pemanenan prey-predator proporsional dan respons fungsional Holling tipe II dinyatakan dalam bentuk persamaan (3) yang ditunjukkan dalam diagram kompartemen  pada Gambar 3.
� 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 = 𝑟𝑟𝑟𝑟 � 1 − 𝑑𝑑 𝑘𝑘 � − 𝑎𝑎𝑑𝑑𝑥𝑥 𝑛𝑛+𝑑𝑑 − 𝑞𝑞𝑒𝑒 1 𝑟𝑟 , 𝑑𝑑𝑥𝑥 𝑑𝑑𝑑𝑑 = 𝑠𝑠 �−𝑠𝑠 + 𝑎𝑎𝑎𝑎𝑑𝑑 𝑛𝑛+𝑑𝑑 � − 𝑞𝑞𝑒𝑒 2 𝑠𝑠 . (3)
Persamaan 3 ditambahkan dengan persamaan diferensial bioekonomi, maka dapat dinyatakan dalam
⎩ ⎨ ⎧ 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 = 𝑟𝑟𝑟𝑟 � 1 − 𝑑𝑑 𝑘𝑘 � − 𝑎𝑎𝑑𝑑𝑥𝑥 𝑛𝑛+𝑑𝑑 − 𝑞𝑞𝑒𝑒 1 𝑟𝑟 , 𝑑𝑑𝑥𝑥 𝑑𝑑𝑑𝑑 = 𝑠𝑠 �−𝑠𝑠 + 𝑎𝑎𝑎𝑎𝑑𝑑 𝑛𝑛+𝑑𝑑 � − 𝑞𝑞𝑒𝑒 2 𝑠𝑠 , 𝑚𝑚 = 𝑒𝑒 1 ( 𝑝𝑝 1 𝑟𝑟 − 𝑐𝑐 ) + 𝑒𝑒 2 ( 𝑝𝑝 2 𝑠𝑠 − 𝑐𝑐 ). (4)
dengan : 𝑟𝑟 = 𝑟𝑟 ( 𝑡𝑡 ) : total populasi prey pada saat 𝑡𝑡 , 𝑠𝑠 = 𝑠𝑠 ( 𝑡𝑡 ) : total populasi predator pada saat 𝑡𝑡 ,
𝑟𝑟 : laju pertumbuhan prey,
𝑠𝑠 : laju kematian predator tanpa adanya makanan,
𝑘𝑘 : daya dukung prey,
𝑛𝑛 :  konstanta saturasi prey dan predator maksimal secara terpisah,
𝑎𝑎 : koefisien predasi,
𝑏𝑏 : laju pertumbuhan predator,
𝑞𝑞 : koefisien daya tangkap prey-predator,
𝑒𝑒 1 , 2 : laju pemanenan pada prey dan predator, 𝑝𝑝 1 : harga per unit prey yang dipanen,
𝑝𝑝 2 : harga per unit predator yang dipanen,
𝑐𝑐 : biaya yang dibutuhkan setiap panen,
𝑚𝑚 : keuntungan ekonomi.

## B. Titik Keseimbangan
Titik kesetimbangan model prey-predator adalah titik yang diperoleh ketika model berada pada keadaan setimbang, keadaan dimana perubahan banyaknya individu dari setiap populasi sepanjang waktu adalah nol. Pada model prey- predator tipe Holling II dengan faktor pemanenan proporsional akan memiliki titik setimbang jika memenuhi:
𝑑𝑑𝑟𝑟 𝑑𝑑𝑠𝑠 = 𝑓𝑓 1 ( 𝑟𝑟 , 𝑠𝑠 ), 𝑑𝑑𝑎𝑎𝑛𝑛 𝑑𝑑𝑠𝑠 𝑑𝑑𝑡𝑡 = 𝑓𝑓 2 ( 𝑟𝑟 , 𝑠𝑠 ).
Dengan ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) adalah titik kesetimbangan, maka 𝑓𝑓 1 ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) = 0 , 𝑓𝑓 2 ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) = 0 .
Pada model ini terdapat dua titik setimbang yang  non negatif, yaitu titik setimbang ketika prey dan predator tidak ada, dan titik setimbang ketika prey dan predator hidup berdampingan.
Titik kesetimbangan ketika prey dan predator tidak ada yaitu 𝑟𝑟 = 0, 𝑠𝑠 = 0 atau ditulis dengan 𝐸𝐸 0 ( 𝑟𝑟 0 , 𝑠𝑠 0 ) = (0,0) . Titik kesetimbangan ketika kedua populasi hidup berdampingan dinotasikan dengan 𝐸𝐸 1 ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) = ( 𝑟𝑟 1 , 𝑠𝑠 1 ) . Dengan menggunakan syarat 𝑓𝑓 1 ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) = 0, 𝑓𝑓 2 ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) = 0, 𝑟𝑟 ≠ 0, 𝑠𝑠 ≠ 0 , diperoleh 𝐸𝐸 1 ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) = � 𝑛𝑛 ( 𝑠𝑠+𝑞𝑞𝑒𝑒 2 ) 𝑎𝑎𝑎𝑎−𝑠𝑠−𝑞𝑞𝑒𝑒 2 , 𝑎𝑎𝑛𝑛 ( 𝑎𝑎𝑎𝑎𝑘𝑘 ( 𝑟𝑟−𝑞𝑞𝑒𝑒 1 ) −𝑘𝑘𝑟𝑟 ( 𝑠𝑠+𝑞𝑞𝑒𝑒 2 ) +𝑘𝑘𝑞𝑞𝑒𝑒 1 ( 𝑠𝑠+𝑞𝑞𝑒𝑒 2 ) −𝑛𝑛𝑟𝑟 ( 𝑠𝑠+𝑞𝑞𝑒𝑒 2 )) ( 𝑎𝑎𝑎𝑎−𝑠𝑠−𝑞𝑞𝑒𝑒 2 ) 2 𝑘𝑘 � . Karena 𝑟𝑟 > 0, 𝑠𝑠 ≥ 0 , maka titik setimbang 𝐸𝐸 2 ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) ada jika 𝑎𝑎𝑏𝑏 > 𝑠𝑠 + 𝑞𝑞𝑒𝑒 2 yang bermakna titik kesetimbangan saat kedua populasi hidup berdampingan ada ketika laju pertumbuhan predator lebih besar daripada laju kematian dan laju pemanenan pada predator.
Dengan memasukkan nilai parameter yaitu 𝑟𝑟 = 2,09; 𝑠𝑠 = 0,5; 𝑘𝑘 = 1,5; 𝑛𝑛 = 0,9; 𝑎𝑎 = 0,5; 𝑏𝑏 = 2, 𝑞𝑞 = 0,003; 𝑒𝑒 1 = 1 𝑒𝑒 2 = 10 , maka diperoleh titik setimbang 𝐸𝐸 1 ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) = (1,0149; 2,5771) dan terbukti bahwa titik setimbang 𝐸𝐸 1 ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) ada karena memenuhi syarat yaitu 𝑛𝑛 ( 𝑠𝑠+𝑞𝑞𝑒𝑒 2 ) 𝑎𝑎𝑎𝑎−𝑠𝑠−𝑞𝑞𝑒𝑒 2 = 1,0149 > 0 atau 𝑎𝑎𝑏𝑏 − 𝑠𝑠 − 𝑞𝑞𝑒𝑒 2 = 0,47 > 0 .
C. Analisis Kestabilan Lokal dari Model Sistem Budidaya Ikan Bandeng di sekitar Titik Kesetimbangan
Untuk menentukan sifat kestabilan asimtotis lokal dari titik setimbang maka perlu dicari terlebih dahulu nilai eigen-nya. Karena model berbentuk nonliner maka perlu dilakukan pelinieran dengan menggunakan matriks Jacobian.
Model prey-predator pada sistem budidaya ikan bandeng dinyatakan sebagaimana pada persamaan (5).
𝑓𝑓 ( 𝑟𝑟 , 𝑠𝑠 ) = � 𝑓𝑓 1 ( 𝑟𝑟 , 𝑠𝑠 ) = 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 = 𝑟𝑟𝑟𝑟 � 1 − 𝑑𝑑 𝑘𝑘 � − 𝑎𝑎𝑑𝑑𝑥𝑥 𝑛𝑛+𝑑𝑑 − 𝑞𝑞𝑒𝑒 1 𝑟𝑟 𝑓𝑓 2 ( 𝑟𝑟 , 𝑠𝑠 ) = 𝑑𝑑𝑥𝑥 𝑑𝑑𝑑𝑑 = 𝑠𝑠 �−𝑠𝑠 + 𝑎𝑎𝑎𝑎𝑑𝑑 𝑛𝑛+𝑑𝑑 � − 𝑞𝑞𝑒𝑒 2 𝑠𝑠 (5)
Maka, matriks Jacobian dari sistem persamaan 5 adalah:
𝐽𝐽 = ⎣ ⎢ ⎢ ⎢ ⎡ 𝜕𝜕𝑓𝑓 1 𝜕𝜕𝑟𝑟 𝜕𝜕𝑓𝑓 1 𝜕𝜕𝑠𝑠 𝜕𝜕𝑓𝑓 2 𝜕𝜕𝑟𝑟 𝜕𝜕𝑓𝑓 2 𝜕𝜕𝑠𝑠 ⎦ ⎥ ⎥ ⎥ ⎤ = ⎣ ⎢ ⎢ ⎡𝑟𝑟 − 2 𝑟𝑟𝑟𝑟 𝑘𝑘 − 𝑎𝑎𝑛𝑛𝑠𝑠 ( 𝑛𝑛 + 𝑟𝑟 ) 2 − 𝑞𝑞𝑒𝑒 1 − 𝑎𝑎𝑟𝑟 𝑛𝑛 + 𝑟𝑟 𝑎𝑎𝑏𝑏𝑛𝑛𝑠𝑠
( 𝑛𝑛 + 𝑟𝑟 ) 2 −𝑠𝑠 ( 𝑛𝑛 + 𝑟𝑟 ) + 𝑎𝑎𝑏𝑏𝑟𝑟 − 𝑞𝑞𝑒𝑒 2 𝑛𝑛 + 𝑟𝑟 ⎦ ⎥ ⎥ ⎤
Selanjutnya dianalisis kestabilan di setiap titik setimbang dari model prey-predator sistem budidaya ikan bandeng.
Pada titik kesetimbangan 𝐸𝐸 0 ( 𝑟𝑟 0 , 𝑠𝑠 0 ) = (0,0) atau ketika kedua populasi punah tidak dianalisis karena tidak ada fenomena yang terjadi pada saat itu sehingga tidak ada yang bisa diamati.
Analisis Kestabilan di sekitar Titik kesetimbangan ketika kedua populasi hidup berdampingan dalam satu tambak yaitu pada saat 𝑟𝑟 ≠ 0 dan 𝑠𝑠 ≠ 0 . Apabila dimisalkan
𝛽𝛽 = 𝑛𝑛 ( 𝑠𝑠 + 𝑞𝑞𝑒𝑒 2 ) 𝑎𝑎𝑏𝑏 − 𝑠𝑠 − 𝑞𝑞𝑒𝑒 2
dengan syarat 𝑎𝑎𝑏𝑏 − 𝑠𝑠 − 𝑞𝑞𝑒𝑒 2 > 0 , maka
𝑟𝑟 ∗ = 𝛽𝛽 dan 𝑠𝑠 ∗ = ( 𝑛𝑛+𝛽𝛽 ) �𝑟𝑟− 𝑟𝑟𝑟𝑟 𝑘𝑘 −𝑞𝑞𝑒𝑒 1 � 𝑎𝑎 .
Matriks Jacobian dari titik setimbang 𝐸𝐸 1 ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) adalah:
𝐽𝐽 1 = � 𝑎𝑎 11 𝑎𝑎 12 𝑎𝑎 21 𝑎𝑎 22 �
dengan
𝑎𝑎 11 = 𝑟𝑟𝛽𝛽 ( 𝑘𝑘 − 2 𝑛𝑛 − 2 𝛽𝛽 ) + 𝛽𝛽 ( 𝑛𝑛 − 𝑞𝑞𝑒𝑒 1 𝑘𝑘 ) 𝑘𝑘 ( 𝑛𝑛 + 𝛽𝛽 ) , 𝑎𝑎 12 = − 𝑎𝑎𝛽𝛽 𝑛𝑛 + 𝛽𝛽 , 𝑎𝑎 21 = 𝑏𝑏𝑛𝑛 ( 𝑘𝑘𝑟𝑟 − 𝑟𝑟𝛽𝛽 − 𝑞𝑞𝑒𝑒 1 𝑘𝑘 ) 𝑘𝑘 ( 𝑛𝑛 + 𝛽𝛽 ) , 𝑎𝑎 22 = −𝑠𝑠𝑛𝑛 + 𝛽𝛽 ( 𝑎𝑎𝑏𝑏 − 𝑠𝑠 ) − 𝑞𝑞𝑒𝑒 2 𝑛𝑛 + 𝛽𝛽 ,
Sehingga diperoleh persamaan karakteristik yaitu
𝜆𝜆 2 + 𝐴𝐴𝜆𝜆 + 𝐵𝐵 = 0
dengan
𝐴𝐴 = − � 𝑟𝑟𝛽𝛽 ( 𝑘𝑘− 2 𝑛𝑛− 2 𝛽𝛽 ) + 𝛽𝛽 ( 𝑛𝑛−𝑞𝑞𝑒𝑒 1 𝑘𝑘 ) 𝑘𝑘 ( 𝑛𝑛 + 𝛽𝛽 ) + −𝑠𝑠𝑛𝑛 + 𝛽𝛽 ( 𝑎𝑎𝑏𝑏−𝑠𝑠 ) −𝑞𝑞𝑒𝑒 2 𝑛𝑛 + 𝛽𝛽 � 𝐵𝐵 = � 𝑟𝑟𝛽𝛽 ( 𝑘𝑘− 2 𝑛𝑛− 2 𝛽𝛽 ) + 𝛽𝛽 � 𝑛𝑛−𝑞𝑞𝑒𝑒 1 𝑘𝑘 � 𝑘𝑘 ( 𝑛𝑛 + 𝛽𝛽 ) � � −𝑠𝑠𝑛𝑛 + 𝛽𝛽 ( 𝑎𝑎𝑏𝑏−𝑠𝑠 ) −𝑞𝑞𝑒𝑒 2 𝑛𝑛 + 𝛽𝛽 � − � − 𝑎𝑎𝛽𝛽 𝑛𝑛 + 𝛽𝛽 � � 𝑏𝑏𝑛𝑛 � 𝑘𝑘𝑟𝑟−𝑟𝑟𝛽𝛽−𝑞𝑞𝑒𝑒 1 𝑘𝑘 � 𝑘𝑘 ( 𝑛𝑛 + 𝛽𝛽 ) �
Dari persamaan karakteristik diatas diperoleh dua nilai eigen sebagai berikut:
𝜆𝜆 1 = −𝐴𝐴+�𝐴𝐴 2 −4𝐵𝐵 2 dan 𝜆𝜆 2 = −𝐴𝐴−�𝐴𝐴 2 −4𝐵𝐵 2
Dari hasil nilai eigen diatas maka akan muncul dua kondisi yaitu saat nilai eigen berbentuk real dan saat nilai eigen berbentuk bilangan kompleks. Syarat agar titik setimbang stabil asimtotis adalah kedua nilai eigen bagian realnya bernilai negatif.
1. Jika 𝐴𝐴 2 − 4 𝐵𝐵 > 0 , maka nilai eigen berbentuk real dan bernilai negatif jika 𝐴𝐴 > 0 dan 𝐴𝐴 2 − 4 𝐵𝐵 ≥ 0 . 2. Jika 𝐴𝐴 2 − 4 𝐵𝐵 < 0 , maka nilai eigen berbentuk bilangan kompleks dengan bagian realnya negatif jika 𝐴𝐴 > 0 dan 𝐵𝐵 > 0 . Dengan memasukkan nilai parameter yaitu 𝑟𝑟 = 2,09; 𝑠𝑠 = 0,5; 𝑘𝑘 = 1,5; 𝑛𝑛 = 0,9; 𝑎𝑎 = 0,5; 𝑏𝑏 = 2, 𝑞𝑞 = 0,003; 𝑒𝑒 1 = 1, 𝑒𝑒 2 = 10, maka diperoleh titik setimbang 𝐸𝐸 1 ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) = (1,0149; 2,5771) . Untuk memastikan apakah sistem stabil di sekitar titik setimbang 𝐸𝐸 1 ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) kita perlu untuk mencari nilai eigen terlebih dahulu.
𝜆𝜆 1 , 2 = −𝐴𝐴 ± √𝐴𝐴 2 − 4 𝐵𝐵 2 dengan:
𝐴𝐴 = 4,9810 > 0 dan 𝐵𝐵 = 0,5408 > 0 dan 𝐴𝐴 2 − 4 𝐵𝐵 = 22,6468 sehingga:
𝜆𝜆 1 = −𝐴𝐴+�𝐴𝐴 2 −4𝐵𝐵 2 = − 0,1111 dan 𝜆𝜆 2 = −𝐴𝐴−�𝐴𝐴 2 −4𝐵𝐵 2 =
− 4,8699 . Karena kedua nilai eigen real dan negatif, maka sistem stabil asimtotis di sekitar titik kesetimbangan 𝐸𝐸 1 ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) . Dari titik kesetimbangan 𝐸𝐸 1 ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) yang dalam kondisi stabil asimtotik disertakan persamaan aljabar bioekonomi untuk menentukan keuntungan ekonomi dari panen yang telah dilakukan. Diberikan persamaan berikut ini:
𝑚𝑚 = 𝑒𝑒 1 ( 𝑝𝑝 1 𝑟𝑟 ∗ − 𝑐𝑐 ) + 𝑒𝑒 2 ( 𝑝𝑝 2 𝑠𝑠 ∗ − 𝑐𝑐 )
Selanjutnya, nilai 𝑟𝑟 ∗ dan 𝑠𝑠 ∗ disubstitusikan dalam persamaan diatas sehingga diperoleh keuntungan ekonomi.
Titik kesetimbangan 𝐸𝐸 1 ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) bergantung pada panen yang dilakukan, sehingga fungsi keuntungan juga bergantung pada kegiatan panen tersebut. Nilai pada kegiatan pemanenan 𝑒𝑒 1 dan 𝑒𝑒 2 memberikan titik kesetimbangan 𝐸𝐸 1 ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) berada dalam kondisi stabil asimtotik dan dapat memberikan keuntungan ekonomi maksimal berdasarkan hasil dari kegiatan panen pada kedua populasi.
Dengan mensubstitusikan nilai parameter yaitu 𝑟𝑟 = 2,09; 𝑠𝑠 = 0,5; 𝑘𝑘 = 1,5; 𝑛𝑛 = 0,9; 𝑎𝑎 = 0,5; 𝑏𝑏 = 2, 𝑞𝑞 = 0,003; 𝑒𝑒 1 = 1, 𝑒𝑒 2 = 10, maka diperoleh titik setimbang 𝐸𝐸 1 ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) = (1,0149; 2,5771) . Dalam sistem budidaya ikan bandeng ini, dilakukan pemanenan pada bandeng sebanyak 1 kali pada saat bandeng berumur 6 bulan, dan ikan lain yang bukan bandeng di berada dalam tambak pun juga dipanen, namun sebelum pemanenan telah dilakukan upaya penangkapan sebanyak 10 kali karena ikan lain bukan bandeng tersebut dapat menghambat pertumbuhan ikan bandeng dan juga masih bernilai ekonomis. Pada hari pertama nener dimasukkan sebanyak 15000 ekor, selama 6 bulan budidaya ikan bandeng, dari survei kondisi di lapangan diperoleh berat bandeng sekitar 200 gram sehingga 1 kg bandeng terdiri dari 5 ekor bandeng. Harga 1 kg bandeng yaitu Rp15.000 sehingga harga setiap ekor bandeng yaitu Rp3000. Ikan lain yang bukan predator juga ikut dipanen dengan kisaran harga setiap ekor ikan yang dipanen adalah Rp700. Total biaya yang dikeluarkan selama 6 bulan budidaya ikan bandeng sebesar Rp14.600.000 dengan rincian biaya Rp6.000.000 untuk biaya membeli nener atau bibit bandeng, Rp600.000 untuk biaya pupuk, Rp5.000.000 untuk biaya pakan selama 2 bulan terakhir sebelum panen, dan Rp3.000.000 untuk biaya panen. Selanjutnya, dengan memasukkan nilai parameter 𝑝𝑝 1 = 3000, 𝑝𝑝 2 = 700, dan 𝑐𝑐 = 14600000 ke dalam titik kesetimbangan 𝐸𝐸 1 ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) dengan total bandeng yang dipanen sebanyak 10149 ekor dengan berat rata-rata 200 gram sehingga diperoleh keuntungan ekonomi 𝑚𝑚 = 𝑒𝑒 1 ( 𝑝𝑝 1 𝑟𝑟 ∗ − 𝑐𝑐 ) + 𝑒𝑒 2 ( 𝑝𝑝 2 𝑠𝑠 ∗ − 𝑐𝑐 ) = 50244000 . Keuntungan ekonomi selama 6 bulan budidaya ikan bandeng yaitu Rp50.244.000.
D. Analisis Kestabilan Global menggunakan Fungsi Lyapunov
Kestabilan global pada sistem untuk titik kesetimbangan 𝐸𝐸 1 ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) diberikan oleh fungsi Lyapunov sebagai berikut [9]:
𝑉𝑉 ( 𝑟𝑟 , 𝑠𝑠 ) = 𝑟𝑟 − 𝑟𝑟 ∗ − 𝑟𝑟 ∗ ln � 𝑟𝑟 𝑟𝑟 ∗ � + 𝑠𝑠 − 𝑠𝑠 ∗ − 𝑠𝑠 ∗ ln � 𝑠𝑠 𝑠𝑠 ∗ �
dengan 𝑉𝑉 dapat dinyatakan sebagai 𝑉𝑉 ( 𝑟𝑟 , 𝑠𝑠 ) = 𝑓𝑓 𝑇𝑇 𝑓𝑓 , dengan 𝑓𝑓 merujuk pada persamaan 5 dan 𝑉𝑉 ( 𝑟𝑟 , 𝑠𝑠 ) positif dengan syarat
𝑟𝑟 − 𝑟𝑟 ∗ − 𝑟𝑟 ∗ ln � 𝑟𝑟 𝑟𝑟 ∗ � + 𝑠𝑠 − 𝑠𝑠 ∗ − 𝑠𝑠 ∗ ln � 𝑠𝑠 𝑠𝑠 ∗ � > 0,
Turunan terhadap waktu dari 𝑉𝑉 ( 𝑟𝑟 , 𝑠𝑠 ) sistem diatas yaitu
𝑑𝑑𝑉𝑉 𝑑𝑑𝑡𝑡 = 𝑟𝑟 − 𝑟𝑟 ∗ 𝑟𝑟 𝑑𝑑𝑟𝑟 𝑑𝑑𝑡𝑡 + 𝑠𝑠 − 𝑠𝑠 ∗ 𝑠𝑠 𝑑𝑑𝑠𝑠 𝑑𝑑𝑡𝑡 ,
dimana 𝑉𝑉̇ dapat dinyatakan sebagai:
𝑉𝑉̇ = 𝑓𝑓 𝑇𝑇 ̇ 𝑓𝑓 + 𝑓𝑓 𝑇𝑇 𝑓𝑓̇ 𝑉𝑉̇ = ( 𝐽𝐽𝑓𝑓 ) 𝑇𝑇 𝑓𝑓 + 𝑓𝑓 𝑇𝑇 𝐽𝐽𝑓𝑓 𝑉𝑉̇ = ( 𝐽𝐽𝑓𝑓 ) 𝑇𝑇 𝑓𝑓 + 𝑓𝑓 𝑇𝑇 𝐽𝐽𝑓𝑓
𝑉𝑉̇ = 𝑓𝑓 𝑇𝑇 𝐽𝐽 𝑇𝑇 𝑓𝑓 + 𝑓𝑓 𝑇𝑇 𝐽𝐽𝑓𝑓 𝑉𝑉̇ = 𝑓𝑓 𝑇𝑇 ( 𝐽𝐽 𝑇𝑇 + 𝐽𝐽 ) 𝑓𝑓 sehingga, 𝑉𝑉̇ dapat disederhanakan menjadi
𝑉𝑉̇ = 𝑓𝑓𝑓𝑓𝑓𝑓 𝑇𝑇
dengan 𝑓𝑓 adalah 𝑓𝑓 = 𝐽𝐽 𝑇𝑇 + 𝐽𝐽 , sehingga:
𝑓𝑓 = − 2 � 2𝑟𝑟𝑑𝑑 𝑘𝑘 + 𝑎𝑎𝑛𝑛𝑥𝑥 ( 𝑛𝑛+𝑑𝑑 ) 2 + 𝑒𝑒 1 − 𝑟𝑟 1 2 � 𝑎𝑎𝑑𝑑 ( 𝑛𝑛+𝑑𝑑 ) −𝑎𝑎𝑎𝑎𝑛𝑛𝑥𝑥
( 𝑛𝑛+𝑑𝑑 ) 2 � 1 2 � 𝑎𝑎𝑑𝑑 ( 𝑛𝑛+𝑑𝑑 ) −𝑎𝑎𝑎𝑎𝑛𝑛𝑥𝑥
( 𝑛𝑛+𝑑𝑑 ) 2 � 𝑠𝑠 ( 𝑛𝑛+𝑑𝑑 ) +𝑞𝑞𝑒𝑒 2 −𝑎𝑎𝑎𝑎𝑑𝑑 ( 𝑛𝑛+𝑑𝑑 ) �
maka, nilai 𝑉𝑉̇ adalah
𝑉𝑉̇ = − 2 �𝑟𝑟 − 𝑟𝑟 ∗ 𝑠𝑠 − 𝑠𝑠 ∗ � � 2𝑟𝑟𝑑𝑑 𝑘𝑘 + 𝑎𝑎𝑛𝑛𝑥𝑥 ( 𝑛𝑛+𝑑𝑑 ) 2 + 𝑒𝑒 1 − 𝑟𝑟 1 2 � 𝑎𝑎𝑑𝑑 ( 𝑛𝑛+𝑑𝑑 ) −𝑎𝑎𝑎𝑎𝑛𝑛𝑥𝑥 ( 𝑛𝑛+𝑑𝑑 ) 2 � 1 2 � 𝑎𝑎𝑑𝑑 ( 𝑛𝑛+𝑑𝑑 ) −𝑎𝑎𝑎𝑎𝑛𝑛𝑥𝑥 ( 𝑛𝑛+𝑑𝑑 ) 2 � 𝑠𝑠 ( 𝑛𝑛+𝑑𝑑 ) +𝑞𝑞𝑒𝑒 2 −𝑎𝑎𝑎𝑎𝑑𝑑 ( 𝑛𝑛+𝑑𝑑 ) � .
𝑉𝑉̇ definit negatif, jika matriks 𝑓𝑓 definit positif [5]. Matriks 𝑓𝑓 definit positif jika:
1. 𝑎𝑎 11 > 0 ,
2𝑟𝑟𝑑𝑑 𝑘𝑘 + 𝑎𝑎𝑛𝑛𝑥𝑥 ( 𝑛𝑛+𝑑𝑑 ) 2 + 𝑒𝑒 1 − 𝑟𝑟 > 0 .
2. 𝑎𝑎 11 𝑎𝑎 22 − 𝑎𝑎 12 𝑎𝑎 21 > 0,
� 2𝑟𝑟𝑑𝑑 𝑘𝑘 + 𝑎𝑎𝑛𝑛𝑥𝑥 ( 𝑛𝑛+𝑑𝑑 ) 2 + 𝑒𝑒 1 − 𝑟𝑟� � 𝑠𝑠 ( 𝑛𝑛+𝑑𝑑 ) +𝑞𝑞𝑒𝑒 2 −𝑎𝑎𝑎𝑎𝑑𝑑 ( 𝑛𝑛+𝑑𝑑 ) � − � 1 2 � 𝑎𝑎𝑑𝑑 ( 𝑛𝑛+𝑑𝑑 ) −𝑎𝑎𝑎𝑎𝑛𝑛𝑥𝑥
( 𝑛𝑛+𝑑𝑑 ) 2 �� 2 > 0 .
Dengan memasukkan nilai parameter yaitu 𝑟𝑟 = 2,09; 𝑠𝑠 = 0,5; 𝑘𝑘 = 1,5; 𝑛𝑛 = 0,9; 𝑎𝑎 = 0,5; 𝑏𝑏 = 2, 𝑞𝑞 = 0,003; 𝑒𝑒 1 = 1, 𝑒𝑒 2 = 10, maka:
𝑎𝑎 11 = 2𝑟𝑟𝑑𝑑 𝑘𝑘 + 𝑎𝑎𝑛𝑛𝑥𝑥 ( 𝑛𝑛+𝑑𝑑 ) 2 + 𝑒𝑒 1 − 𝑟𝑟 = 2,0544 > 0 dan � 2𝑟𝑟𝑑𝑑 𝑘𝑘 + 𝑎𝑎𝑛𝑛𝑥𝑥 ( 𝑛𝑛+𝑑𝑑 ) 2 + 𝑒𝑒 1 − 𝑟𝑟� � 𝑠𝑠 ( 𝑛𝑛+𝑑𝑑 ) +𝑞𝑞𝑒𝑒 2 −𝑎𝑎𝑎𝑎𝑑𝑑 ( 𝑛𝑛+𝑑𝑑 ) � − � 1 2 � 𝑎𝑎𝑑𝑑 ( 𝑛𝑛+𝑑𝑑 ) −𝑎𝑎𝑎𝑎𝑛𝑛𝑥𝑥 ( 𝑛𝑛+𝑑𝑑 ) 2 �� 2 =
0.0632 > 0 .
Dari uraian diatas, diketahui bahwa 𝑉𝑉 ( 𝑟𝑟 , 𝑠𝑠 ) > 0 dan 𝑉𝑉̇ < 0 sehingga memenuhi kondisi pada fungsi lyapunov dan dapat dikatakan bahwa sistem stabil asimtotik global, karena memenuhi kondisi 1 dan 2.
E. Eksistensi Bifurkasi pada Model Sistem Budidaya Ikan Bandeng
Untuk persamaan diferensial bioekonomi pada model prey-predator dimana rate pemanenan pada prey ( 𝑒𝑒 1 ) diasumsikan sebagai parameter bifurkasi, akan diselidiki adanya bifurkasi dalam model.
Diberikan model prey-predator dengan pemanenan prey- predator proporsional dan respons fungsional Holling tipe II dengan penambahan persamaan diferensial bioekonomi sebagai berikut:
𝑑𝑑𝑟𝑟 𝑑𝑑𝑡𝑡 = 𝑟𝑟𝑟𝑟 � 1 − 𝑟𝑟 𝑘𝑘� − 𝑎𝑎𝑟𝑟𝑠𝑠 𝑛𝑛 + 𝑟𝑟 − 𝑞𝑞𝑒𝑒 1 𝑟𝑟 ,
𝑑𝑑𝑠𝑠 𝑑𝑑𝑡𝑡 = 𝑠𝑠 �−𝑠𝑠 + 𝑎𝑎𝑏𝑏𝑟𝑟 𝑛𝑛 + 𝑟𝑟� − 𝑞𝑞𝑒𝑒 2 𝑠𝑠 ,
𝑚𝑚 = 𝑒𝑒 1 ( 𝑝𝑝 1 𝑟𝑟 − 𝑐𝑐 ) + 𝑒𝑒 2 ( 𝑝𝑝 2 𝑠𝑠 − 𝑐𝑐 ).
Karena 𝑚𝑚 merupakan variabel yang bergantung pada 𝑟𝑟̇ ( 𝑡𝑡 ) dan 𝑠𝑠̇ ( 𝑡𝑡 ) yang mengalami perubahan setiap waktu maka bisa ditulis menjadi
𝑑𝑑𝑚𝑚 𝑑𝑑𝑡𝑡 = 𝑒𝑒 1 �𝑝𝑝 1 � 𝑑𝑑𝑟𝑟 𝑑𝑑𝑡𝑡� − 𝑐𝑐� + 𝑒𝑒 2 �𝑝𝑝 2 � 𝑑𝑑𝑠𝑠 𝑑𝑑𝑡𝑡� − 𝑐𝑐� sehingga model menjadi:
𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 = 𝑟𝑟𝑟𝑟 � 1 − 𝑑𝑑 𝑘𝑘 � − 𝑎𝑎𝑑𝑑𝑥𝑥 𝑛𝑛+𝑑𝑑 − 𝑞𝑞𝑒𝑒 1 𝑟𝑟 , 𝑑𝑑𝑥𝑥 𝑑𝑑𝑑𝑑 = 𝑠𝑠 �−𝑠𝑠 + 𝑎𝑎𝑎𝑎𝑑𝑑 𝑛𝑛+𝑑𝑑 � − 𝑞𝑞𝑒𝑒 2 𝑠𝑠 , 𝑑𝑑𝑚𝑚 𝑑𝑑𝑑𝑑 = 𝑒𝑒 1 𝑝𝑝 1 �𝑟𝑟𝑟𝑟 − 𝑟𝑟𝑑𝑑 2 𝑘𝑘 − 𝑎𝑎𝑑𝑑𝑥𝑥 𝑛𝑛+𝑑𝑑 − 𝑞𝑞𝑒𝑒 1 𝑟𝑟� + 𝑒𝑒 2 𝑝𝑝 2 � 𝑎𝑎𝑎𝑎𝑑𝑑𝑥𝑥 𝑛𝑛+𝑑𝑑 − 𝑠𝑠𝑠𝑠 − 𝑞𝑞𝑒𝑒 2 𝑠𝑠� − ( 𝑒𝑒 1 𝑐𝑐 + 𝑒𝑒 2 𝑐𝑐 ) (6)
maka, matriks Jacobian dari sistem diatas adalah:
𝐽𝐽 = ⎣ ⎢ ⎢ ⎢ ⎡ 𝑟𝑟 − 2𝑟𝑟𝑑𝑑 𝑘𝑘 − 𝑎𝑎𝑛𝑛𝑥𝑥 ( 𝑛𝑛+𝑑𝑑 ) 2 − 𝑞𝑞𝑒𝑒 1 − 𝑎𝑎𝑑𝑑 𝑛𝑛+𝑑𝑑 0 𝑎𝑎𝑎𝑎𝑛𝑛𝑥𝑥 ( 𝑛𝑛+𝑑𝑑 ) 2 −𝑠𝑠 ( 𝑛𝑛+𝑑𝑑 ) +𝑎𝑎𝑎𝑎𝑑𝑑−𝑞𝑞𝑒𝑒 2 𝑛𝑛+𝑑𝑑 0 𝑒𝑒 1 𝑝𝑝 1 𝑟𝑟 − 2𝑟𝑟𝑑𝑑𝑒𝑒 1 𝑝𝑝 1 𝑘𝑘 − 𝑒𝑒 1 𝑝𝑝 1 𝑎𝑎𝑛𝑛𝑥𝑥
( 𝑛𝑛+𝑑𝑑 ) 2 − 𝑒𝑒 1 𝑝𝑝 1 𝑞𝑞 + 𝑒𝑒 2 𝑝𝑝 2 𝑎𝑎𝑎𝑎𝑛𝑛𝑥𝑥 ( 𝑛𝑛+𝑑𝑑 ) 2 − 𝑒𝑒 1 𝑝𝑝 1 𝑎𝑎𝑑𝑑 𝑛𝑛+𝑑𝑑 − 𝑒𝑒 2 𝑝𝑝 2 𝑎𝑎𝑎𝑎𝑑𝑑 𝑛𝑛+𝑑𝑑 − 𝑒𝑒 2 𝑝𝑝 2 𝑠𝑠 − 𝑒𝑒 2 2 𝑝𝑝 2 𝑞𝑞 0 ⎦ ⎥ ⎥ ⎥ ⎤
sehingga  dari matriks Jacobian diatas  diperoleh nilai karakteristik yaitu:
𝑎𝑎 0 𝜆𝜆 3 + 𝑎𝑎 1 𝜆𝜆 2 + 𝑎𝑎 2 𝜆𝜆 + 𝑎𝑎 3 = 0 dengan:
𝑎𝑎 0 = 1 ,
𝑎𝑎 1 = − ��𝑟𝑟 − 2𝑟𝑟𝑑𝑑 𝑘𝑘 − 𝑎𝑎𝑛𝑛𝑥𝑥 ( 𝑛𝑛+𝑑𝑑 ) 2 − 𝑞𝑞𝑒𝑒 1 � + � −𝑠𝑠 ( 𝑛𝑛+𝑑𝑑 ) +𝑎𝑎𝑎𝑎𝑑𝑑−𝑞𝑞𝑒𝑒 2 𝑛𝑛+𝑑𝑑 �� ,
𝑎𝑎 2 = �𝑟𝑟 − 2𝑟𝑟𝑑𝑑 𝑘𝑘 − 𝑎𝑎𝑛𝑛𝑥𝑥 ( 𝑛𝑛+𝑑𝑑 ) 2 − 𝑞𝑞𝑒𝑒 1 � � −𝑠𝑠 ( 𝑛𝑛+𝑑𝑑 ) +𝑎𝑎𝑎𝑎𝑑𝑑−𝑞𝑞𝑒𝑒 2 𝑛𝑛+𝑑𝑑 � ,
𝑎𝑎 3 = 0 .
Kestabilan sistem diatas dapat dianalis dengan kriteria Routh-Hurwitz. Kriteria kestabilan Routh-Hurwitz digunakan untuk mengecek langsung kestabilan melalui koefisien tanpa menghitung akar - akar dari polinomial yang ada. Dari persamaan karakteristik
𝑎𝑎 0 𝜆𝜆 3 + 𝑎𝑎 1 𝜆𝜆 2 + 𝑎𝑎 2 𝜆𝜆 + 𝑎𝑎 3 = 0
Kestabilan sistem dianalisis dengan kriteria Routh- Hurwitz dengan menyusunnya sebagai berikut:
𝑠𝑠 3 𝑎𝑎 0 𝑎𝑎 2
𝑠𝑠 2 𝑎𝑎 1 𝑎𝑎 3 𝑠𝑠 1 𝑎𝑎 1 𝑎𝑎 2 − 𝑎𝑎 0 𝑎𝑎 3 𝑎𝑎 1 𝑠𝑠 0 𝑎𝑎 3
Syarat agar persamaan diferensial stabil adalah semua variabel pada kolom pertama bertanda positif. Maka, sistem stabil jika:
1. 𝑎𝑎 1 , 𝑎𝑎 2 > 0 , 2. 𝑎𝑎 1 𝑎𝑎 2 −𝑎𝑎 0 𝑎𝑎 3 𝑎𝑎 1 > 0 atau 𝑎𝑎 1 𝑎𝑎 2 > 0 .
Dari uraian sebelumnya, dijelaskan bahwa persamaan 6 akan bersifat stabil jika 𝑎𝑎 1 > 0 , 𝑎𝑎 2 > 0 dan 𝑎𝑎 1 𝑎𝑎 2 > 0 . Karena pada titik ini berlaku polinomial yaitu 𝑎𝑎 0 𝑠𝑠 3 + 𝑎𝑎 1 𝑠𝑠 2 + 𝑎𝑎 2 𝑠𝑠 + 𝑎𝑎 3 = 0 , maka berdasarkan teorema bifurkasi Hopf pada titik tetap tersebut terjadi jika 𝑎𝑎 1 > 0, 𝑎𝑎 2 > 0 dan 𝑎𝑎 1 𝑎𝑎 2 = 0 sehingga 𝑑𝑑 𝑑𝑑𝑒𝑒 2 ( 𝑎𝑎 1 𝑎𝑎 2 ) ≠ 0 .

## F. Simulasi Mode
Simulasi dilakukan dengan memasukkan data dan nilai parameter yang diperoleh kedalam persamaan model matematika Prey-Predator bioekonomi dengan pemanenan Prey-Predator proporsional dan fungsi Holling tipe II. Simulasi ini dilakukan dengan metode runge-kutta orde 4 dengan bantuan software MATLAB R2021a. Nilai awal dan parameter dapat dilihat pada Tabel 1.
Pada simulasi ini diberikan nilai awal yaitu 𝑟𝑟 = 1,5 , dan 𝑠𝑠 = 0,0025 dengan waktu akhir adalah 180 hari. Berdasarkan nilai parameter dan nilai awal, diperoleh titik kesetimbangan 𝐸𝐸 1 ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) = (1,0149; 2,5771) . Dinamika populasi Prey-Predator digambarkan dalam kurva bidang solusi untuk memudahkan kita melihat kestabilan populasi Prey-Predator ( 𝑟𝑟 , 𝑠𝑠 ) dalam jangka waktu 𝑡𝑡 . Berikut ini ditunjukkan grafik bidang solusi populasi prey dan predator dalam waktu 180 hari.
Gambar 4 merupakan grafik bidang solusi populasi prey dan predator dengan diperoleh titik kesetimbangan 𝐸𝐸 1 = (1,0149, 2,5771) . Berdasarkan gambar tersebut, ditunjukkan bahwa model stabil asimtotik pada titik kesetimbangan 𝐸𝐸 1 . Kestabilan dinamika populasi prey dan predator adalah simpul stabil. Gambar 4 menunjukkan pada saat 𝑡𝑡 = 0 sampai 𝑡𝑡 = 30 , 15000 ekor bibit bandeng atau yang biasa disebut nener dimasukkan dalam tambak dan masih belum ada predator yang masuk karena tambak sudah dibersihkan dan diberi pupuk sebelum dimasuki nener. Populasi bandeng ini cenderung stabil dan jika kita perbesar gambarnya populasi bandeng sedikit menurun karena nener yang tidak mampu beradaptasi dengan lingkungan baru mengalami kematian serta adanya keterbatasan daya dukung lingkungannya karena tidak ada suplai makanan dari petani saat bandeng berumur 0-30 hari. Pada hari ke-31, dari 15000 ekor bandeng yang dimasukkan, masih bertahan 14970 ekor bandeng yang nantinya akan bertumbuh semakin besar dalam tambak. Benih bandeng mencari makanan sendiri dalam tambak berupa gambut-gambut yang berada di dasar tambak. Setelah 𝑡𝑡 = 31 , ikan lain yang bukan bandeng mulai masuk ke tambak karena adanya air pasang dari laut. Sampai di hari ke- 70, terdapat 14821 ekor bandeng yang masih bertahan, namun akibat adanya interaksi dengan predator, persediaan makanan yang semakin berkurang, serta tempat yang semakin sempit karena adanya ikan lain yang bukan bandeng,
populasi bandeng ini mengalami penurunan yang drastis mulai hari ke-80 sampai hari ke-130 hingga tersisa 10220 ekor yang sudah mulai membesar. Pada hari ke 131-180 populasi prey cenderung stabil karena tubuh ikan mulai membesar hingga di waktu akhir budidaya atau saat pemanenan didapatkan 10149 ekor dengan berat rata-rata 200 gram. Pada hari ke-151 sampai hari ke-180 ini ikan bandeng diberikan pakan  tambahan setiap harinya sebanyak 2 kali. Dalam sistem budidaya ikan bandeng ini, dilakukan pemanenan pada bandeng sebanyak 1 kali pada saat bandeng berumur 6 bulan, dan ikan lain yang bukan bandeng di berada dalam tambak pun juga dipanen, namun sebelum pemanenan telah dilakukan upaya penangkapan sebanyak 10 kali karena ikan lain bukan bandeng tersebut dapat menghambat pertumbuhan ikan bandeng dan juga masih bernilai ekonomis. Selama 6 bulan budidaya ikan bandeng, dari survei kondisi di lapangan diperoleh berat bandeng sekitar 200 gram sehingga 1 kg bandeng terdiri dari 5 ekor bandeng. Harga 1 kg bandeng yaitu Rp15.000 sehingga harga setiap ekor bandeng yaitu Rp3000. Ikan lain yang bukan predator juga ikut dipanen dengan kisaran harga setiap ekor ikan yang dipanen adalah Rp700. Total biaya yang dikeluarkan selama 6 bulan budidaya ikan bandeng sebesar Rp14.600.000 dengan rincian biaya Rp6.000.000 untuk biaya membeli nener atau bibit bandeng, Rp600.000 untuk biaya pupuk, Rp5.000.000 untuk biaya pakan selama 2 bulan terakhir sebelum panen, dan Rp3.000.000 untuk biaya panen. Selanjutnya, dengan memasukkan nilai parameter 𝑝𝑝 1 = 3000, 𝑝𝑝 2 = 700, dan 𝑐𝑐 = 14600000 ke dalam titik kesetimbangan 𝐸𝐸 1 ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) dengan total bandeng yang dipanen sebanyak 10149 ekor dengan berat rata-rata 200 gram sehingga diperoleh keuntungan ekonomi 𝑚𝑚 = 𝑒𝑒 1 ( 𝑝𝑝 1 𝑟𝑟 ∗ − 𝑐𝑐 ) + 𝑒𝑒 2 ( 𝑝𝑝 2 𝑠𝑠 ∗ − 𝑐𝑐 ) = 50244000 . Keuntungan ekonomi selama 6 bulan budidaya ikan bandeng yaitu Rp50.244.000.
Mensubstitusikan nilai parameter yaitu 𝑟𝑟 = 2,09; 𝑠𝑠 = 0,5; 𝑘𝑘 = 1,5; 𝑛𝑛 = 0,9; 𝑎𝑎 = 0,5; 𝑏𝑏 = 2, 𝑞𝑞 = 0,003; 𝑒𝑒 2 = 10 , dengan merubah nilai 𝑒𝑒 1 sebanyak 3 kali yaitu 𝑒𝑒 1 = 1, 𝑒𝑒 2 = 3 , dan 𝑒𝑒 3 = 5 akan ditunjukkan eksistensi bifurkasi pada model.
Ketika 𝑒𝑒 1 = 1 diperoleh nilai 𝑎𝑎 1 = 1,0431 > 0 , 𝑎𝑎 2 = 0,1525 > 0 , dan 𝑎𝑎 1 𝑎𝑎 2 = 0,1590 > 0 . Sehingga dapat ditunjukkan bahwa sistem stabil karena memenuhi syarat i dan ii. Namun karena 𝑎𝑎 1 𝑎𝑎 2 ≠ 0 maka sistem ini tidak terjadi bifurkasi.
Ketika 𝑒𝑒 1 = 3 diperoleh nilai 𝑎𝑎 1 = 1,0463 > 0 , 𝑎𝑎 2 = 0,1509 > 0 , dan 𝑎𝑎 1 𝑎𝑎 2 = 0,1579 > 0 . Sehingga dapat ditunjukkan bahwa  sistem stabil karena memenuhi syarat i dan ii. Namun karena 𝑎𝑎 1 𝑎𝑎 2 ≠ 0 maka sistem ini tidak terjadi bifurkasi.
Ketika 𝑒𝑒 1 = 5 diperoleh nilai 𝑎𝑎 1 = 1,0495 > 0 , 𝑎𝑎 2 = 0,1494 > 0 , dan 𝑎𝑎 1 𝑎𝑎 2 = 0,1568 > 0 . Sehingga dapat ditunjukkan bahwa sistem stabil karena memenuhi syarat i dan ii. Namun karena 𝑎𝑎 1 𝑎𝑎 2 ≠ 0 maka sistem ini tidak terjadi bifurkasi.
Sehingga pada parameter rate pemanenan pada prey atau 𝑒𝑒 1 yang diubah sebanyak 3 kali tidak mempengaruhi stabilitas, sehingga eksistensi bifurkasi hopf tidak terjadi.
Dalam sistem  budidaya ikan bandeng sendiri, sistem ini dikatakan stabil jika populasi prey dan predator selalu ada atau tidak pernah habis selama 6 bulan waktu budidaya tersebut. Pada sistem budidaya ikan bandeng di Kecamatan
Manyar, Kabupaten Gresik diasumsikan bahwa  preynya adalah ikan bandeng, sedangkan predatornya yaitu ikan jenis lain yang bukan bandeng seperti kepiting, belut, udang, ikan laut, ikan patin yang mana predator ini masuk ke tambak saat adanya pasang laut. Keberadaan ikan jenis lain ini dapat menghambat pertumbuhan ikan bandeng karena mengurangi persediaan makanan yang ada di dasar tambak, juga mempersempit tempat hidup dari ikan bandeng. Dari simulasi model juga diperoleh keuntungan ekonomi yang optimal ketika prey dan predator sama-sama dipanen.

## IV. KESIMPULAN
Berdasarkan analisis model prey-predator bioekonomi dengan pemanenan proporsional dan fungsi Holling tipe II yang diaplikasikan pada sistem budidaya ikan bandeng di kecamatan Manyar, Kabupaten Gresik dapat diperoleh kesimpulan sebagai berikut: (1)  Dengan memasukkan nilai parameter yaitu 𝑟𝑟 = 2,09; 𝑠𝑠 = 0,5; 𝑘𝑘 = 1,5; 𝑛𝑛 = 0,9; 𝑎𝑎 = 0,5; 𝑏𝑏 = 2, 𝑞𝑞 = 0,003; 𝑒𝑒 1 = 1, 𝑒𝑒 2 = 10 ,  dapat diketahui bahwa sistem stabil asimtotis lokal maupun global pada titik kesetimbangan 𝐸𝐸 1 ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) = (1,0149; 2,5771) .  Dalam sistem budidaya ikan bandeng yaitu  sistem ini dikatakan stabil jika populasi prey dan predator selalu ada atau tidak pernah habis selama 6 bulan waktu budidaya tersebut. Sistem ini juga stabil apabila laju pertumbuhan populasi prey ini lebih besar daripada koefisien predasi oleh predator. Dari titik kesetimbangan 𝐸𝐸 1 ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) yang dalam kondisi stabil asimtotik dapat diperoleh keuntungan ekonomi yang optimal dari panen yang telah dilakukan. Titik kesetimbangan 𝐸𝐸 1 ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) bergantung pada panen yang dilakukan, sehingga fungsi keuntungan juga bergantung pada kegiatan panen tersebut. Nilai pada kegiatan pemanenan 𝑒𝑒 1 dan 𝑒𝑒 2 memberikan titik kesetimbangan 𝐸𝐸 1 ( 𝑟𝑟 ∗ , 𝑠𝑠 ∗ ) berada dalam kondisi stabil asimtotik dan dapat memberikan keuntungan ekonomi optimal berdasarkan hasil dari kegiatan panen pada kedua populasi. Dengan mensubstitusikan nilai parameter diperoleh keuntungan ekonomi selama 6 bulan budidaya ikan bandeng yaitu Rp50.244.000. (2) Dengan mengubah nilai parameter rate pemanenan pada prey ( 𝑒𝑒 1 ) sebanyak 3 kali yaitu 𝑒𝑒 1 = 1 , 𝑒𝑒 1 = 3 , dan 𝑒𝑒 1 = 5 , stabilitas dari sistem tidak berubah atau tetap stabil. Hal ini menunjukkan bahwa dalam sistem ini tidak terjadi bifurkasi. (3) Dari analisis simulasi
yang telah dilakukan, interpretasi ke dalam sistem budidaya ikan bandeng yaitu sistem ini dikatakan stabil jika populasi prey dan predator selalu ada atau tidak pernah habis selama 6 bulan waktu budidaya tersebut. Pada sistem budidaya ikan bandeng di Kecamatan Manyar, Kabupaten Gresik diasumsikan bahwa preynya adalah ikan bandeng, sedangkan predatornya yaitu ikan jenis lain yang bukan bandeng seperti kepiting, belut, udang, ikan laut, ikan patin yang mana predator ini masuk ke tambak saat adanya pasang laut. Keberadaan ikan jenis lain ini dapat menghambat pertumbuhan ikan bandeng karena mengurangi persediaan makanan yang ada di dasar tambak, juga mempersempit tempat hidup dari ikan bandeng. Dari simulasi model juga diperoleh keuntungan ekonomi yang optimal ketika prey dan predator sama-sama dipanen.

## DAFTAR PUSTAKA
[1] L. R. Sari, “Model matematika pemanenan ikan dengan kebijakan panen selektif,” Limits: Journal of Mathematics and Its Applications , vol. 17, no. 1, p. 83, 2020, doi: 10.12962/limits.v17i1.6753.
[2] N. Bellomo, M. Delitala, and P. D. Torino, Lecture Notes on Mathematical Modelling From Applied Sciences to Complex Systems , vol. 8. Roma: Società Italiana di Matematica Applicata e Industriale, 2010.
[3] M. Li, B. Chen, and H. Ye, “A bioeconomic differential algebraic predator–prey model with nonlinear prey harvesting,” Appl Math Model , vol. 42, pp. 17–28, Feb. 2017, doi: https://doi.org/10.1016/j.apm.2016.09.029.
[4] Y. Pei, M. Chen, X. Liang, and C. Li, “Model-based on fishery management systems with selective harvest policies,” Math Comput Simul , vol. 156, pp. 377–395, Feb. 2019, doi: https://doi.org/10.1016/j.matcom.2018.08.009.
[5] R. Sundari and E. Apriliani, “Konstruksi fungsi lyapunov  untuk menentukan kestabilan,” Jurnal Sains dan Seni ITS , vol. 6, no. 1, Mar. 2017, doi: 10.12962/j23373520.v6i1.22862.
[6] Y. A. Kuznetsov, Elements of Applied Bifurcation Theory, Second Edition , 2nd ed., vol. 112. New York: Springer, 1998.
[7] M. S. W. Sunaryo, Z. Salleh, and M. Mamat, “Mathematical model of three species food chain with holling type-III functional response,” International Journal of Pure and Applied Mathematics , vol. 89, no. 5, pp. 647–657, Apr. 2013, doi: 10.12732/ijpam.v89i5.1.
[8] M. R. Alydrus, “Analisis Dinamik Model Predator Prey Leslie-Gower dengan Fungsi Respon Sokol-Howel,” Matematika, Institut Teknologi Sepuluh Nopember, Surabaya, 2015.
[9] Sutimin, S. Khabibah, and D. Anis Munawwaroh, “Analysis of a fishery model with two competing prey species in the presence of a predator species for optimal harvesting,” E3S Web Conf. , vol. 31, 2018.