PEMANENAN OPTIMAL BERDASARKAN SISTEM KOMPETISI PADA TIGA SPESIES IKAN
Oleh: WAHYU HARTONO G05400048
DEPARTEMEN MATEMATIKA FAKULTAS MATEMATIKA DAN ILMU PENGETAHUAN ALAM INSTITUT PERTANIAN BOGOR BOGOR 2007
ABSTRACT WAHYU HARTONO. An Optimal Harvesting Based on Competitive System in Three Species Fishery. Supervised by FARIDA HANUM and SISWANDI The problem of fish harvest policy which give maximum profit and continuity (the fish population exctinction is not happen) is an important things for the fishery industries. Many scientists try to solve this problem with many ways and one of them is to make a mathematical model and then solve the model with appropriate method. In this paper, a mathematical model is proposed and analysed to study the dynamics of oneprey two-predators system with harvesting toward the predators by the fisherman. Criteria for global stability of the positive equilibria are obtained by the Lyapunov method. The problem of optimal harvest policy is solved by using Pontryagin’s maximal principle. Finally, computer simulations are performed to investigate the dynamics of the system and optimal behavior of harvesting by choosing the value of the system’s parameters. The global behavior of the system has been proved by suitable Lyapunov function that has been constructed, and assuming some parametric restrictions. It has been found that the system is globally asymptotically stable in certain parametric regions. With the given value of parameters, it has been obtained the optimal harvesting for the predators.
RINGKASAN WAHYU HARTONO. Pemanenan Optimal Berdasarkan Sistem Kompetisi pada Tiga Spesies Ikan. Dibimbing oleh FARIDA HANUM dan SISWANDI Permasalahan kebijakan pemanenan ikan yang memberikan keuntungan maksimum dan berkelanjutan (tidak terjadi kepunahan dari populasi ikan yang dipanen) adalah hal yang sangat penting bagi industri perikanan. Para ilmuwan berusaha menyelesaikan permasalahan tersebut dengan berbagai cara, salah satunya adalah dengan memodelkan permasalahan tersebut secara matematis kemudian menyelesaikannya dengan menggunakan metode yang sesuai untuk model tersebut. Dalam tulisan ini, model matematika diajukan dan dianalisis untuk mempelajari kedinamikan sistem satu-mangsa dua-pemangsa dari tiga spesies ikan dengan pemanenan terhadap kedua pemangsa oleh para nelayan. Kriteria kestabilan global untuk titik ekuilibrium positif diperoleh dengan metode Lyapunov. Masalah kebijakan pemanenan optimal diselesaikan dengan menggunakan prinsip maksimum Pontryagin. Akhirnya, simulasi komputer ditampilkan untuk menyelidiki kedinamikan dari sistem dan perilaku optimal pemanenan dengan memberikan nilai parameter-parameter sistemnya. Perilaku global dari sistem telah dibuktikan secara detail dengan menggunakan fungsi Lyapunov dan beberapa batasan parameter. Sistem interaksi tiga spesies ikan seperti yang dimodelkan pada tulisan ini bersifat stabil asimtotik global pada beberapa wilayah parameter. Dengan memberikan nilai-nilai parameternya, diperoleh level hasil pemanenan optimal kedua predator (ikan-ikan yang dipanen).
PEMANENAN OPTIMAL BERDASARKAN SISTEM KOMPETISI PADA TIGA SPESIES IKAN
Skripsi sebagai salah satu syarat untuk memperoleh gelar Sarjana Sains pada Fakultas Matematika dan Ilmu Pengetahuan Alam Institut Pertanian Bogor
WAHYU HARTONO G05400048
DEPARTEMEN MATEMATIKA FAKULTAS MATEMATIKA DAN ILMU PENGETAHUAN ALAM INSTITUT PERTANIAN BOGOR BOGOR 2007
Judul
:
Nama NRP
: :
Pemanenan Optimal Berdasarkan Sistem Kompetisi pada Tiga Spesies Ikan Wahyu Hartono G05400048
Menyetujui: Pembimbing I,
Pembimbing II,
Dra. Farida Hanum, M.Si. NIP 131 956 709
Drs. Siswandi, M.Si. NIP 131 957 320
Mengetahui: Dekan Fakultas Matematika dan Ilmu Pengetahuan Alam Institut Pertanian Bogor
Prof. Dr. Ir. Yonny Koesmaryono, MS NIP 131 473 999
Tanggal Lulus:
PRAKATA Alhamdulillah, puji dan syukur penulis panjatkan kepada Allah SWT yang senantiasa melimpahkan rahmat, karunia, serta nikmat-Nya. Sehingga penulis dapat menyelesaikan skripsi yang merupakan salah satu syarat kelulusan program sarjana pada Departemen Matematika, Fakultas Matematika dan Ilmu Pengetahuan Alam, Institut Pertanian Bogor. Dalam penyelesaian tugas akhir ini, penulis banyak dibantu oleh berbagai pihak. Terima kasih penulis haturkan untuk: 1. Ayah dan Umi atas semua kasih sayang, semangat, bekal, doa, ketabahan, dan kesabaran yang luar biasa yang diberikan selama penulis melaksanakan studi di Institut Pertanian Bogor. 2. Ibu Farida Hanum selaku pembimbing I yang telah dengan sabar membantu, membimbing, dan berbagi ilmu pengetahuannya kepada penulis. 3. Bapak Siswandi selaku pembimbing II, atas pengertian, pengalaman dan ilmu pengetahuan yang telah diberikan kepada penulis. 4. Bapak Ali Kusnanto atas kesediaannya menjadi penguji luar dalam sidang tugas akhir penulis. 5. Bapak H. Harsono beserta Ibu Hj. Ety Riani yang telah memberikan banyak sekali bantuan (tempat kost, beasiswa, dll) ketika penulis sedang menghadapi saat - saat kritis dalam hidup. 6. Rian, Yudi, dan Anton atas kesediaannya menjadi pembahas dalam seminar penulis. 7. Mbak Suci, Mbak Murni, Mas Eko dan Mbak Lala atas dukungan dan kesabarannya. 8. Ust. Dadang, Ust. Fatih Karim, Ust. Abdullah yang telah memberikan banyak pencerahan. 9. Lukman, Erus, Aldoko, Alam, Ari, Mas Hepy, dan Kak Ferry atas bantuan dan kesabarannya mendengarkan curhat penulis. 10. Nita, Cecep, serta rekan-rekan di www.sosmath.com atas diskusi yang telah memberikan banyak pemahaman. 11. Rekan-rekan mahasiswa Matematika 37, 38, 39, 40. 12. Seluruh dosen dan staf Departemen. 13. Semua pihak yang tidak dapat disebutkan satu persatu. Akhir kata, semoga tulisan ini dapat memberikan manfaat untuk kita semua. Kritik dan masukan yang membangun sangat penulis harapkan untuk kemajuan tulisan ini. Semoga Allah SWT senantiasa melimpahkan karunia-Nya kepada kita semua. Amin. Bogor, Mei 2007
Wahyu Hartono
RIWAYAT HIDUP Penulis dilahirkan di Jakarta pada tanggal 17 Mei 1981 sebagai anak kelima dari enam bersaudara. Ayah penulis bernama Muhammad Nurpi dan Ibu bernama Tarni. Pada tahun 1993 penulis menyelesaikan pendidikan di SD Negeri Paseban 10 Petang Jakarta Pusat. Penulis melanjutkan pendidikan di SLTP Negeri 216 Jakarta Pusat pada tahun yang sama. Pada tahun 1996 penulis melanjutkan pendidikan di SMU Negeri 77 Cempaka Putih Jakarta Pusat. Pada tahun 1999, penulis diterima di Universitas Indonesia melalui jalur Ujian Masuk Perguruan Tinggi Negeri (UMPTN) di Jurusan Biologi, Fakultas Matematika dan Ilmu Pengetahuan Alam. Pada tahun 2000, penulis diterima di Institut Pertanian Bogor melalui jalur Ujian Masuk Perguruan Tinggi Negeri (UMPTN) di Departemen Matematika, Fakultas Matematika dan Ilmu Pengetahuan Alam IPB. Selama mengikuti kegiatan perkuliahan, penulis aktif di Lembaga Dakwah Kampus Badan Kerohanian Islam Mahasiswa (LDK BKIM) IPB dan pada tahun 2004 menjabat sebagai kepala Departemen Syi’ar dan Pers BKIM IPB. Penulis juga aktif sebagai panitia pada berbagai acara antara lain “Diskusi Nasional Umat Islam” yang diselenggarakan oleh BKIM IPB pada tahun 2001, “Seminar Sains, Islam, dan Peradaban” yang diselenggarakan oleh Gugus Mahasiswa Matematika (GUMATIKA) IPB pada tahun 2004, asisten mata kuliah Kalkulus pada tahun 2004, ketua panitia Tablig Akbar dengan tema: ”Bersatu, Bergerak, Menuju Bogor Bersyariah” pada tahun 2005, panitia penggalangan ormas-ormas Islam se-Bogor untuk menolak kedatangan presiden Amerika Serikat George W. Bush ke kota Bogor pada tahun 2006, ketua panitia diskusi remaja dengan tema: ”Pergaulan Bebas; Potret Buram remaja Saat Ini” pada tahun 2007, serta berbagai kegiatan keislaman di kota Bogor selama tahun 2002 sampai 2007. Selama tahun 2001 sampai 2007, penulis aktif mengajar matematika pada lembaga pendidikan di kota Bogor serta pengajar les privat matematika dari rumah ke rumah.
DAFTAR ISI Halaman DAFTAR GAMBAR .....................................................................................................................viii DAFTAR LAMPIRAN..................................................................................................................viii I PENDAHULUAN 1.1 Latar Belakang ................................................................................................................... 1 1.2 Tujuan ................................................................................................................................ 1 1.3 Metode dan Sistematika Penulisan .................................................................................... 1 II LANDASAN TEORI 2.1 Sistem Persamaan Diferensial, nilai eigen, dan vektor eigen ............................................. 1 2.2 Titik tetap, Pelinearan dan Kestabilan................................................................................ 2 2.3 Kestabilan Lyapunov.......................................................................................................... 5 2.4 Keterbatasan, Kekontinuan, Kemotonan Fungsi dan Kekonvergenan Integral Takwajar .. 6 2.5 Solusi Persamaan Diferensial Biasa dengan Metode Operator D....................................... 7 2.6 Kontrol Optimum ............................................................................................................... 8 III MODEL-MODEL DASAR 3.1 Model Pertumbuhan Logistik............................................................................................. 9 3.2 Persamaan Lotka-Volterra.................................................................................................. 9 3.3 Fungsi Pemanenan............................................................................................................ 10 3.4 Model Keuntungan Maksimum........................................................................................ 11 3.5 Teori Modal...................................................................................................................... 11 IV PEMBAHASAN 4.1 Laju Pertumbuhan Mangsa............................................................................................... 11 4.2 Laju Pertumbuhan Pemangsa I......................................................................................... 12 4.3 Laju Pertumbuhan Pemangsa II ....................................................................................... 12 V ANALISIS KESTABILAN 5.1 Penentuan Titik Tetap ...................................................................................................... 13 5.2 Konstruksi Matriks Jacobi................................................................................................ 14 5.3 Analisis Kestabilan Titik Tetap ........................................................................................ 16 VI HASIL PEMANENAN MAKSIMUM ..................................................................................... 24 VII KEBIJAKAN PEMANENAN OPTIMAL 7.1 Fungsi Objektif Usaha Pemanenan .................................................................................. 25 7.2 Penentuan Solusi Optimal Usaha Pemanenan .................................................................. 25 7.3 Contoh Penerapan pada Industri Perikanan...................................................................... 27 VIII SIMPULAN ........................................................................................................................... 28 DAFTAR PUSTAKA ..................................................................................................................... 29 LAMPIRAN.................................................................................................................................... 30
vii
DAFTAR GAMBAR Halaman 1 Konsep Kestabilan ......................................................................................................................... 2 2 Attractor (stabil)............................................................................................................................. 4 3 Attractor (stabil)............................................................................................................................. 4 4 Repellor (takstabil)......................................................................................................................... 5 5 Repellor (takstabil)......................................................................................................................... 5 6 Titik Sadel ...................................................................................................................................... 5 7 Titik Sadel ...................................................................................................................................... 5 8 Diagram kedinamikan interaksi antarspesies ............................................................................... 13 9 Kurva solusi hasil pemanenan maksimum ................................................................................... 24 10 Dinamika mangsa di E 4 ............................................................................................................ 28 11 Dinamika pemangsa I di E 4 ...................................................................................................... 28 12 Dinamika pemangsa II di E 4 .................................................................................................... 28
DAFTAR LAMPIRAN Halaman 1 Titik Tetap.................................................................................................................................... 31 2 Nilai Eigen ................................................................................................................................... 32 3 Evaluasi Nilai-nilai Parameter...................................................................................................... 36 4 Kestabilan Global......................................................................................................................... 39 5 Populasi Optimal.......................................................................................................................... 41 6 Kurva MSY .................................................................................................................................. 54 7 Dinamika Sistem .......................................................................................................................... 56 8 Syarat Teorema 7 ......................................................................................................................... 59
viii
I PENDAHULUAN 1.1 Latar Belakang Permasalahan kebijakan pemanenan ikan yang memberikan keuntungan maksimum dan berkelanjutan (tidak terjadi kepunahan dari populasi ikan yang dipanen) adalah hal yang sangat penting bagi industri perikanan. Para ilmuwan berusaha menyelesaikan permasalahan tersebut dengan berbagai cara, salah satunya adalah dengan memodelkan permasalahan tersebut secara matematis kemudian menyelesaikannya dengan menggunakan metode yang sesuai untuk model tersebut. Dalam tulisan ini akan dipelajari dinamika populasi dari model pertumbuhan dua spesies ikan yang berkompetisi yang bergantung tidak hanya pada sumberdaya eksternal tetapi juga bergantung pada spesies ketiga, yaitu sumberdaya yang dapat memperbaharui dirinya sendiri. Selain itu, pertumbuhan dari dua spesies ikan yang berkompetisi juga dipengaruhi oleh usaha pemanenan terhadap kedua spesies ikan tersebut. Permasalahan yang ditampilkan dalam tulisan ini diformulasikan dalam bentuk kontrol optimum untuk unit waktu takterbatas (infinite horizon). Calon solusi optimal diperoleh dengan menerapkan Prinsip Maksimum Pontryagin. Karya ilmiah ini merupakan rekonstruksi dari tulisan Chattopadhyay dan kawan-kawan yang berjudul A resource based competitive system in three species fishery.
1.2 Tujuan Tujuan penulisan karya ilmiah ini adalah: 1. melakukan analisis kestabilan titik tetap dan mengamati dinamika populasi ikan, 2. menganalisis model kebijakan pemanenan ikan yang dapat memberikan keuntungan maksimum dan berkelanjutan. 1.3 Metode dan Sistematika Penulisan Metode penulisan karya ilmiah ini adalah studi literatur. Karya ilmiah ini terdiri atas delapan bagian. Bagian pertama menjelaskan latar belakang masalah, tujuan, dan sistematika penulisan. Bagian kedua menyajikan landasan teori yang membahas sistem dinamika, teori modal, dan kontrol optimum. Bagian ketiga menyajikan modelmodel dasar yang membangun model yang akan dianalisis. Bagian keempat menyajikan pembahasan masalah dengan merekonstruksi model yang akan dianalisis. Bagian kelima menyajikan analisis kestabilan. Bagian keenam membahas masalah pemanenan maksimum. Bagian ketujuh menyajikan formulasi kebijakan pemanenan optimal yang dibuat dalam bentuk formulasi masalah kontrol optimum dan pembahasan penentuan calon solusi optimal dengan menerapkan prinsip maksimum Pontryagin. Bagian terakhir menyajikan kesimpulan tentang hasil yang diperoleh dari pembahasan sebelumnya.
II LANDASAN TEORI Beberapa landasan teori yang dibahas pada bab ini meliputi sistem persamaan diferensial (SPD), titik tetap, pelinearan dan konsep kestabilan, serta teori kontrol optimum yang disarikan dari berbagai sumber pustaka. 2.1 Sistem persamaan diferensial (SPD), nilai eigen, dan vektor eigen Definisi 1. (SPD Mandiri) Misalkan suatu sistem persamaan diferensial (SPD) dinyatakan sebagai dx = x = f ( x), x ∈ n dt dengan f fungsi kontinu dari x dan mempunyai turunan parsial kontinu. SPD tersebut disebut sistem persamaan diferensial
mandiri jika tidak mengandung t secara eksplisit di dalamnya. (Kreyszig, 1983) Definisi 2. (SPD Linear Orde Satu) Jika suatu sistem persamaan diferensial (SPD) dinyatakan sebagai x = Ax + b, x(0) = x0 , x ∈ R n
(2.1.1)
dengan A adalah matriks koefisien berukuran n × n , maka sistem tersebut dinamakan SPD linear orde 1 dengan kondisi awal x(0) = x0 . Sistem (2.1.1) disebut homogen jika b = 0 dan takhomogen jika b ≠ 0 . (Fisher, 1990)
1
2
Definisi 3. (Vektor Eigen dan Nilai Eigen) Jika A adalah suatu matriks n × n , maka n disebut suatu vektor taknol x pada vektor eigen dari A jika A x adalah suatu penggandaan skalar dari x ; yaitu, Ax = λx (2.1.2)
untuk suatu skalar λ . Skalar λ disebut nilai eigen dari A , dan x disebut suatu vektor eigen dari A yang berpadanan dengan λ . (Anton, 2000) Untuk mencari nilai λ dari matriks A , maka persamaan (2.1.2) dapat ditulis kembali sebagai ( A − λI ) x = 0 (2.1.3) dengan I matriks identitas. Persamaan (2.1.3) mempunyai solusi taknol jika dan hanya jika p (λ ) = det( A − λ I ) = A − λ I = 0 (2.1.4) Persamaan (2.1.4) disebut karakteristik dari matriks A .
persamaan (Anton, 2000)
2.2 Titik tetap, pelinearan dan kestabilan Definisi 4. (Titik Tetap dan Titik Biasa) Misalkan diberikan SPD taklinear berikut x = f ( x) : U ⊂
n
→
Definisi 6. (Titik Tetap Stabil Asimtotik Lokal) Titik x∗ dikatakan titik tetap stabil asimtotik lokal jika titik x∗ stabil dan terdapat ε > 0 sedemikian sehingga jika x∗ − x0 < ε
maka lim x ( t ) = x∗ , dengan x ( 0 ) = x0 . t →∞
(Szidarovszky & Bahill, 1998) Definisi 7. (Titik Tetap Stabil Asimtotik Global) Titik x∗ dikatakan titik tetap stabil asimtotik global jika titik x∗ stabil dan ∗ n x0 ∈ Ω ⊆ R , lim x ( t ) = x , dengan t →∞
x ( 0 ) = x0 .
(Szidarovszky & Bahill, 1998)
ε0
Takstabil Stabil
ε1 ε x
∗
Stabil asimtotik lokal
(2.2.1)
n
Titik x∗ dengan f ( x∗ ) = 0 disebut titik kritis atau titik tetap atau titik keseimbangan. Titik x pada bidang fase (2.2.1) yang bukan titik tetap, yaitu f ( x ) ≠ 0 , disebut titik biasa atau titik regular.
Stabil asimtotik global Gambar 1 Konsep Kestabilan dari Definisi 5, 6, dan 7.
(Tu, 1994) Definisi 5. (Titik Tetap Stabil) Misalkan x∗ adalah titik tetap sebuah SPD mandiri dan x ( t ) adalah solusi yang
memenuhi kondisi awal x ( 0 ) = x0 dengan ∗
∗
x0 ≠ x . Titik x dikatakan titik tetap stabil jika terdapat ε 0 > 0 yang memenuhi sifat berikut: untuk setiap ε1 , dengan 0 < ε1 < ε 0 , terdapat
ε > 0 sedemikian sehingga jika x − x0 < ε ∗
maka x∗ − x ( t ) < ε1 , untuk setiap t > t0 . (Szidarovszky & Bahill, 1998) Catatan: Norm x = x12 + x22 + ... + xn2
Definisi 8. (Titik Tetap Takstabil) Misalkan x∗ adalah titik tetap sebuah SPD mandiri dan x ( t ) adalah solusi yang
memenuhi kondisi awal x ( 0 ) = x0 dengan x0 ≠ x∗ . Titik x∗ dikatakan titik tetap takstabil jika terdapat ε 0 > 0 yang memenuhi sifat berikut: untuk sembarang ε , 0 < ε < ε 0 , terdapat x0
sedemikian sehingga jika x∗ − x0 < ε maka x∗ − x ( t ) ≥ ε 0 , untuk setiap t > t0 .
(Velhurst, 1990)
3
λ dengan bagian-bagian real nol, berarti Re ( λ ) ≠ 0 .
Definisi 9. (Titik Tetap Simple) Suatu titik tetap SPD taklinear dikatakan simple jika pada sistem linearnya ( Ax ) tidak
mempunyai nilai jika det A ≠ 0 .
eigen
nol,
(Tu, 1994) Teorema 1. (Teorema Taylor untuk n variabel) Misalkan f : Ω ⊂ n → dan segmen garis yang menghubungkan x0 dan x berada pada Ω.
yaitu
(Tu, 1994) Definisi 10. Titik Tetap Simple Hyperbolic Titik tetap simple dikatakan hyperbolic jika A pada Ax tidak memiliki nilai-nilai eigen
(
Jika x = ( x1 ,..., xn ) dan x0 = x1( ) ,..., xn ( 0
0)
)
maka
(
n
f ( x ) = f ( x0 ) + ∑ f i ( x 0 ) xi − xi ( i =1
(
1 n + ∑ f ijk ( x 0 ) xi − xi ( 0) 3! i , j , k =1
(x
i2
− xi2 (
0)
)( x
im
− xim (
0)
)( x
)+ R
m
j
0)
) + 2!1 ∑∑ f
− x j ( 0)
n
n
i =1 j =1
)( x
k
ij
( x0 ) ( xi − xi (0) ) ( x j − x j (0) )
)
− xk ( 0) + ... +
(
)
n 1 f i ,i ,...,i ( x0 ) xi1 − xi1 ( 0) ... ∑ m ! i1 ,i2 ,...,im =1 1 2 m
(x)
( 2.2.2 )
dengan Rm ( x ) =
(
n 1 0 f i1 ,i2 ,...,im+1 ( x 0 + ch ) xi1 − xi1 ( ) ∑ ( m + 1)! i1 ,i2 ,...,im+1 =1
)( x
i2
− xi2 (
0)
) ...( x
im+1
− xim+1 (
0)
)
(2.2.3)
untuk nilai c pada ( 0,1) dan h = x − x0 . (Grossman, 1995) Dengan memerhatikan sistem taklinear x = f ( x) pada (2.2.1) dan menggunakan
dengan
∗
ekspansi Taylor di sekitar titik tetap x , misalkan di titik asal ( 0, 0 ) untuk penyederhanaan, diperoleh x = Ax + ϕ ( x ) (taklinear)
(2.2.4)
dengan
lim
A ≡ Df ( x 1
⎢ ∂x1 ≡ ⎢⎢ ⎢ ∂f ( x∗ ) ⎢ n ⎢⎣ ∂x1 dan ϕ ( x)
lim ϕ ( x ) = 0 .
b≡
∂f1 ≡ a12 ∂x2
c≡
∂f 2 ≡ a21 ∂x1
d≡
∂f 2 ≡ a22 ∂x2
ϕ1 ( x1 , x2 )
r →0
)≡ ⎡ ∂f ( x ) ⎢
∂f1 ≡ a11 ∂x1
r
= lim
ϕ 2 ( x1 , x2 )
r →0
r
dengan
∗
∗
x →0
dan
a≡
Ax
∂f1 ( x∗ ) ⎤ ⎥ ∂xn ⎥ ⎡ a11 a1n ⎤ ⎥ ⎥≡⎢ ⎥ ⎥ ⎢ ann ⎥⎦ ∂f n ( x∗ ) ⎥ ⎢⎣ an1 ⎥ ∂xn ⎥⎦ sedemikian sehingga
pada
(2.2.4)
disebut
pelinearan dari (2.2.4) atau sistem yang dilinearkan yaitu x = Ax. (linear) (2.2.5) Untuk sistem pada bidang, U ⊂ R 2 diperoleh x = f ( x ) = Ax + ϕ ( x ) sebagai x1 = f1 ( x ) = ax1 + bx2 + ϕ1 ( x1 , x2 ) x2 = f 2 ( x ) = cx1 + dx2 + ϕ 2 ( x1 , x2 )
r = x12 + x22 ,
yaitu ϕ1 dan ϕ2 menjadi sangat kecil dan bisa diabaikan. Secara umum, untuk setiap titik tetap ∗ ∗ ( x , y ) ≠ ( 0, 0 ) , dapat didefinisikan variabelvariabel baru ξ ≡ ( x − x∗ ) dan η ≡ ( y − y ∗ ) ,
sehingga x = f ( x ) adalah ⎛ x ⎞ ⎛ ξ ⎞ ⎡ a b ⎤ ⎡ x − x∗ ⎤ ⎡ϕ1 (ξ ,η ) ⎤ +⎢ ⎥ ( 2.2.6 ) ⎜ ⎟≡⎜ ⎟=⎢ ⎥⎢ ∗⎥ ⎝ y ⎠ ⎝η ⎠ ⎣ c d ⎦ ⎣ y − y ⎦ ⎣ϕ 2 (ξ ,η ) ⎦
Persamaan (2.2.6) hanya menggambarkan perpindahan titik tetap dari ( 0, 0 ) ke titik
tetap taknol ( x∗ , y ∗ ) .
(Tu, 1994)
4
Matriks A pada (2.2.4) disebut sebagai matriks Jacobi atau biasa juga disebut sebagai matriks variasi. (Abell & Braselton, 2004) Teorema 2. (Teorema Pelinearan dari Hartman & Grobman) Misalkan sistem dinamik taklinear x = f ( x) pada (2.2.1) mempunyai titik tetap simple hyperbolic x∗ , dan misalkan titik tetap tersebut berada pada ( 0, 0 ) untuk penyederhanaan, maka pada U di sekitar x∗ ∈ n , bidang fase dari sistem taklinear (2.2.4) dan sistem hasil pelinearannya (2.2.5) adalah ekuivalen. (Tu, 1994) Misalkan λ j = α j ± i β j adalah nilai eigen dari matriks Jacobi A dengan Re λ j = α j dan Im ( λ j ) = β j . Kestabilan
( )
titik tetap simple hyperbolic mempunyai beberapa perilaku sebagai berikut: 1. Stabil Attractor, jika : a. Re λ j = α j < 0 dengan Im ( λ j ) = 0
( )
(2.2.7) untuk setiap j = 1, 2,3 . Lihat Gambar 2. b. Terdapat Re ( λ j ) = α j < 0 dengan Im ( λ j ) = 0 dan Re ( λk ) = α k < 0
dengan Im ( λk ) ≠ 0 untuk j = 1, 2,3 dan j≠k. Lihat Gambar 3.
(2.2.8)
2. Takstabil Repellor, jika : a. Re ( λ j ) = α j > 0 dengan Im ( λ j ) = 0 untuk setiap j = 1, 2,3. Lihat Gambar 4.
Gambar 2 Attractor (stabil).
b. Terdapat Re ( λ j ) = α j > 0 dengan Im ( λ j ) = 0 dan Re ( λk ) = αk > 0
dengan Im ( λk ) ≠ 0 untuk j = 1, 2,3 dan j≠k. Lihat Gambar 5.
(2.2.10)
3. Titik Sadel, jika : a. Terdapat Re λ j = α j > 0 dengan
( )
Im ( λ j ) = 0 dan Re ( λk ) = αk < 0
dengan Im ( λk ) ≠ 0 untuk j = 1, 2,3 dan j≠k. (2.2.11) Lihat Gambar 6. b. Terdapat Re λ j = α j > 0 dengan
( )
Im ( λ j ) = 0 dan Re ( λk ) = αk < 0
dengan Im ( λk ) = 0 untuk j = 1, 2,3 dan j≠k. (2.2.12) Lihat Gambar 7. c. Terdapat Re λ j = α j < 0 dengan
( )
Im ( λ j ) = 0 dan Re ( λk ) = α k > 0
dengan Im ( λk ) ≠ 0 untuk j = 1, 2,3 dan j≠k. (2.2.13) Lihat Gambar 8. d. Terdapat Re λ j = α j < 0 dengan
( )
Im ( λ j ) = 0 dan Re ( λk ) = αk > 0
dengan Im ( λk ) = 0 untuk j = 1, 2, 3 dan j≠k. Lihat Gambar 9.
(2.2.14) (Tu, 1994)
(2.2.9)
Gambar 3 Attractor (stabil).
5
Gambar 4 Repellor (takstabil).
Gambar 7 Titik Sadel.
Gambar 5 Repellor (takstabil).
Gambar 8 Titik Sadel.
Gambar 6 Titik Sadel.
Gambar 9 Titik Sadel.
Gambar 2 sampai 9 merupakan aliran-aliran hyperbolic dalam 3 dimensi.
2.3 Kestabilan Liapunov Dalam Matematika, kestabilan Liapunov sering digunakan untuk mempelajari masalah dinamika sistem. Berikut ini adalah definisi-definisi formal yang digunakan untuk membangun kestabilan Liapunov. Definisi 11. (Matriks Definit Positif) Suatu matriks simetriks real A disebut definit positif jika xT Ax > 0 untuk semua x taknol dalam n . (Leon, 2001)
Teorema 3 Misalkan A adalah matriks simetrik real berorde n × n . Maka A adalah definit positif jka dan hanya jika semua nilai-nilai eigennya adalah positif. (Leon, 2001) Definisi 12. (Fungsi Liapunov) Suatu fungsi bernilai real dan terturunkan V ( x ) pada titik x di sekitar x , untuk penyederhanaan misalkan x = 0 , sedemikian sehingga V (x) ≥ 0 , V ( 0) = 0 yaitu
6
V ( x ) = 0 jika dan hanya jika x = x ( = 0 ) ,
disebut fungsi Liapunov. (Tu, 1994) Definisi 13. (Fungsi Definit Positif/Negatif dan Semidefinit Positif/Negatif) Misalkan V (x) adalah fungsi bernilai real dengan x adalah ( x1 , x2 ,..., xn ) . V ( x ) adalah definit positif (negatif) pada daerah S di sekitar titik asal jika V ( x ) > 0 (V ( x ) < 0) untuk semua x ≠ 0 pada S, dan V(0) = 0 .
2.4 Keterbatasan, Kekontinuan, Kemonotonan Fungsi dan Kekonvergenan Integral Takwajar Definisi 14. (Fungsi Terbatas) Suatu fungsi f : A → B disebut terbatas jika terdapat bilangan real M > 0 sedemikian sehingga f ( x ) ≤ M untuk setiap x ∈ A . (Zipse, 1987) Definisi 15. ( Fungsi Kontinu) Suatu fungsi f dikatakan kontinu di c jika
V ( x ) adalah semidefinit positif (negatif)
pada daerah S di sekitar titik asal jika V ( x ) ≥ 0 (V ( x ) ≤ 0) untuk semua x ≠ 0 pada S, dan V(0) = 0 . ( Smith, 1987) Teorema 4. (Stabil Seragam) Misalkan x* ( t ) , t ≥ t0 adalah solusi nol dari sistem x = X(x) , dan X(0) = 0 , maka x* ( t ) adalah stabil seragam untuk
t ≥ t0 jika
terdapat fungsi bernilai real V ( x ) pada daerah di sekitar x = 0 dengan : (i) V ( x ) dan turunan-turunan parsialnya kontinu; ( ii ) V ( x ) adalah definit positif; ( iii ) V ( x ) adalah semidefinit negatif. Teorema 5. (Stabil Asimtotik Lokal dan Global) Misalkan semua kondisi pada Teorema 4 dipenuhi, kecuali ( iii ) diganti oleh ( iii)’ V ( x ) adalah definit negatif
2.
f ( c ) ada.
3.
lim f ( x ) = f ( c ) .
x →c
x →c
Definisi 16. (Kontinu Bagian demi Bagian) Fungsi yang mempunyai sejumlah berhingga titik ketakkontinuan disebut fungsi yang kontinu bagian demi bagian (piecewise continuity). (Stewart, 1991) Definisi 17. (Kemonotonan) Misalkan f terdefinisi pada selang I (terbuka, tertutup, atau tak satupun). Fungsi f dikatakan i. naik pada I jika untuk setiap pasang bilangan x1 dan x2 dalam I , x1 < x2 ⇒ f ( x1 ) < f ( x2 )
ii.
turun pada I jika untuk setiap pasang bilangan x1 dan x2 dalam I , x1 < x2 ⇒ f ( x1 ) > f ( x2 )
Catatan:
disebut fungsi Liapunov lemah, sedangkan fungsi yang memenuhi Teorema 6 disebut fungsi Liapunov kuat. ( Smith, 1987)
lim f ( x ) ada.
Jika salah satu dari ketiga syarat tersebut takterpenuhi maka f dikatakan takkontinu. Titik dimana fungsi f takkontinu disebut titik ketakkontinuan (discontinuity) dari f . (Stewart, 1991)
maka solusi nol-nya adalah stabil asimtotik.
1. Daerah di sekitar titik asal disebut domain kestabilan asimtotik, atau domain atraksi. 2. Ketika domain atraksinya adalah seluruh bidang x , maka sistemnya disebut stabil asimtotik global. 3. Fungsi V ( x ) yang memenuhi Teorema 5
1.
(Purcell, 1999) Teorema 6 Misalkan fungsi f kontinu pada [ a, b ] dan terturunkan pada ( a, b ) . f ′ ( x ) > 0 untuk setiap x pada
i.
Jika
ii.
( a, b ) , maka f Jika f ′ ( x ) < 0 ( a, b ) , maka f
naik pada [ a, b ] . untuk setiap x pada turun pada [ a, b ] . (Purcell, 1999)
7
Definisi 18. (Integral Takwajar) t
Jika
∫ f ( x ) dx ada untuk setiap
t ≥ a , maka
a
∞
∫ a
t
f ( x ) dx = lim ∫ f ( x ) dx t →∞
(2.4.1)
a
dengan syarat limitnya ada (sebagai bilangan berhingga). Jika limitnya ada maka integral pada (2.4.1) disebut konvergen. (Stewart, 1991) 2.5 Solusi Persamaan Diferensial Biasa dengan Metode Operator D Persamaan diferensial berorde n dengan koefisien-koefisien konstan dapat dituliskan sebagai an y n + an −1 y n −1 + … + a1 y + a0 = Q . (2.5.1) d sehingga persamaan (2.5.1) dx dapat dituliskan kembali sebagai persamaan (2.5.2) berikut ( an D n + an −1 D n −1 + … + a1 D + a0 ) y = Q.
Misalkan D ≡
dengan Q adalah fungsi dari x . Solusi umum dari persamaan (2.5.2) berbentuk y = CF + PI dengan CF (Complementary Function) mengandung n konstanta sembarang dan PI (Particular Integral) yang tidak mengandung konstanta sembarang. Persamaan (2.5.2) dapat ditulis dalam bentuk f ( D) y = Q , (2.5.3) dengan f ( D ) adalah polinomial berderajat n . Untuk kasus Q = 0 , dengan menyusun persamaan bantu (Auxiliary Equation) diperoleh solusi dari persamaan (2.5.3) yang biasa disebut CF . Sedangkan untuk 1 Q adalah menentukan PI dimisalkan f ( D) fungsi dari x
(Q ≠ 0)
sedemikian sehingga
⎛ 1 ⎞ f ( D ) ⎜⎜ Q ⎟⎟ = Q , ⎝ f ( D) ⎠
(2.5.4)
maka y=
1 Q f ( D)
sebagai operator yang dioperasikan pada Q , untuk menghasilkan PI . Untuk kasus Q = eα x dengan α adalah konstanta, aturan turunan memberikan hasil sebagai berikut DQ = Deα x = α eα x (2.5.6)
( D + a ) eα x = Deα x + aeα x = α eα x + aeα x = (α + a ) eα x .
(2.5.7)
Dari (2.5.6) dan (2.5.7) dapat disimpulkan f ( D ) eα x = f (α ) eα x , sehingga 1 1 eα x = eα x , f (α ) ≠ 0 . f ( D) f (α ) (Bajpai, 1970) Contoh. Carilah solusi umum dari persamaan diferensial 3 y ′′ − y ′ − 4 y = 5e 2 x . (2.5.8) Jawab. CF diperoleh dari 3 y ′′ − y ′ − 4 y = 0 persamaan bantunya adalah 3β 2 − β − 4 = 0 4 , adalah akar-akar 3 persamaan bantu. Oleh karena itu diperoleh
dan
β1 = −1, β 2 =
CF = A1e β1 + A2 e β2 4
x
= A1e − x + A2 e 3 ,
konstanta
sembarang.
dengan Misalkan
A1 , A2 D≡
d , dx
maka persamaan (2.5.8) menjadi ( 3D 2 − D − 4 ) y = 5e2 x . PI diperoleh dari 5e2 x y= ( 3D 2 − D − 4 )
= 5.
1 3( 2) − ( 2) − 4 2
e2 x
5 2x e . 6 Jadi solusi umum dari persamaan diferensial (2.5.8) adalah y = CF + PI =
(2.5.5)
adalah solusi PI dari persamaan (2.5.3). 1 pada persamaan (2.5.5) bertugas f ( D)
4 x 5 = A1e − x + A2 e 3 + e 2 x , dengan A1 , A2 6 konstanta sembarang.
8
∞
2. 6 Kontrol Optimum Dalam suatu masalah ekonomi yang berkembang menurut waktu, sistem pada waktu t dapat diungkapkan dengan peubah keadaan (state variabel) x1 ( t ) ,..., xn (t ) atau dalam bentuk vektor x ( t ) ∈
. Dengan nilai
t yang berbeda, vektor x ( t ) menempati
posisi yang berbeda di ruang . Dinamika sistem dapat dinyatakan secara matematik melalui persamaan diferensial: (2.6.1) x = f ( x ( t ) ,U ( t ) , t ) dx . Misalkan diketahui keadaan dt sistem pada waktu t0 , x ( t0 ) = x0 ∈ . Jika
dengan x =
dipilih
peubah
kontrol
U (t ) ∈
yang
terdefinisi untuk t ≥ t0 , maka diperoleh persamaan diferensial orde satu dengan peubah taktentu x ( t ) . Karena x0 diberikan maka (2.6.1) mempunyai solusi tunggal. Solusi yang diperoleh merupakan respon terhadap kontrol U yang dilambangkan dengan xU ( t ) . Dengan memiliki fungsi kontrol yang sesuai, berbagai solusi dapat diperoleh. Agar solusi yang diperoleh adalah solusi yang diinginkan, diperlukan adanya kriteria bagi solusi yang diinginkan, artinya untuk setiap kontrol U ( t ) dan responnya state
maks u (.) ∫ f 0 ( t , x, u ) dt , terhadap kendala x = f ( x, u ) v.e.
J (U ) =
∫ f ( x (t ) ,U (t ) , t ) dt 0
(2.6.2)
t0
dengan f 0 fungsi yang diberikan. T tidak
harus fixed (ditentukan) dan x (T ) mempunyai
kondisi tertentu. Di antara semua fungsi/peubah kontrol yang diperoleh ditentukan salah satu sehingga J menjadi maksimum. Kontrol yang bersifat demikian disebut kontrol optimum. (Tu, 1993) Berikut ini akan dibahas prinsip maksimum Pontryagin untuk masalah kontrol optimum dengan horizon waktu takberhingga. Teorema 7. Misalkan diberikan seperti berikut
suatu
permasalahan
(virtually
everywhere),
u ∈ U , x ( 0 ) = x0 ,
(2.6.4)
lim xi ( t ) = x ,
(2.6.5)
i
t →∞
i = 1,… , m .
U adalah himpunan yang berada di
m
,
∂fi , ∂x i = 0,… , n , kontinu, demikian juga dengan ∂f 0 , dan v.e. (virtually everywhere) berarti ∂t terdefinisi di mana-mana, kecuali di sejumlah berhingga titik. Diasumsikan bahwa ( x∗ ( t ) , u ∗ ( t ) ) fi :
n
→ ,
i = 0,… , n ,
dan
fi
adalah solusi optimal di antara semua pasangan ( x (.) , u (.) ) yang memenuhi (2.6.4), (2.6.5), dan integral pada (2.6.3) konvergen, dengan u (.) kontinu bagian demi bagian pada setiap interval. Diasumsikan juga bahwa ∞
f 0 (t + σ , x ∗ (t ) , u ∗ (t )) d (t )
∫ 0
setiap σ
ada
( −ε , ε ) ,
pada interval
untuk
sehingga
untuk fungsi α ( t ) yang kontinu bagian demi bagian, berlaku
∂f 0 ( t + σ , x∗ ( t ) , u ∗ ( t ) ) ∂t
x ( t ) dihubungkan dengan fungsi T
(2.6.3)
0
≤ α ( t ) untuk setiap ∞
δ ∈ ( −ε , ε ) , t ∈ [ 0, ∞ ) dan
∫
α (t )dt < ∞ .
0
Maka terdapat pasangan
( p , p (t )) ≠ 0 0
( p , p (t )) , 0
p0 ≥ 0,
untuk setiap t , sedemikian
sehingga p ( t ) adalah solusi kontinu dari persamaan adjoint p = − H x ( t , x∗ ( t ) , u ∗ ( t ) , p ( t ) ) v.e.,
(2.6.6)
dengan H = p0 f 0 ( t , x, u ) + p ( t ) f ( x, u ) . Selanjutnya maks u H ( t , x∗ ( t ) , u ( t ) , p ( t ) )
= H ( t , x∗ ( t ) , u ∗ ( t ) , p ( t ) ) v.e.
dan lim H ( t , x∗ ( t ) , u ∗ ( t ) , p ( t ) ) = 0. t →∞
(2.6.7) (2.6.8)
(Seierstad, 2002)
III MODEL-MODEL DASAR
Analisis model kebijakan pemanenan optimal dalam perikanan, memerlukan modelmodel dasar yang mendukung dibangunnya model kebijakan pemanenan tersebut. Oleh karena itu, dalam bab ini akan ditampilkan model-model dasar dan teori pembangun dari model kebijakan pemanenan optimal perikanan yang menjadi topik utama dalam tulisan ini. 3.1 Model Pertumbuhan Logistik Menurut Shone (1997), jika jumlah populasi pada waktu t disebut stok, tingkat stok ini akan berubah bergantung pada perbedaan antara arus masuk dan arus keluar. Contoh dari arus masuk adalah kelahiran dan imigrasi, sedangkan arus keluar adalah kematian dan emigrasi. Pada kasus populasi ikan terdapat faktor penangkapan ikan yang dilakukan oleh manusia dalam suatu periode tertentu yang akan mempengaruhi tingkat stok. Dari faktor ’alamiah’ yang mempengaruhi populasi didapatkan persamaan
Perubahan netto dalam populasi = arus masuk – arus keluar = (kelahiran + imigrasi) – (kematian + emigrasi) Kelahiran dan kematian disebut sebagai faktor internal dari populasi, sedangkan imigrasi dan emigrasi disebut sebagai faktor eksternal. Persamaannya menjadi Perubahan netto dalam populasi = perubahan internal + perubahan eksternal = (kelahiran – kematian)+ migrasi dengan migrasi adalah imigrasi dikurangi emigrasi. Misalkan n ( t ) merupakan variabel yang memberikan kontribusi pada perubahan internal. Ukuran populasi pada suatu waktu adalah x(t ) , yang menotasikan banyaknya individu pada waktu t , sehingga perubahan internal dari populasi adalah n ( t )x(t ) . Misalkan juga m(t ) menotasikan migrasi yang terjadi pada suatu interval waktu yang sama dengan pengukuran perubahan internal, dan diukur pada waktu t , maka m(t ) menotasikan perubahan eksternal, sehingga perubahan populasi dx(t ) / dt diberikan oleh
dx = n ( t ) x(t ) + m(t ) (3.1.1) dt Jika diasumsikan tidak ada migrasi maka m(t ) = 0 untuk setiap t . Jika diasumsikan juga ada pengurangan dalam proses pertumbuhan populasi yang proporsional terhadap ukuran populasi, dengan kata lain, laju pertumbuhan r direduksi oleh faktor ax(t ) maka variabel yang memberikan kontribusi pada perubahan internal pada populasi menjadi n(t ) = r − ax(t ) . Dari asumsi tentang migrasi dan perubahan internal, maka persamaan (3.1.1) dapat dituliskan sebagai dx = (r − ax(t )) x(t )dt r ⎛ x(t ) ⎞ = r ⎜1 − ⎟ , dengan K = a (3.1.2) K ⎝ ⎠ yang merupakan persamaan logistik dengan r , K adalah parameter positif, r merupakan laju pertumbuhan intrinsik yang berhubungan dengan kelahiran dan kematian alami populasi. Laju kelahiran per kapita adalah ⎛ x(t ) ⎞ yang bergantung pada x(t ) . r ⎜1 − K ⎟⎠ ⎝
r menyatakan daya muat a lingkungan (carrying capacity) yang menyatakan kapasitas maksimum populasi dalam lingkungan tersebut. Hal ini berarti jika di dalam populasi ada x individu, maka lingkungan masih bisa memuat K − x individu.
Parameter K =
3.2 Persamaan Lotka -Volterra 3.2.1 Persaingan Interspesifik Jika ada 2 spesies populasi y dan z dalam satu lingkungan yang berinteraksi secara negatif (terjadi persaingan), maka persamaannya adalah ry ( K y − y − α yz z ) dy (3.2.1.1) = y dt Ky rz ( K z − z − α zy y ) dz = z dt Kz
α yz =
a yz a yy
,
(3.2.1.2) (3.2.1.3)
9
10
α zy =
azy
azz dengan
(3.2.1.4)
α yz = koefisien interaksi dari populasi
z
pada populasi y , α zy = koefisien interaksi dari populasi y pada populasi z , a yz = laju pertumbuhan per kapita populasi y terkait dengan penambahan individu dari populasi z , a zy =laju pertumbuhan per kapita
dari satu dari
populasi z terkait dengan penambahan satu individu dari populasi y , dalam populasi y a yy =laju intrinsik r terkait dengan penambahan satu individu dari populasi y , azz =laju intrinsik r dalam populasi z terkait dengan penambahan satu individu dari populasi z , ry =laju intrinsik populasi y , rz =laju intrinsik populasi z , K y =daya muat lingkungan untuk populasi y, K z =daya muat lingkungan untuk populasi z . (Vandermeer, 1981)
3.2.2 Mangsa-Pemangsa Jika ada 2 spesies populasi x dan y dalam satu lingkungan dengan x adalah populasi mangsa dan y adalah populasi pemangsa dengan asumsi bahwa satu-satunya sumber kematian bagi populasi mangsa x adalah dimakan oleh pemangsa y dan satusatunya sumber reproduksi bagi populasi pemangsa y adalah dengan memakan mangsa x , maka interaksi antara kedua populasi digambarkan dalam persamaan dx (3.2.2.1) = bx − axy y x dt dy (3.2.2.2) = a yx x − d y y dt yang dikenal sebagai persamaan mangsapemangsa Lotka –Volterra dengan : x =banyaknya populasi mangsa, y =banyaknya populasi pemangsa, bx =laju kelahiran mangsa x tanpa kehadiran pemangsa y ,
a xy y =laju
kematian
mangsa
dimangsa oleh pemangsa y , a yx x =laju kelahiran pemangsa
karena
x y
untuk
setiap mangsa x yang dimangsa, d y =laju kematian pemangsa y karena ketiadaan mangsa x . (Vandermeer, 1981) 3.3 Fungsi Pemanenan Menurut (Shone, 1997) dan (Jaeger, 1997), fungsi pemanenan ikan adalah bentuk dari fungsi produksi yang berarti penangkapan pada waktu t , dinotasikan dengan output T sebagai fungsi dari input-input. Diasumsikan inputnya ada dua jenis. Pertama, stok ketersediaan ikan untuk ditangkap, s (t ) ; dan kedua, usaha yang dilakukan oleh nelayan, U (t ) . Fungsi pemanenan ikan yang digunakan dalam literatur adalah T = ηUs , η > 0 dengan η menotasikan efisiensi teknis dari
pemanenan/koefisien ketertangkapan, U (t ) adalah banyaknya perjalanan yang dilakukan oleh kapal penangkap ikan (atau banyaknya hari dioperasikannya jaring-jaring ikan). 3.4 Model Keuntungan Maksimum Menurut Nicholson (2002) dan Nababan (1988), sebuah perusahaan yang diasumsikan memaksimumkan keuntungan, akan memilih input dan outputnya dengan tujuan tunggal untuk mencapai keuntungan ekonomi maksimum. Para manajer perusahaan berusaha untuk menjadikan selisih antara total pendapatan dan total biaya ekonomi sebesar mungkin dengan membuat keputusan secara marjinal. Dalam kegiatannya, perusahaan menjual tingkat output tertentu s , dengan harga pasar p per unit. Maka penerimaan total/Total Revenue ( TR ) dapat diketahui, yaitu TR ( s ) = p ( s ) s .
Dalam memproduksi s , total biaya ekonomi yang ditanggung dapat dinyatakan dengan Total Cost (TC (s )) . Selisih antara pendapatan dan biaya disebut laba ekonomi π yang bergantung pada jumlah output yang diproduksi, yaitu π ( s ) = TR ( s ) − TC ( s ) .
11
Kondisi yang diperlukan untuk nilai s yang memaksimumkan laba ditentukan dengan dπ dTR dTC = π , (s) = − =0 ds ds ds sehingga kondisi agar laba maksimum adalah dTR dTC = . ds ds Persamaan ini menyatakan bahwa marginal revenue (MR) sama dengan Marginal Cost (MC), yaitu: MR = MC 3.5 Teori Modal Apabila modal awal P diinvestasikan dan pada suku bunga tahunan δ dimajemukan k kali per tahun, maka setelah t tahun nilainya akan bertambah secara eksponensial seperti formula berikut : kt
⎛ δ⎞ B ( t ) = P ⎜1 + ⎟ . (3.5.1) ⎝ k⎠ Jika suku bunga dimajemukan secara kontinu (bukan tahunan, bulanan, atau harian, tetapi setiap saat) dengan kata lain k → ∞ , maka formula (3.5.1) akan berubah dengan perhitungan sebagai berikut: k maka n → ∞ dan k = nδ misal n =
δ
kt
⎛ δ⎞ ⎛ 1⎞ P ⎜1 + ⎟ = P ⎜1 + ⎟ ⎝ k⎠ ⎝ n⎠
nδ t
δt
⎡⎛ 1 ⎞ n ⎤ = P ⎢⎜ 1 + ⎟ ⎥ . ⎢⎣⎝ n ⎠ ⎥⎦
Karena n → ∞ , maka δt
⎡⎛ 1 ⎞ n ⎤ B ( t ) = lim P ⎢⎜ 1 + ⎟ ⎥ = Peδ t (3.5.2) n →∞ ⎢⎝ n ⎠ ⎥⎦ ⎣ dengan B ( t ) adalah future value, waktu t ≥ 0
dan e ≈ 2.718 . Present value dari penerimaan sebesar B di masa yang akan datang adalah P = B ( t ) e −δ t (3.5.3) Total present value yang kontinu dari deretan pendapatan P > 0 , dengan 0 ≤ t ≤ T dapat diformulasikan sebagai berikut T
∫
P = B ( t ) e −δ t dt
(3.5.4)
0
Untuk kasus T → ∞ , persamaan (3.5.4) dimodifikasi menjadi: T
P = lim
T →∞
∫ B (t ) e
−δ t
dt
0
∞
∫
= B ( t ) e−δ t dt
(3.5.5)
0
sehingga
(Hoffmann and Bradley, 1989)
IV PEMBAHASAN Model satu mangsa dua pemangsa dengan dua pemangsa terjadi hubungan kompetisi 4.1 Laju pertumbuhan mangsa Misalkan x merupakan banyaknya individu mangsa pada waktu t , y merupakan banyaknya individu pemangsa I pada waktu t , dan z merupakan banyaknya individu pemangsa II pada waktu t , dengan individu mangsa adalah makanan satu-satunya bagi kedua pemangsa, dan hadirnya para nelayan yang memangsa pemangsa I dan II maka interaksi antarspesies dapat dimodelkan sebagai berikut : Dari persamaan (3.1.2), laju pertumbuhan per kapita dari mangsa dapat dinyatakan sebagai berikut: 1 dx x⎞ ⎛ (4.1.1) = p ⎜1 − ⎟ x dt ⎝ K⎠ dengan p merupakan laju pertumbuhan intrinsik yang berhubungan dengan kelahiran dan kematian alami populasi. Parameter K
menyatakan daya muat lingkungan (carrying capacity) yang menyatakan kapasitas maksimum populasi x dalam lingkungan tersebut. Penambahan setiap satu individu pemangsa I ke dalam populasi akan mengakibatkan bertambahnya laju pemangsaan pemangsa I, sehingga laju pertumbuhan per kapita mangsa akan menurun, dengan a adalah besarnya laju pemangsaan dari pemangsa I dan y adalah banyaknya individu pemangsa I, maka persamaan (4.1.1) dimodifikasi menjadi: 1 dx x⎞ ⎛ = p ⎜1 − ⎟ − ay (4.1.2) x dt K ⎝ ⎠ yang sesuai dengan persamaan (3.2.2.1). Selanjutnya dengan cara yang sama, penambahan setiap satu individu pemangsa II ke dalam populasi mengakibatkan bertambahnya laju pemangsaan pemangsa II yang akan menyebabkan menurunnya laju pertumbuhan per kapita dari mangsa sehingga persamaan (4.1.2) dimodifikasi menjadi:
12
1 dx x⎞ ⎛ (4.1.3) = p ⎜1 − ⎟ − ay − bz x dt K ⎝ ⎠ dengan b adalah besarnya laju pemangsaan dari pemangsa II dan z adalah banyaknya individu pemangsa II.
=
Dari persamaan (3.2.1.1), ry ( K y − y − α yz z ) Ky ry ( K y − y ) Ky
−
ryα yz z Ky
.
Dengan menyubstitusikan α yz =
(4.2.1) a yz
dari
a yy
persamaan (3.2.1.3) ke persamaan (4.2.1) diperoleh ⎛ Ky dy y ⎞ ry a yz = ry ⎜ − z. ⎟− ⎜ K y K y ⎟ K y a yy y dt ⎝ ⎠
maka
persamaan
(4.2.5) menjadi: dy y⎞ ⎛ = q ⎜ 1 − ⎟ + cx − dz − vU . y dt L ⎝ ⎠
(4.2.6)
4.3 Laju pertumbuhan pemangsa II
4.2 Laju pertumbuhan pemangsa I
dy = y dt
q = ry , L = K y , d = a yz ,
(4.2.2)
Kemudian dengan menyubstitusi K y =
ry a yy
dari persamaan (3.1.2) ke dalam persamaan (4.2.2) diperoleh ⎛ dy y ⎞ = ry ⎜ 1 − (4.2.3) ⎟−a z ⎜ K y ⎟ yz y dt ⎝ ⎠ Penambahan setiap satu individu mangsa ke dalam populasi mengakibatkan bertambahnya laju pemangsaan dari pemangsa I sehingga energi yang terkandung pada mangsa akan dikonversi menjadi energi pemangsa I, dan hal ini akan mengakibatkan bertambahnya pertumbuhan per kapita pemangsa I, sehingga persamaannya dimodifikasi menjadi: ⎛ dy y ⎞ = ry ⎜ 1 − (4.2.4) ⎟ − a yz z + cx ⎜ ⎟ y dt ⎝ Ky ⎠ dengan c adalah faktor konversi yang menetapkan banyaknya pemangsa I baru yang dilahirkan untuk tiap mangsa yang ditangkap. Munculnya para nelayan yang memanen pemangsa I akan mengurangi pertumbuhan per kapita pemangsa I, sehingga persamaannya dimodifikasi menjadi ⎛ dy y ⎞ (4.2.5) = ry ⎜ 1 − ⎟ − a yz z + cx − vU ⎜ ⎟ y dt ⎝ Ky ⎠ dengan v adalah koefisien ketertangkapan pemangsa I dan U adalah usaha penangkapan yang dilakukan oleh para nelayan. Misalkan,
Dari persamaan (3.2.1.2), rz ( K z − z − α zy y ) dz = z dt Kz r ( K − z ) rzα zy y = z z − Kz Kz Dengan
menyubstitusi
α zy =
(4.3.1) azy
dari azz ke persamaan (4.3.1)
persamaan (3.2.1.4) diperoleh ⎛K dz z ⎞ rz azy y. = rz ⎜ z − ⎟− z dt ⎝ K z K z ⎠ K z a zz
(4.3.2)
rz azz dari persamaan (3.1.2) ke dalam persamaan (4.3.2) diperoleh ⎛ dz z ⎞ (4.3.3) = rz ⎜1 − ⎟ − azy y. z dt K z ⎠ ⎝ Penambahan setiap satu individu mangsa ke dalam populasi akan mengakibatkan bertambahnya laju pemangsaan dari pemangsa II sehingga energi yang terkandung pada mangsa akan dikonversi menjadi energi pemangsa II, dan mengakibatkan bertambahnya pertumbuhan per kapita pemangsa II, sehingga persamaannya menjadi: ⎛ dz z ⎞ = rz ⎜ 1 − (4.3.4) ⎟ − azy y + ex z dt ⎝ Kz ⎠ dengan e adalah faktor konversi yang menyatakan banyaknya pemangsa II baru yang dilahirkan untuk tiap mangsa yang ditangkap. Munculnya para nelayan yang memanen pemangsa II akan mengurangi pertumbuhan per kapita pemangsa II, sehingga persamaannya dimodifikasi menjadi ⎛ dz z ⎞ (4.3.5) = rz ⎜ 1 − ⎟ − azy y + ex − wU z dt ⎝ Kz ⎠
Kemudian dengan menyubstitusi K z =
dengan w adalah koefisien ketertangkapan pemangsa II dan U adalah usaha penangkapan yang dilakukan oleh para nelayan. Misalkan, r = rz , M = K z , f = a zy , maka persamaan (4.3.5) dimodifikasi menjadi:
13
dz z ⎛ = r ⎜1 − z dt M ⎝
⎞ ⎟ + ex − fy − wU ⎠
(4.3.6)
U
U
v
c
w d
a
y
x
z
b
f y⎞ ⎛ q ⎜1 − ⎟ ⎝ L⎠
e
x⎞ ⎛ p ⎜1 − ⎟ ⎝ K⎠
z ⎞ ⎛ r ⎜1 − ⎟ ⎝ M⎠
Gambar 10 Diagram kedinamikan interaksi antarspesies.
Persamaan (4.1.3), (4.2.6), (4.3.6) dapat disusun menjadi suatu sistem persamaan diferensial taklinear
⎫ ⎪ ⎪ ⎪⎪ ⎬ (4.3.7) ⎪ ⎪ dz z ⎞ ⎛ = rz ⎜ 1 − ⎟ + ezx − fzy − wUz ⎪ dt ⎝ M⎠ ⎭⎪ dengan dx = laju pertumbuhan individu mangsa, dt dy = laju pertumbuhan individu pemangsa I, dt dz = laju pertumbuhan individu pemangs II, dt K , L, M > 0 , p, q, r > 0 , 0 ≤ U ≤ U maks , dx x⎞ ⎛ = px ⎜1 − ⎟ − axy − bxz dt K ⎝ ⎠ dy y ⎛ ⎞ = qy ⎜ 1 − ⎟ + cxy − dyz − vUy dt ⎝ L⎠
0 < a < 1 , 0 < b < 1, 0 < c < 1 , 0 < d < 1, 0 < e < 1 , 0 < f < 1, 0 < v < 1 , 0 < w < 1 , dan x, y, z ≥ 0 .
V ANALISIS KESTABILAN
Pada bab ini akan dianalisis kestabilan dari tiap titik tetap pada sistem (4.3.7). Sebagai langkah awal akan ditentukan titik tetapnya, kemudian sistem ini akan dilinearkan dengan mengonstruksi matriks Jacobinya dan selanjutnya menganalisis kestabilannya dengan memeriksa nilai eigen.
Dari (5.1.1) diperoleh
5.1 Penentuan titik tetap :
⎫ ⎪⎪ ⎛ ⎛ ⎞ ⎬ y∗ ⎞ ∗ ∗ ⎜ q ⎜1 − ⎟ + cx − dz − vU ⎟ = 0 ⎪ ⎜ ⎟ ⎜ ⎟ L ⎠ ⎪⎭ ⎝ ⎝ ⎠
Titik
tetap
(x , y , z ) ∗
∗
∗
dengan membuat: dx dy dz = 0, = 0, =0. dt dt dt Dari persamaan (4.3.7) maka ⎛ ⎛ x∗ ⎞ ⎞ x∗ ⎜ p ⎜1 − ⎟ − ay∗ − bz∗ ⎟ = 0 ⎜ ⎜ ⎟ K ⎟⎠ ⎝ ⎝ ⎠
(5.1.4)
Dari (5.1.2) didapat : y∗ = 0 atau
diperoleh
(5.1.5)
Dari (5.1.3) didapat : ⎫ ⎪⎪ ⎛ ⎛ ⎞ ⎬ z∗ ⎞ ⎜ r ⎜1 − ⎟ + ex∗ − fy∗ − wU ⎟ = 0 ⎪ ⎜ ⎜ M⎟ ⎟ ⎠ ⎝ ⎝ ⎠ ⎭⎪ z ∗ = 0 atau
(5.1.1)
⎛ ⎛ ⎞ y∗ ⎞ y∗ ⎜ q ⎜ 1 − ⎟ + cx∗ − dz ∗ − vU ⎟ = 0 ⎜ ⎜ ⎟ L ⎟⎠ ⎝ ⎝ ⎠
(5.1.2)
⎛ ⎛ z∗ z∗ ⎜ r ⎜1 − ⎜ ⎜ M ⎝ ⎝
(5.1.3)
⎞ ⎞ ∗ ∗ ⎟⎟ + ex − fy − wU ⎟ = 0 ⎟ ⎠ ⎠
⎫ ⎪⎪ ⎛ ⎛ x∗ ⎞ ⎞ ⎬ ⎜ p ⎜1 − ⎟ − ay∗ − bz ∗ ⎟ = 0 ⎪ ⎜ ⎜ ⎟ K ⎟⎠ ⎝ ⎝ ⎠ ⎭⎪ x∗ = 0 atau
(5.1.6)
Dari persamaan (5.1.4) - (5.1.6) diperoleh titik tetap sebagai berikut :
14
E0 = (0, 0, 0) E1 = ( x1 , 0, 0)
E2 = ( x2 , y2 , 0 ) E3 = ( x3 , 0, z3 ) E4 = ( x4 , y4 , z4 ) E5 = ( 0, y5 , 0 ) E6 = ( 0, y6 , z6 ) E7 = ( 0, 0, z7 )
x4 =
z6 =
(5.1.7)
y2 =
x3 = z3 =
K ( pq + aL ( − q + Uv ) ) acKL + pq Lp ( cK + q − Uv ) acKL + pq
K ( pr + bM ( − r + Uw ) ) beKM + pr Mp ( eK + r − Uw ) beKM + pr
(
)
− LM (adeK + bcfK + dfp ) + beKMq + (acKL + pq )r
z4 =
y6 =
x2 =
K r ( − aLq + pq + aLUv ) − dLM ( fp − ar + aUw ) + bM ( − qr + fL ( q − Uv ) + qUw )
y4 =
y5 = L −
dengan x1 = K
⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎭
L ( pr ( cK + q − Uv ) − dMp ( eK + r − Uw ) + bKM ( eq − cr − eUv + cUw ) ) − LM (adeK + bcfK + dfp ) + beKMq + (acKL + pq)r
M ( aeKL ( q − Uv ) + cKL ( fp − ar + aUw ) − p ( eKq − fLq + qr + fLUv − qUw ) ) M ( adeKL + dfLp + bK ( cfL − eq ) ) − ( acKL + pq ) r
5.2 Konstruksi Matriks Jacobi
LvU q
L ( −qr + rUv + dM ( r − Uw ) ) dfLM − qr M ( fL(q − Uv) + q ( −r + Uw ) ) dfLM − qr
MUw r (lihat Lampiran 1). Diasumsikan titik tetap E1 sampai E7 ada, z7 = M −
yaitu jika nilai-nilai x1 , x2 , x3 , x4 , y2 , y4 , y5 , y6 , z3 , z4 , z6 , z7 positif. Penentuan titik keseimbangan ini diperoleh dengan menggunakan Software Mathematica 5.1.
Misalkan SPD taklinear pada persamaan (2.2.1) dimodifikasi untuk n = 3 sebagai berikut dx ⎫ = P ( x, y , z ) ⎪ dt ⎪ dy ⎪ (5.2.1) = Q ( x, y , z ) ⎬ dt ⎪ ⎪ dz = R ( x, y , z ) ⎪ dt ⎭ dengan turunan-turunan parsial fungsi P, Q dan R kontinu di R3 . Dengan menggunakan perluasan Taylor (Teorema 1) di sekitar titik
(
)
tetap x∗ , y ∗ , z ∗ , diperoleh
15
P ( x, y , z ) = P ( x ∗ , y ∗ , z ∗ ) +
∂P ( x∗ , y ∗ , z ∗ ) ∂P ( x∗ , y ∗ , z ∗ ) ∂P ( x∗ , y ∗ , z ∗ ) x − x∗ + y − y∗ + z − z ∗ + ϕ1 ( x, y, z ) ∂x ∂y ∂z
Q ( x, y , z ) = Q ( x ∗ , y ∗ , z ∗ ) +
∂Q ( x∗ , y ∗ , z ∗ ) ∂Q ( x∗ , y ∗ , z ∗ ) ∂Q ( x∗ , y ∗ , z ∗ ) x − x∗ + y − y∗ + z − z ∗ + ϕ 2 ( x, y , z ) ∂x ∂y ∂z
R ( x, y , z ) = R ( x ∗ , y ∗ , z ∗ ) +
∂R ( x∗ , y ∗ , z ∗ ) ∂R ( x∗ , y ∗ , z ∗ ) ∂R ( x∗ , y ∗ , z ∗ ) x − x∗ + y − y∗ + z − z ∗ + ϕ 3 ( x, y , z ) ∂x ∂y ∂z
(
)
(
(
)
(
(
)
(
)
)
)
(
)
(
)
(
)
(5.2.2) ⎡ϕ1 ( x, y, z ) ⎤ ⎢ ⎥ ϕ ( Z ) = ⎢ϕ2 ( x, y, z ) ⎥ ⎢ ⎥ ⎣⎢ϕ3 ( x, y, z ) ⎦⎥
dengan lim ∗
( v ) →( v
ϕ i ( x, y , z )
)
(
x − x∗
) ( 2
+ y − y∗
) ( 2
+ z − z∗
)
2
=0
Karena suku ϕ ( Z ) pada sistem (5.2.3) sangat
i = 1, 2,3 dan
(
)
v = ( x, y , z ) , v ∗ = x ∗ , y ∗ , z ∗ .
Jika didefinisikan variabel baru ξ1 = x − x∗
ξ2 = y − y∗
ξ3 = z − z ∗ sehingga
(
)
(
)
(
)
∗ d ξ1 d x − x = = P ( x, y , z ) − P ( x ∗ , y ∗ , z ∗ ) dt dt ∗ dξ2 d y − y = = Q ( x, y , z ) − Q ( x ∗ , y ∗ , z ∗ ) dt dt ∗ d ξ3 d z − z = = R ( x, y , z ) − R ( x ∗ , y ∗ , z ∗ ) dt dt
maka SPD (5.2.2) dapat direduksi menjadi dZ = AZ + ϕ ( Z ) (5.2.3) dt dengan ⎡ Px ⎢ A = ⎢Qx ⎢ Rx ⎣
Py Qy Ry
Pz ⎤ ⎥ Qz ⎥ Rz ⎥⎦
⎡ x − x∗ ⎤ ⎢ ⎥ Z = ⎢ y − y∗ ⎥ ⎢ z − z∗ ⎥ ⎣ ⎦
kecil dibandingkan suku linear AZ di sekitar titik tetap, maka ϕ ( Z ) pada sistem (5.2.3) dapat diabaikan sehingga diperoleh dZ = AZ . (5.2.4) dt Persamaan (5.2.4) dapat dianggap sebagai hampiran yang memenuhi sistem taklinier (5.2.3). Matriks koefisien A untuk persamaan (5.2.3) merupakan matriks Jacobi yang dievaluasi pada titik tetap
(x , y , z ), ∗
∗
∗
dan
Px , Py , Pz , Qx , Qy , Qz , Rx , Ry , dan Rz adalah turunan parsial fungsi yang bersesuaian dengan peubah yang dimaksud. Bentuk persamaan (5.2.4) disebut pelinearan dari sistem persamaan diferensial (5.2.3). Dengan melakukan pelinearan pada sistem persamaan di sekitar titik tetapnya, maka diperoleh matriks Jacobi dari sistem (4.3.7) sebagai berikut: ⎡ ∂P ∂P ∂P ⎤ ⎢ ⎥ ⎢ ∂x ∂y ∂z ⎥ ⎢ ∂Q ∂Q ∂Q ⎥ (5.2.5) A x∗ , y ∗ , z ∗ = ⎢ ⎥ ⎢ ∂x ∂y ∂z ⎥ ⎢ ∂R ∂R ∂R ⎥ ⎢ ⎥ ⎣ ∂x ∂y ∂z ⎦ atau ⎡ A11 A12 A13 ⎤ ∗ ∗ ∗ (5.2.6) A x , y , z = ⎢⎢ A21 A22 A23 ⎥⎥ ⎢⎣ A31 A32 A33 ⎥⎦
(
)
(
)
dengan A11 = p −
2 px − ay − bz K
16
A12 = −ax
5.3.1 Kestabilan titik tetap E0 Dengan menggunakan Software Mathematica 5.1 (lihat Lampiran 2) diperoleh nilai eigen sebagai berikut: λ1 = p, λ2 = q − Uv, λ3 = r − Uw Karena dari asumsi p > 0 , maka λ1 > 0 sehingga titik tetap E0 bersifat tidak stabil.
A13 = −bx A21 = cy A22 = q − Uv + cx − A23 = − dy
2qy − dz L
A31 = ez A32 = − fz 2rz A33 = r − Uw + ex − fy − M Matriks Jacobi di titik tetapnya diperoleh dengan menyubstitusi E0 , E1 , E2 , E3 , E4 , E5 , E6 , dan E7 pada persamaan (5.1.7) ke dalam matriks A pada persamaan (5.2.6) sehingga dapat disusun persamaan karakteristik dan diperoleh nilai-nilai eigen pada masing-masing titik tetapnya (lihat Lampiran 2).
5.3 Analisis kestabilan titik tetap
Kondisi kestabilan dari titik tetap E0 , E1 , E2 , E3 , E5 , E6 , dan E7 diperoleh dengan memeriksa nilai-nilai eigen dari matriks A , kemudian disesuaikan dengan kondisi kestabilan pada sub bab 2.2. Khusus untuk titik tetap E4 , akan ditunjukkan bersifat stabil asimtotik global.
5.3.2 Kestabilan titik tetap E1 Dengan Software Mathematica 5.1 (lihat Lampiran 2) diperoleh nilai eigen sebagai berikut: λ1 = − p, λ2 = cK + q − Uv, λ3 = eK + r − Uw . Karena λ1 < 0, maka titik tetap E1 bersifat stabil jika λ2 < 0 , λ3 < 0 .
Dari
λ2 = cK + q − Uv < 0
diperoleh
cK + q , karena λ3 = eK + r − Uw < 0 , v eK + r , Sehingga titik maka diperoleh U > w tetap E1 bersifat stabil jika U >
⎧ cK + q eK + r ⎫ U > max ⎨ , ⎬. w ⎭ ⎩ v Berikut ini nilai-nilai eigen untuk titik tetap E2 , E3 , E5 , E6 , dan E7 yang ditentukan dari Software Mathematica 5.1 (lihat Lampiran 2).
Nilai eigen di titik tetap E2 λ1 =
1 (−cKp + aLp ( q − Uv ) − pq ( p + q − Uv ) 2 ( acKL + pq ) +
λ2 = − +
( p ( p ( q (cK − aL + p + q ) + ( aL − q )Uv ) − 4 ( acKL + pq )( cK + q − Uv ) ( pq + aL ( −q + Uv )))) ) 2
1 (cKp − aLp ( q − Uv ) + pq ( p + q − Uv ) 2 ( acKL + pq )
( p ( p ( q (cK − aL + p + q ) + ( aL − q )Uv ) − 4 ( acKL + pq )( cK + q − Uv ) ( pq + aL ( −q + Uv )))) )
λ3 = r +
2
−cfKLp + aeKL ( − q + Uv ) + p ( eKq − fLq + fLUv ) acKL + pq
− Uw
17
Nilai eigen di titik tetap E3 Out[367]=
:l1= Bq +
l2= B
p r Hc K - U v L - d M p He K + r - U wL - b K M Hc r + e U v - c U wL beK M+ pr
2 Hb e K M + p rL 1
H- e K p r + b M p Hr - U wL -
p r H p + r - U wL +
l3= B-
F,
,
H p H p Hr He K - b M + p + rL + Hb M - rL U wL2 -
4 Hb e K M + p rL He K + r - U wL H p r + b M H- r + U wLLLLLF ,
He K p r - b M p Hr - U wL + p r H p + r - U wL + 2 Hb e K M + p rL , H p H p Hr He K - b M + p + rL + Hb M - rL U wL2 -
1
4 Hb e K M + p rL He K + r - U wL H p r + b M H- r + U wLLLLLF>
Uvy Uvy i zzF , l3 = Br + f L ijj - 1 + zz - U wF> :l1= @- q + U v D , l2 = B p + a L jj - 1 + q { q { k k
Nilai eigen di titik tetap E5 Out[371]=
Nilai eigen di titik tetap E6 1 Out[373]= :l1= B Hr H- p q + a L Hq - U v LL + d f L M-qr
l2= B
d L M H f p - a r + a U wL + b M H- f L q + q r + f L U v - q U wLLF , 1
H- f L r Hq - U v L + d M q H- r + U wL +
, q r Hq + r - U Hv + wLL + H4 Hd f L M - q rL H-q r + r U v + d M Hr - U wLL H f L Hq - U v L + q H- r + U wLL +
2d f L M-2qr
H f L r Hq - U v L + q Hd M Hr - U wL - r Hq + r - U Hv + wLLLL LLF , 2
l3= B
1
H- f L r Hq - U v L + d M q H- r + U wL +
, q r Hq + r - U Hv + wLL - H4 Hd f L M - q rL H-q r + r U v + d M Hr - U wLL H f L Hq - U v L + q H- r + U wLL +
2d f L M-2qr
H f L r Hq - U v L + q Hd M Hr - U wL - r Hq + r - U Hv + wLLLL LLF> 2
Nilai eigen di titik tetap E7 Out[375]=
:l1= @- r + U wD , l2 = B p + b M J- 1 +
5.3.3 Kestabilan titik tetap E2 Sesuai dengan asumsi, titik tetap E2 ada jika x2 > 0 dan y2 > 0 . x2 =
K ( pq + aL ( − q + Uv ) )
⇔U >
acKL + pq
( aL − p ) q v
.
>0
Uw r
NF , l3 = Bq - U v + d M J- 1 +
y2 =
Lp ( cK + q − Uv ) acKL + pq
Uw r
NF>
>0
cK + q . v Jadi titik tetap E2 ada jika ⇔U <
( aL − p ) q v
cK + q v
(5.3.3.1)
18
Titik tetap E2 bersifat stabil jika Re ( λ1 ) < 0 , Re ( λ2 ) < 0 , dan Re ( λ3 ) < 0 . Re ( λ1 ) = −
1 cKpq − aLp ( q − Uv ) + pq ( p + q − Uv ) < 0 2 ( acKL + pq )
(
)
⇔ ( cKpq − aLp ( q − Uv ) + pq ( p + q − Uv ) ) > 0 ⇔ cKpq − aLpq + p 2 q + pq 2 + aLpUv − pqUv > 0
jika cKpq + pq 2 − pqUv > 0 dan −aLpq + p 2 q + aLpUv > 0 . cK + q v aL − p ) q ( −aLpq + p 2 q + aLpUv > 0 ⇔ U > v ( aL − p ) q cK + q sehingga Re ( λ1 ) < 0 jika
0 ⇔ U <
Re ( λ2 ) =
(5.3.3.2)
1 −cKpq + aLp ( q − Uv ) − pq ( p + q − Uv ) < 0 2 ( acKL + pq )
(
)
( ) ⇔ ( cKpq − aLp ( q − Uv ) + pq ( p + q − Uv ) ) > 0 ⇔ −cKpq + aLp ( q − Uv ) − pq ( p + q − Uv ) < 0
⇔ cKpq − aLpq + p 2 q + pq 2 + aLpUv − pqUv > 0. Jadi Re ( λ2 ) = Re ( λ1 ) . Re ( λ3 ) = r +
−cfKLp + aeKL ( − q + Uv ) + p ( eKq − fLq + fLUv )
acKL + pq Jika kedua ruas dikalikan dengan acKL + pq diperoleh
− Uw < 0
( acKL + pq ) r − cfKLp + aeKL ( −q + Uv ) + p ( eKq − fLq + fLUv ) − Uw ( acKL + pq ) < 0 ⇔ −cfKLp − aeKLq + eKpq − fLpq + acKLr + pqr + U ( aeKLv + fLpv − acKLw − pqw ) < 0 ⇔ −cfKLp − aeKLq + eKpq − fLpq + acKLr + pqr < U ( − aeKLv − fLpv + acKLw + pqw ) . Jika −aeKLv − fLpv + acKLw + pqw > 0, maka U>
−cfKLp − aeKLq + eKpq − fLpq + acKLr + pqr ( −aeKLv − fLpv + acKLw + pqw )
⇔U >
eKpq + ( acKL + pq ) r − cfpKL − aeqKL − fpqL
(5.3.3.3)
w ( acKL + pq ) − v ( aeKL + fpL )
Jadi titik tetap E 2 bersifat stabil jika U>
eKpq + ( acKL + pq ) r − cfpKL − aeqKL − fpqL
(5.3.3.4)
w ( acKL + pq ) − v ( aeKL + fpL )
dengan w ( acKL + pq ) − v ( aeKL + fpL ) > 0 5.3.4 Kestabilan titik tetap E3 Titik tetap E3 ada jika x3 > 0 dan z3 > 0 .
⇔U >
x3 =
z3 =
K ( pr + bM ( − r + Uw ) ) beKM + pr
⇔ pr + bM ( − r + Uw ) > 0
>0
( bM − p ) r
bMw Mp ( eK + r − Uw ) beKM + pr
⇔ ( eK + r − Uw ) > 0
>0
19
( bM − p ) r
eK + r w Jadi titik tetap E3 ada jika ⇔U <
bMw
eK + r . w
(5.3.4.1)
Titik tetap E3 bersifat stabil jika Re ( λ1 ) < 0 , Re ( λ2 ) < 0 , dan Re ( λ3 ) < 0 . Re ( λ1 ) = q +
pr ( cK − Uv ) − dMp ( eK + r − Uw ) − bKM ( cr + eUv − cUw ) beKM + pr
<0.
Jika kedua ruas dikalikan dengan beKM + pr diperoleh q ( beKM + pr ) + pr ( cK − Uv ) − dMp ( eK + r − Uw) − bKM ( cr + eUv − cUw ) < 0
⇔ −deKMp + beKMq − bcKMr + cKpr − dMpr + pqr + U (−beKMv − prv + bcKMw + dMpw) < 0 ⇔ −deKMp + beKMq − bcKMr + cKpr − dMpr + pqr < U (beKMv + prv − bcKMw − dMpw) .
Jika beKMv + prv − bcKMw − dMpw > 0, maka U>
− deKMp + beKMq − bcKMr + cKpr − dMpr + pqr (beKMv + prv − bcKMw − dMpw)
⇔U >
q ( beKM + pr ) + cKpr − M ( bcrK + depK + dpr ) v ( beKM + pr ) − Mw ( bcK − dp )
Re ( λ2 ) = −
(5.3.4.2)
1 ( eKpr − bMp ( r − Uw) + pr ( p + r − Uw) ) < 0 2 ( beKM + pr )
⇔ ( eKpr − bMp ( r − Uw ) + pr ( p + r − Uw ) ) > 0 ⇔ eKpr − bMpr + p 2 r + pr 2 + bMpUw − prUw > 0
jika
eKpr + pr 2 − prUw > 0 dan −bMpr + p 2 r + bMpUw > 0
eK + r w ( bM − p ) r −bMpr + p 2 r + bMpUw > 0 ⇔ U > bMw ( bM − p ) r eK + r sehingga Re ( λ2 ) < 0 jika 0 ⇔ U <
Re ( λ3 ) =
(5.3.4.3)
1 ( −eKpr + bMp ( r − Uw) − pr ( p + r − Uw ) ) < 0 2 ( beKM + pr )
⇔ ( −eKpr + bMp ( r − Uw ) − pr ( p + r − Uw ) ) < 0
⇔ ( eKpr − bMp ( r − Uw ) + pr ( p + r − Uw ) ) > 0 ⇔ eKpr − bMpr + p 2 r + pr 2 + bMpUw − prUw > 0 Perhatikan bahwa Re ( λ3 ) = Re ( λ2 ) .
Jadi titik tetap E3 bersifat stabil jika U>
q ( beKM + pr ) + cKpr − M ( bcrK + depK + dpr )
dengan
v ( beKM + pr ) − Mw ( bcK − dp )
(5.3.4.4)
v ( beKM + pr ) − Mw ( bcK − dp ) > 0 .
Dengan menggunakan cara yang sama, batasbatas kestabilan di titik tetap E5 , E6 , dan E7 dapat ditentukan.
5.3.5 Kestabilan titik tetap E4 Titik tetap E4 ada jika x4 > 0, y4 > 0, dan z4 > 0 .
20
Teorema 8 Titik tetap interior E4 adalah stabil global asimtotik jika i. a = c, b = e; 4rq 2 > (d + f ) LM (Chattopadhyay et al., 1996)
ii.
Bukti. Pilih fungsi Liapunov V ( x, y, z ) = V1 ( x ) + V2 ( y ) + V3 ( z )
x
∗
) = p − ln (1 + p ) . (
V1 x∗ + ε
bahwa
x
∗
) >0 ⇔
p > ln (1 + p )
⇔ e p > (1 + p ) ⇔ e p − p −1 > 0 .
Bukti. Misalkan
f ( p) = e p − p −1
⎛ y ⎞ + ( y − y ) − y ln ⎜⎜ ∗ ⎟⎟ ⎝y ⎠
merupakan fungsi naik untuk p > 0 , sehingga
⎛ z ⎞ + ( z − z ) − z ln ⎜ ∗ ⎟ ⎝z ⎠
e − p − 1 > e − 0 − 1 = 0 , jika
p>0,
f ′( p) > 0 .
∗
maka
Dari
f ( p ) > f (0)
∗
V ( x, y, z ) adalah nol. Untuk ( x, y, z ) > 0 akan ditunjukkan i. Jika x > x* , y > y* , z > z* maka
ep >1,
Teorema
jika
dan
0
p
dengan x, y, z > 0 dan x∗ , y∗ , z ∗ adalah konstanta positif. Untuk x = x* , y = y* , z = z* maka
berarti
(
V1 x∗ + ε
) > 0 . Akibatnya V ( x ) > 0 .
x Dengan cara yang sama dapat ditunjukkan bahwa V2 ( y ) > 0 dan V3 ( z ) > 0 untuk
y > y* , z > z * .
x < x* . Misalkan terdapat
Untuk
ii.
(
)
x = x∗ − ε .
x = x∗ + ε . ⎛ x ⎞ V1 x∗ + ε = ( x − x∗ ) − x∗ ln ⎜ ∗ ⎟ ⎝x ⎠
(
)
⎛ x∗ + ε = ε − x ln ⎜ ∗ ⎜ x ⎝ ∗
⎞ ⎟⎟ ⎠
ε ⎞ ⎛ = ε − x∗ ln ⎜ 1 + ∗ ⎟ . x ⎠ ⎝ Jika kedua ruas dibagi dengan x∗ , diperoleh
(
x∗
Misalkan
)
= p=
ε x∗
ε x
∗
ε ⎞ ⎛ − ln ⎜ 1 + ∗ ⎟ . x ⎝ ⎠ , maka
p>0
sehingga ⎛ x∗ − ε − x∗ ln ⎜ ∗ ⎜ x ⎝
Dengan kata lain, akan ditunjukkan V ( x, y, z ) adalah fungsi definit positif untuk ( x, y , z ) ≠ ( x ∗ , y ∗ , z ∗ ) sehingga
Ini
1
∗
V1 x∗ − ε = ( x∗ − ε − x∗ )
ε >0
jika
p>0.
Jika x < x* , y < y* , z < z* maka V ( x, y, z ) adalah positif .
x > x* . Misalkan terdapat
f ( p)
hanya
ε >0
Untuk
sehingga
6,
V ( x, y, z ) adalah positif .
V1 x∗ + ε
akan
f ′( p) = e p −1.
f ( p) > 0 .
ditunjukkan
maka
Karena
∗
i.
Akan dibuktikan
⎛ x ⎞ = ( x − x∗ ) − x∗ ln ⎜ ∗ ⎟ ⎝x ⎠ ∗
ii.
(
V1 x∗ + ε
⎞ ⎟⎟ ⎠
ε ⎞ ⎛ = −ε − x∗ ln ⎜1 − ∗ ⎟ . ⎝ x ⎠ Jika kedua ruas dibagi dengan x∗ diperoleh
(
V1 x∗ − ε x p=
∗
ε
) =− ε
ε ⎞ ⎛ − ln ⎜ 1 − ∗ ⎟ . x ⎝ x ⎠
Misalkan
∗
(
)
= x x − p − ln (1 − p ) . Fungsi ln (1 − p ) terdefinisi ∗
, maka p > 0 sehingga
V1 x∗ − ε ∗
jika 1 − p > 0 ⇔ p < 1 . Karena p > 0 , maka 0 < p < 1 . Akan dibuktikan
(
V1 x∗ − ε x
∗
) > 0 ⇔ − p > ln (1 − p ) ⇔ e− p > 1 − p .
sehingga
Bukti. Misalkan ditunjukkan
f ( p ) = e− p + p − 1 f ( p) > 0 .
maka
akan
f ′ ( p ) = −e − p + 1 .
21
e− p < 1 ,
Nilai
p>0,
untuk
akibatnya
f ′ ( p ) > 0 untuk 0 < p < 1 . Dari Teorema 6, f ( p)
dapat ditunjukkan bahwa V2 ( y ) > 0
dan
V3 ( z ) > 0 untuk y < y , z < z . *
*
untuk
Dari (i) dan (ii) terbukti bahwa V ( x, y, z ) > 0
0 < p < 1 , sehingga f ( p ) > f ( 0 ) jika dan
untuk x, y, z > 0 dengan x∗ , y∗ , z∗ adalah konstanta positif, sehingga V ( x, y , z ) adalah fungsi definit positif. Turunan pertama dari fungsi V terhadap waktu dapat dicari dengan cara sebagai berikut
merupakan
hanya
fungsi
naik
e p + p − 1 > e0 + 0 − 1 = 0
jika
ep > 1− p .
Ini
berarti
(
∗
V1 x − ε ∗
⇔
) >0.
x Akibatnya V1 ( x ) > 0 . Dengan cara yang sama
⎛ y ⎞ ⎛ x ⎞ ⎛ z ⎞ V ( x, y, z ) = ( x − x∗ ) − x∗ ln ⎜ ∗ ⎟ + ( y − y∗ ) − y∗ ln ⎜⎜ ∗ ⎟⎟ + ( z − z ∗ ) − z ∗ ln ⎜ ∗ ⎟ ⎝z ⎠ ⎝x ⎠ ⎝y ⎠ dV ∂V dx ∂V dy ∂V dz = ⋅ + ⋅ + ⋅ V′ = ∂x dt ∂y dt ∂z dt dt ⎛ x∗ ⎞ dx ⎛ y ∗ ⎞ dy ⎛ z ∗ ⎞ dz = ⎜⎜ 1 − ⎟⎟ + ⎜⎜ 1 − ⎟⎟ + ⎜ 1 − ⎟⎟ y ⎠ dt ⎜⎝ z ⎠ dt x ⎠ dt ⎝ ⎝ ∗ ⎛ x − x ⎞⎛ ⎛ ⎞ ⎛ y − y∗ ⎞ ⎛ ⎛ ⎞ x⎞ y⎞ = ⎜⎜ ⎟⎟ ⎜ px ⎜ 1 − ⎟ − axy − bxz ⎟ + ⎜⎜ ⎟⎟ ⎜ qy ⎜ 1 − ⎟ + cxy − dyz − vUy ⎟ x K y L ⎝ ⎠ ⎝ ⎠ ⎠ ⎝ ⎠ ⎝ ⎠⎝ ⎠⎝ ⎛ z − z∗ ⎞ ⎛ ⎛ ⎞ z ⎞ + ⎜⎜ ⎟⎟ ⎜ rz ⎜1 − ⎟ + exz − fyz − wUz ⎟ z M ⎠ ⎠ ⎝ ⎠⎝ ⎝ x⎞ y⎞ ⎛ ⎛ ⎞ ⎛ ⎛ ⎞ = x − x∗ ⎜ p ⎜1 − ⎟ − ay − bz ⎟ + y − y∗ ⎜ q ⎜1 − ⎟ + cx − dz − vU ⎟ K L ⎠ ⎠ ⎝ ⎝ ⎠ ⎝ ⎝ ⎠
(
(
)
⎛ ⎛ z + z − z∗ ⎜ r ⎜1 − M ⎝ ⎝
(
)
)
⎞ ⎞ ⎟ + ex − fy − wU ⎟ ⎠ ⎠
Di titik keseimbangannya ( x∗ , y∗ , z ∗ ) : berlaku ⎛ ⎛ x∗ ⎞ ⎞ x∗ ⎜ p ⎜1 − ⎟ − ay ∗ − bz ∗ ⎟ = 0 ⎟ ⎜ ⎜ ⎟ K⎠ ⎝ ⎝ ⎠ ⎛ ⎛ ⎞ y∗ ⎞ y ∗ ⎜ q ⎜ 1 − ⎟ + cx∗ − dz ∗ − vU ⎟ = 0 ⎜ ⎟ ⎜ ⎟ L ⎠ ⎝ ⎝ ⎠
⎛ ⎛ z∗ z∗ ⎜ r ⎜1 − ⎜ ⎜ M ⎝ ⎝
⎞ ⎞ ∗ ∗ ⎟⎟ + ex − fy − wU ⎟ = 0 ⎟ ⎠ ⎠
Di titik E4 , ( x∗ , y ∗ , z ∗ ) > ( 0, 0, 0 ) sehingga ⎛ ⎛ x∗ ⎞ ⎞ ⎜ p ⎜ 1 − ⎟ − ay ∗ − bz ∗ ⎟ = 0 ⎜ ⎟ ⎜ ⎟ K⎠ ⎝ ⎝ ⎠ ⎛ ⎛ ⎞ y∗ ⎞ ⎜ q ⎜1 − + cx∗ − dz ∗ − vU ⎟ = 0 ⎟ ⎜ ⎜ ⎟ L ⎟⎠ ⎝ ⎝ ⎠
⎛ ⎛ z∗ ⎜ r ⎜1 − ⎜ ⎜ M ⎝ ⎝
⎞ ⎞ ∗ ∗ ⎟⎟ + ex − fy − wU ⎟ = 0 ⎟ ⎠ ⎠
22
V ′ dapat dituliskan sebagai ⎛ ⎛ ⎛ ⎛ x∗ ⎞ ⎞⎞ x⎞ V ′ = x − x∗ ⎜ p ⎜1 − ⎟ − ay − bz − ⎜ p ⎜⎜ 1 − ⎟⎟ − ay∗ − bz ∗ ⎟ ⎟ ⎜ ⎟⎟ ⎜ ⎝ K⎠ K⎠ ⎝ ⎝ ⎠⎠ ⎝
(
)
⎛ ⎛ ⎛ ⎛ y∗ ⎞ ⎞⎞ y⎞ + y − y∗ ⎜ q ⎜1 − ⎟ + cx − dz − vU − ⎜ q ⎜⎜1 − ⎟⎟ + cx∗ − dz∗ − vU ⎟ ⎟ ⎜ ⎟⎟ ⎜ ⎝ L⎠ L⎠ ⎝ ⎝ ⎠⎠ ⎝ ∗ ⎛ ⎛ ⎛ ⎞⎞ z ⎞ z ⎞ ⎛ + z − z ∗ ⎜ r ⎜1 − ⎟ + ex − fy − wU − ⎜ r ⎜⎜1 − ⎟⎟ + ex∗ − fy∗ − wU ⎟ ⎟ ⎜ ⎟⎟ ⎜ ⎝ M⎠ ⎝ ⎝ M⎠ ⎠⎠ ⎝
(
)
(
)
p q ( x − x∗ ) − a( y − y ∗ ) − b( z − z ∗ )] + ( y − y ∗ )[− ( y − y ∗ ) + c( x − x∗ ) − d ( z − z ∗ )] L K r + ( z − z ∗ )[− ( z − z ∗ ) + e( x − x∗ ) − f ( y − y∗ )] M p q ∗ 2 = −[ ( x − x ) + ( x − x∗ )( y − y∗ )(a − c) + ( y − y ∗ ) 2 + ( y − y∗ )( z − z ∗ )(d + f ) L K r + ( z − z ∗ ) 2 + ( z − z ∗ )( x − x∗ )(b − e)] M
= ( x − x∗ )[−
Persamaan V ′ di atas dapat ditulis dalam bentuk kuadratik
a−c b−e ⎛ p ⎜ K 2 2 ⎜ a−c q d+ f ⎜ A=⎜ L 2 2 ⎜ − + b e d f r ⎜ ⎜ 2 M 2 ⎝ Jika dipilih a = c, b = e
V ′ = − xT Ax
dengan ⎡ x − x∗ ⎤ ⎢ ⎥ x = ⎢ y − y∗ ⎥ ⎢ ∗⎥ ⎢⎣ z − z ⎥⎦
⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
maka matriks A menjadi : ⎛ p ⎜K ⎜ A = ⎜⎜ 0 ⎜ ⎜0 ⎜ ⎝
xT = [( x − x∗ ) ( y − y ∗ ) ( z − z ∗ )] dan
0 q L d+ f 2
⎞ ⎟ ⎟ d+ f ⎟ 2 ⎟ ⎟ r ⎟ M ⎟⎠ 0
Nilai-nilai eigen dari matriks A adalah : p λ1 = , K 2
λ2 = −
)
(
)
8LM 2
λ3 = −
(
−4Mq − 4Lr − (4Mq + 4Lr ) + 16LM d 2 LM + 2dfLM + f 2 LM − 4qr
−4Mq − 4Lr + (4Mq + 4Lr ) + 16LM d 2 LM + 2dfLM + f 2 LM − 4qr 8LM
,
.
Lihat Lampiran 4.
λ1 > 0 (benar, karena dari asumsi p, K > 0 ), λ2 > 0 ⇔ −
−4Mq − 4 Lr −
( 4Mq + 4 Lr )2 + 16 LM ( d 2 LM + 2dfLM + f 2 LM − 4qr ) 8 LM
>0
23
( 4Mq + 4 Lr )2 + 16 LM ( d 2 LM + 2dfLM + f 2 LM − 4qr ) > 0
⇔ 4Mq + 4 Lr +
⇔ ⎛⎜ − ⎝
( 4Mq + 4 Lr )2 + 16 LM ( d 2 LM + 2dfLM + f 2 LM − 4qr ) ⎞⎟
(
⎠
)
2
< ( 4 Mq + 4 Lr )
⇔ ( 4 Mq + 4 Lr ) + 16 LM d LM + 2dfLM + f LM − 4qr < ( 4Mq + 4 Lr ) 2
(
2
2
)
2
2
⇔ 16 LM d 2 LM + 2dfLM + f 2 LM − 4qr < 0
(
)
⇔ LM d 2 + 2df + f 2 < 4qr ⇔ LM ( d + f ) < 4qr 2
⇔
4rq 2 > (d + f ) , LM
λ3 > 0 ⇔ −
−4 Mq − 4 Lr +
f 2 LM − 4qr
8 LM
⇔ −4 Mq − 4 Lr + ⇔
( 4Mq + 4 Lr )2 + 16 LM ( d 2 LM + 2dfLM +
( 4 Mq + 4 Lr )
2
(
(
)
+ 16 LM d 2 LM + 2dfLM + f 2 LM − 4qr < 0
( 4Mq + 4 Lr )2 + 16 LM ( d 2 LM + 2dfLM +
)
f 2 LM − 4qr < 4 Mq + 4 Lr
)
⇔ ( 4 Mq + 4 Lr ) + 16 LM d LM + 2dfLM + f LM − 4qr < ( 4 Mq + 4 Lr ) 2
(
) >0
2
2
)
2
⇔ 16 LM d 2 LM + 2dfLM + f 2 LM − 4qr < 0
(
)
⇔ LM d 2 + 2df + f 2 < 4qr ⇔ LM ( d + f ) < 4qr 2
⇔
4rq 2 > (d + f ) . LM
Dari Teorema 3, A adalah matriks definit positif. V ′ < 0 , jika matriks A adalah definit positif. Akhirnya dengan memakai Teorema 5, titik tetap E4 adalah stabil asimtotik global jika i. a = c, b = e; 4rq 2 ii. > (d + f ) LM Berikut ini diberikan contoh titik-titik tetap sistem persamaan diferensial (4.3.7) beserta kestabilannya. Misalkan nilai-nilai parameter yang diberikan: p = 2.5,
q = 2,
L = 220000, −6
r = 2.5,
K = 350000,
M = 225000,
a = 2.5 x10 −6 ,
−6
b = 1.2 x10 , c = 2.5 x10 , d = 1.55 x10 −5 , e = 1.2 x10 −6 , f = 3.75 x10 −6 , v = 1.47 x10 −4 ,
dan w = 2.57 x10 −4 yang memenuhi syarat Teorema 8, maka kestabilan titik-titik tetapnya adalah:
1. 2. 3. 4. 5. 6. 7. 8.
E0 {0, 0, 0} takstabil (repellor), E1{350000, 0, 0} takstabil (titik sadel), E2 {291154., 168130., 0} takstabil (titik sadel), E3{337777., 0, 72753.7} takstabil (titik sadel), E4 {301407., 126253., 26215.4} stabil asimtotik global, E5 {0, 88062.9, 0} takstabil (titik sadel), E6 {0, 61748.9, 15433.5} takstabil (titik sadel), dan E7 {0, 0, 36273.7} takstabil (titik sadel).
Lihat Lampiran 3.
VI HASIL PEMANENAN MAKSIMUM
Menurut Clark (1976), manajemen sumberdaya yang dapat diperbaharui (seperti perikanan-pen) didasarkan pada konsep maximum sustainable yield (MSY)/hasil pemanenan maksimum. Konsep tersebut didasarkan pada model pertumbuhan biologi yang mengasumsikan jika banyaknya stok dalam populasi lebih rendah dibandingkan dengan suatu level stok K , maka akan terdapat surplus produksi yang dapat dipanen terus menerus tanpa menyebabkan kepunahan populasi. Jika surplus tersebut tidak dipanen maka akan menyebabkan peningkatan level stok menuju daya muat lingkungan K , dimana surplus produksi menjadi nol. Karena surplus produksi adalah sama dengan hasil pemanenan pada setiap level populasi, maka MSY dicapai ketika surplus produksi adalah yang paling besar pada setiap
T=
( L ( beKM + pr ) v
2
level populasi (yaitu pada level dimana laju pertumbuhan populasi adalah maksimal). Ekuilibrium biologis dari model pemanenan perikanan seperti pada (4.3.7) diperoleh dengan menentukan dy dt = 0 dan dz dt = 0 yang berarti bahwa tingkat pertumbuhan bersih dari persediaan populasi akan sama dengan nol, bila tingkat pertumbuhan alami sama dengan hasil pemanenan dari industri perikanan. Ekuilibrium hasil pemanenan T diperoleh dengan menyubstitusi y4 dan z4 persamaan (5.1.7) ke dalam fungsi pemanenan T = vUy4 + wUz4 , (6.1) sehingga diperoleh persamaan berikut
− LM ( bcK + aeK + ( d + f ) p ) vw + M ( acKL + pq ) w2
M ( adeKL + dfLp + bK ( cfL − eq ) ) − ( acKL + pq ) r
)U
2
+
( L ( − p ( cK + q ) r + dMp ( eK + r ) + bKM ( −eq + cr )) vM ( cKL ( fp − ar ) + q ( aeKL − p ( eK − fL + r ) )) w) U M ( adeKL + dfLp + bK ( cfL − eq ) ) − ( acKL + pq ) r
(6.2) Persamaan (6.2) adalah suatu persamaan kuadrat yang memiliki grafik yang berbentuk parabola. Nilai U memaksimumkan T diperoleh dengan
dT = 0 , sehingga diperoleh nilai dU U seperti berikut:
menentukan
U = −dLMp (eK + r )v + L ( p (cK + q ) r + bKM (eq − cr ))v + M (cKL (−fp + ar ) + q (−aeKL + p (eK − fL + r )))w 2 L (beKM + pr )v 2 − 2 LM (bcK + aeK + (d + f ) p )vw + 2 M (acKL + pq )w
2
(6.3) Solusi inilah yang dikenal dengan solusi hasil pemanenan maksimum/maximum sustainable yield (MSY). T TMSY 200000 150000 100000 50000
UMSY 2500
5000
7500
10000 12500 15000
2UMSY U
Gambar 11 kurva MSY dengan menggunakan nilainilai parameter pada halaman 23.
Dari Gambar 11 dapat dilihat hubungan antara T (TMSY ) dan U (U MSY ) , yaitu kenaikan level usaha akan meningkatkan hasil pemanenan hingga suatu level maksimum yaitu pada U MSY , kemudian akan menurun menuju nol yaitu pada 2U MSY . Sedangkan pada saat U ≥ 2U MSY maka hasil pemanenan adalah nol. Selanjutnya dengan menyubstitusi U MSY , y4 dan z4 ke dalam persamaan TMSY = vU MSY y4 + wU MSY z4
(6.4)
akan diperoleh hasil pemanenan maksimum tanpa menggerakkan populasi menuju 24
25
kepunahan. Dengan menyubstitusi nilai-nilai parameternya, dapat diperiksa bahwa
U MSY = U 4 .
VII KEBIJAKAN PEMANENAN OPTIMAL 7.1 Fungsi Objektif Usaha Pemanenan Dalam kegiatannya, perusahaan menjual hasil pemanenan/tingkat output tertentu T yang terdiri dari dua macam yaitu T1 = vUy dan T2 = wUz , dengan harga pasar berturut-turut h1 dan h2 per unit yang diasumsikan konstan. Maka penerimaan total (TR) dapat diketahui, yaitu TR = h1T1 + h2T2 . = h1vUy + h2 wUz . Untuk memproduksi T1 dan T2 , total biaya ekonomi yang ditanggung dapat dinyatakan TC yang merupakan biaya dengan pemanenan per unit usaha B dikalikan dengan tingkat usaha pemanenan U , sehingga TC = BU . Selisih antara total pendapatan dan total biaya disebut laba ekonomi π (sewa ekonomi) yang bergantung pada jumlah output yang diproduksi, yaitu π = TR − TC = h1vUy + h2 wUz − BU = (h1vy + h2 wz − B )U . (7.1.1) Total laba ekonomi di masa yang akan datang (dalam kasus T = ∞ ) yang dihitung dengan present value adalah ∞
π 0 = ∫ π ( y , z , U ) e −δ t dt ∞
(7.1.2)
0
Tujuan usaha pemanenan adalah untuk memaksimumkan present value (nilai sekarang) dari pendapatan pemanenan ikan dengan syarat ikan yang dipanen tidak mengalami kepunahan dan menggunakan sebagai tingkat usaha pemanenan U kontrolnya, sehingga fungsi objektifnya adalah: ∞
Maks ∫ e −δ t (h1vy + h2 wz − B)U dt U
∞
Maks ∫ e−δ t (h1vy + h2 wz − B)U dt U
(7.1.3)
0
terhadap kendala persamaan (4.3.7). Menurut [Tu, 1993] permasalahan di atas adalah masalah kontrol optimum.
(7.2.1)
0
terhadap kendala dx x⎞ ⎛ = px ⎜ 1 − ⎟ − axy − bxz dt K ⎝ ⎠ dy y⎞ ⎛ = qy ⎜ 1 − ⎟ + cxy − dyz − vUy dt ⎝ L⎠ dz z ⎞ ⎛ = r z ⎜1 − ⎟ + exz − fyz − wUz dt M ⎝ ⎠ Permasalahan kontrol optimum ini akan diselesaikan dengan menggunakan Prinsip Maksimum Pontryagin. Untuk kasus horizon dapat waktu takberhingga ( T → ∞ ), ditunjukkan bahwa π = e −δ t (h1vy + h2 wz − B)U memenuhi hipotesis Teorema 7, sehingga Teorema tersebut dapat digunakan (lihat Lampiran 8). Persamaan Hamilton : H = e
0
= ∫ e−δ t (h1vy + h2 wz − B)U dt .
7.2 Penentuan Solusi Optimal Usaha Pemanenan Untuk pemanenan yang optimal dituliskan kembali fungsi objektif pada persamaan (7.1.3) terhadap kendala persamaan (4.3.7) sebagai berikut.
−δ t
( h1 vy + h2 wz − B )U
⎡ p ⎤ + λ1 ( t ) ⎢ − x 2 − axy − bxz + px ⎥ ⎣ K ⎦ ⎡ q ⎤ + λ 2 ( t ) ⎢ − y 2 + cxy − dyz + y ( q − vU ) ⎥ ⎣ L ⎦ ⎡ r ⎤ + λ3 ( t ) ⎢ − z 2 + exz − fyz + z ( r − wU ) ⎥ (7.4) ⎣ M ⎦ dengan λ1 , λ2 , dan λ3 adalah variabel variabel adjoint (co-state). U adalah variabel 0 ≤ U ≤ U maks yang kontrol dengan
mendefinisikan himpunan kontrol Vt = [0, U maks ] . Misalkan U adalah kontrol optimum dengan x, y dan z adalah responrespon yang bersesuaian. Dari Prinsip Maksimum Pontryagin, terdapat variabelvariabel adjoint λ1 , λ2 , dan λ3 untuk setiap t ≥ 0 sedemikian sehingga
26
d λ1 ∂H ⎛ 2 p ⎞ x + ay + bz − p ⎟ λ1 − cyλ2 − ezλ3 =− = dt ∂x ⎜⎝ K ⎠ d λ2 ∂H ⎛ 2q ⎞ y + cx − dz + q ⎟ λ2 + fzλ3 =− = − h1ve−δ t − λ2 v U + axλ1 − ⎜ − dt ∂y ⎝ L ⎠
(
⎫ ⎪ ⎪ ⎪⎪ ⎬ (7.5) ⎪ ⎪ ⎪ ⎭⎪
)
d λ3 ∂H ⎛ 2r ⎞ z + ex − fy + r ⎟ λ3 + dyλ2 =− = − h2 we −δ t − λ3 w U + bxλ1 − ⎜ − dt M ∂z ⎝ ⎠
(
)
Akan ditinjau kondisi optimal pemanenan di titik tetap misalkan E4 . Titik tetap
E4 ( x4 , y4 , z4 ) dipilih karena merupakan titik
interior dan sudah dibuktikan bersifat stabil asimtotik global. Dengan menggunakan d operator diferensial D ≡ , maka sistem dt persamaan (7.5) di titik tetap E4 dapat dituliskan kembali sebagai berikut: px ( D − 4 )λ1 + cy4 λ2 + ez4 λ3 = 0 (7.6) K qy ax4 λ1 − ( D − 4 )λ2 + fz4 λ3 = U 4 h1 v e −δ t L (7.7) bx4 λ1 + dy4 λ2 − ( D −
Sistem persamaan diferensial linear pada persamaan (7.6), (7.7), dan (7.8) dapat diselesaikan dengan menggunakan metode operator. Setelah dilakukan eliminasi terhadap λ2 dan λ3 diperoleh hasil sebagai berikut: (a3 D 3 + a2 D 2 + a1 D + a0 )λ1 = M 1e −δ t dengan a3 = 1 ,
(7.9)
px4 qy4 rz4 ), + + K L M qry4 z4 pqx4 y4 prx4 z4 a1 = + + − dfx4 z4 LM KL KM + bex4 z4 + acx4 y4 , a2 = −(
rz4 )λ3 = U 4 h2 we −δ t M
(7.8)
a0 = −(
pqr dfp beq acr − + + − ade − bcf ) x4 y4 z4 , KLM K L M
M 1 = h1v[−cy4δ + (de −
⎡ ⎤ cr eq ⎞ ⎛ ) y4 z4 ]U 4 + h2 w ⎢ −ez4δ + ⎜ cf − ⎟ y4 z4 ⎥ U 4 , M L ⎝ ⎠ ⎣ ⎦
dan N = −δ 3 + 3a2δ 2 − 3a1δ − a0 ≠ 0 . (Lihat Lampiran 5). Solusi lengkap dari persamaan (7.9) adalah: M (7.10) λ1 = A1e β1 t + A 2 e β 2 t + A3 e β3 t + 1 e −δ t N Ai (i = 1, 2,3) adalah konstanta dengan sembarang, β i (i = 1, 2,3) adalah akar-akar persamaan karakteristik (7.9). Dengan cara yang sama diperoleh: M λ2 = A1e β1 t + A 2 e β2 t + A3 e β3 t + 2 e −δ t (7.11) N M λ3 = A1e β1 t + A 2 e β2 t + A3 e β3 t + 3 e −δ t (7.12) N Sesuai dengan persamaan (2.6.8) Teorema 7, kondisi transversalitas pada t → ∞ , adalah:
lim H ( y ∗ , z ∗ ,U ∗ , λ , t ) = t →∞
lim e −δ t ( h1vy + h2 wz − B ) U t →∞
+
lim λ1 t →∞
dx dt
dy dz + lim λ2 + lim λ3 =0. t →∞ t →∞ dt dt Jika δ > 0 , maka lim e −δ t = 0 . Jika βi < 0 t →∞
atau Ai = 0 ( i = 1, 2,3) , maka lim A1e β1 t + A 2 e β2 t + A3 e β3 t = 0 , t →∞
Akibatnya
kondisi transversalitas dipenuhi jika δ > 0 dan ( β i < 0 atau Ai = 0 , (i = 1, 2,3) ). Untuk penyederhanaan penghitungan, dipilih
27
Ai = 0 (i = 1, 2,3) , maka dari persamaan (7.10), (7.11), dan (7.12) diperoleh M λ1 = 1 e−δ t . (7.13) N
M 2 −δ t e , N M λ3 = 3 e−δ t , N
λ2 =
(7.14) (7.15)
dengan ⎡ ⎤ cr eq ⎞ ⎛ M 1 = h1v[−cy4δ + (de − ) y4 z4 ]U 4 + h2 w ⎢ −ez4δ + ⎜ cf − ⎟ y4 z4 ⎥ U 4 , M L⎠ ⎝ ⎣ ⎦ ⎡ ⎤ px rz pr ⎛ ⎞ ⎛ pf ⎞ + be ⎟) x4 z4 ]U 4 + h2 w ⎢ fz4δ + ⎜ + ae ⎟ x4 z4 ⎥ U 4 , M 2 = − h1v[δ 2 + ( 4 + 4 )δ + ⎜ K M ⎝ KM ⎠ ⎝ K ⎠ ⎣ ⎦ M 3 = h1v[dy4δ + (
⎡ ⎤ qy ⎞ pd ⎛ px ⎛ pq ⎞ + bc) x4 y4 ]U 4 − h2 w ⎢δ 2 + ⎜ 4 + 4 ⎟ δ + ⎜ + ac ⎟ x4 y4 ⎥ U 4 . K K L KL ⎝ ⎠ ⎝ ⎠ ⎣ ⎦
⎫ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎭
(7.16)
Lihat Lampiran 5. Kondisi persamaan Hamilton H maksimum untuk U ∈ Vt adalah ∂H = e −δ t (h1vy + h2 wz − B) − λ2vy − λ3 wz ∂U =0
(7.17)
atau dπ . (7.18) dU Dari kondisi U yang memaksimumkan laba dπ ekonomi π adalah = 0 , maka diperoleh dU M M λ2 vy + λ3 wz = 0 ⇔ 2 vy + 3 wz = 0 N N
λ2 vy + λ3 wz = e −δ t
M 3 wz . (7.19) M 2v Turunan pertama persamaan (7.13), (7.14), dan (7.15) terhadap t adalah: d λ1 −δ M 1 −δ t (7.20) = e , dt N d λ2 −δ M 2 −δ t = e , (7.21) dt N d λ3 −δ M 3 −δ t (7.22) = e . dt N Dengan menyubstitusi persamaan (7.19) (7.20), (7.21), (7.22) ke persamaan (7.5) diperoleh ⇔ y=−
⎞M M wz M wz M −δ M 1 ⎛ 2 p x − a 3 + bz − p ⎟ 1 + c 3 − ez 3 =⎜ N K M v N Nv N 2 ⎝ ⎠ ⎞M M M ⎞ M ⎛ 2q M 3 wz −δ M 2 ⎛ = − ⎜ h1v − 2 v ⎟ U + ax 1 − ⎜ + cx − dz + q ⎟ 2 + fz 3 N N N L M v N N ⎝ ⎠ 2 ⎝ ⎠ ⎞M −δ M 3 M M wz M wz M ⎛ 2r ⎛ ⎞ z + ex + f 3 + r ⎟ 3 − d 3 = − ⎜ h2 w − 3 w ⎟ U + bx 1 − ⎜ − N N N M M v N Nv ⎝ ⎠ 2 ⎝ ⎠
Persamaan (7.23) terdiri dari 3 variabel ( x, z,U ) dan 3 persamaan sehingga bisa dilakukan eliminasi untuk mendapatkan solusi optimal x = xδ , z = zδ , dan U = U δ . Selanjutnya dengan menyubstitusi z = zδ ke persamaan (7.19) diperoleh solusi optimal y = yδ (lihat Lampiran 5).
⎫ ⎪ ⎪ ⎪ ⎪ (7.23) ⎬ ⎪ ⎪ ⎪ ⎪⎭
7.3 Contoh Penerapan pada Industri Perikanan Misalkan dari pengamatan pada suatu perairan diduga parameter – parameter pada persamaan (4.3.7) sebagai berikut : p = 2.5 per tahun q = 2 per tahun r = 2.5 per tahun K = 350000 Metrik Ton (MT) L = 220000 MT M = 225000 MT
28
gambar dinamika mangsa, pemangsa I, dan pemangsa II sebagai berikut
a = 2.5 x10 −6 b = 1.2 x10
−6
c = 2.5 x10 −6 d = 1.55 x10
x 320000
−5
310000
e = 1.2 x10 −6
300000
f = 3.75 x10 −6
290000
−4
v = 1.47 x10 per Standar Pemanenan Harian (SPH) w = 2.57 x10 −4 per SPH Dari data tersebut dapat diperoleh nilai maximum sustainable yield dan level usaha pemanenan maksimum tanpa menggerakkan populasi menuju kepunahan, yaitu : H MSY = 206404 MT/ tahun ,
U MSY = 8159.37 SPH/ tahun . Misalkan diberikan : h1 = $ 700 per MT, h2 = $ 500 per MT, B = $ 2500 per SPH. Untuk δ = 10 % per tahun, persamaan (7.23) menjadi: -168000 + x - 1.74214 z = 0 732063 + U + 0.8447 x - 19.1456 z = 0 -5850.92 + U - 0.00322 x + 0.007603 z = 0
Dengan melakukan substitusi dan eliminasi diperoleh: xδ = 254768 MT, zδ = 49805.6 MT, U δ = 6292.8 SPH/ tahun . Kemudian zδ disubstitusikan ke persamaan (7.19) diperoleh yδ = 106185 MT. Hasil pemanenan optimal total Tδ adalah 178774 MT/ tahun dengan Tδ pemangsa I sebesar 98225.7 MT/ tahun Tδ pemangsa II sebesar dan 80548.1 MT/ tahun , sedangkan keuntungan optimal π δ = $ 9.33 x107 (lihat Lampiran 5). Dengan memasukkan nilai optimal xδ , yδ , zδ sebagai nilai awal, diperoleh
280000 270000 260000 t 5
10
15
20
25
Gambar 12 Dinamika mangsa di E4 . y 130000 120000 110000 100000 90000 t 10
20
30
40
50
Gambar 13 Dinamika pemangsa I di E4 . z 50000 45000 40000 35000 30000 25000 t 10
20
30
40
50
Gambar 14 Dinamika pemangsa II di E4 .
Gambar 12, 13, dan 14 memperlihatkan bahwa pemanenan yang dilakukan akan menggerakkan banyaknya populasi menuju nilai keseimbangannya. Lihat Lampiran 7.
VIII SIMPULAN
Dengan mengonstruksi fungsi Liapunov yang sesuai dan membuat beberapa batasan parameter, sistem interaksi tiga spesies ikan dengan dua spesies yang saling berkompetisi seperti yang dimodelkan pada karya ilmiah ini bersifat stabil asimtotik global pada beberapa wilayah parameter.
Masalah kebijakan pemanenan optimal telah diselesaikan dengan menggunakan Prinsip Maksimum Pontryagin. Dengan memberikan nilai-nilai parameternya, diperoleh level populasi optimal dari pemangsa I dan II adalah 44.65 % dan 35.8 % dari daya muat lingkungan masing-masing populasi.
DAFTAR PUSTAKA
Abell, M.L. and Braselton, J.P. 2004. Differential Equations with Mathematica. 3th Edition. Elsevier Academic Press. San Diego. Anton, H. 2000. Dasar-dasar Aljabar Linear. Jilid 2. Terjemahan Hari Suminto. Interaksara, Batam. Bajpai, A.C., Calus, I.M. and Hyslop, J. 1970. Ordinary Differential Equations. John Wiley & Sons Ltd. England. Chattopadhyay, J., Mukhopadhyay, A., and Tapaswi, P. K.1996. A resource based competitive system in three species fishery. Nonlinear Studies 3(1):73-84. Clark, C.W. 1976. Mathematical Bioeconomics: The Optimal management of renewable resources. John Wiley & Sons, Inc. New York Jaeger, W. K. 1997. Saving Salmon with Fishwheels: A Bioeconomic Analysis. Natural Resources Journal 37:785-808. Fisher, S.D. 1990. Complex Variables. Second Edition. Wadsword & Brooks/ ColeBooks & Software, Pasific Grove. California. Grossman, S.I. 1995. Multivariable Calculus, Linear Algebra, and Differential Equations. Saunders College Publishing. New York. Hoffmann, LD and Bradley, GL. 1989. Calculus for Business, Economics and the Social and Life Sciences. 4th edition. McGraw-Hill. Singapore.
Terjemahan Ign Bayu Mahendra dan Abdul Aziz. Erlangga, Jakarta. Purcell, E.J. dan Varberg, D. 1999. Kalkulus dan Geometri Analitis. Edisi ke-5. Terjemahan I Nyoman Susila dan kawankawan. Erlangga, Jakarta. Seirstad A. 2002. Condition implying the vanishing of the Hamiltonian at the infinite horizon in optimal control problems. Memorandum no. 22. Departement of Economics University of Oslo. http://www.oekonomi.uio.no/memo/. [14 Feb 2007]. Shone, R. 1997. Economic Dynamics : Phase Diagrams and Their Economic Applications. Cambridge : Cambridge University Press. Smith, P. 1987. Nonlinear Ordinary Differential Equations. Clarendon Press. Oxford. Stewart, J. 1991. Single Variable Calculus. 2nd Edition. Brooks/Cole Publishing. Pacific Grove, California. Szidarovszky, F. and Bahill, A. T. 1998. Linear System Theory. 2nd Edition. CRC Press LLC. Florida. Tu, P. N. V. 1993. Introductory Optimization Dynamics: Optimal Control with Economics and Management Applications. Second Revised and Enlarged Edition. Springer Verlag. Berlin.
Kreyszig, E. 1983. Advanced Engineering Mathematics. Fifth Edition. John Wiley and Sons. New York.
Tu, P. N. V. 1994. Dynamical Systems : An Introduction with Applications in Economics and Biology. Second Revised and Enlarged Edition. Springer Verlag. Berlin.
Leon, T.J. 2001. Aljabar Linear dan Aplikasinya. Edisi ke-5. Terjemahan Alit Bondan. Erlangga, Jakarta.
Vandermeer, J. 1981. Elementary Mathematical Ecology. John Wiley & Sons, Inc. New York.
Nababan, M. 1988. Pengantar Matematika untuk Ilmu Ekonomi dan Bisnis. Erlangga, Jakarta.
Velhurst, F. 1990. Nonlinear Differential Equation and Dynamical System. Springer-Verlag. Heidelberg.
Nicholson, W. 2002. Mikroekonomi Intermediate dan Aplikasinya. Edisi ke-8.
Zipse, P. W. 1987. Introduction to Mathematical Analysis. McGraw-Hill Book Company. Singapore.
29
LAMPIRAN
Lampiran 1 (Titik Tetap)
Hasil ini diperoleh dengan menggunakan fasilitas Dynpac yang dapat di download di www.wolframresearch.com kemudian dijalankan pada software Mathematica 5.1. In[328]:= sysid
Mathematica 5.1.0, DynPac 10.69, 5ê 24ê 2007
In[329]:= "Mencari titik keseimbangan "; In[330]:= intreset; plotreset;
In[331]:= setstate@8x, y, z
In[332]:= setparm@8 p, q, r, K, L, M, a, b, c, d, e, f, v, U
In[333]:= slopevec = : p x jj1 k
xy z
y y zy i i z - a x y - b x z , q y jj1 - zz +c x y - d y z - v U y, r z jj1 zz +e x z - f y z - w U z>; { k { k K L M{
In[334]:= eqstates = findpolyeq; In[335]:= E0 = eqstates@@1DD
Out[335]= 80, 0, 0<
In[336]:= E1 = FullSimplify @eqstates@@2DDD
Out[336]= 8K, 0, 0<
In[337]:= E2 = FullSimplify @eqstates@@3DDD Out[337]= :
K H p q +a L H-q +U vLL acKL+p q
,
L p Hc K +q - U vL acKL+p q
, 0>
In[338]:= E3 = FullSimplify @eqstates@@4DDD Out[338]= :
K H p r + b M H- r +U wLL b eKM+p r
, 0,
M p He K + r - U wL b eKM+p r
>
In[339]:= E4 = FullSimplify @eqstates@@5DDD Out[339]= :
K H r H-a L q + p q +a L U vL - d L M H f p - a r +a U wL + b M H-q r + f L Hq - U vL +q U wLL -L M Ha d e K + b c f K +d f pL + b e K M q +Ha c K L + p qL r
L H p r Hc K +q - U vL - d M p He K + r - U wL + b K M He q - c r - e U v +c U wLL -L M Ha d e K + b c f K +d f pL + b e K M q +Ha c K L + p qL r
,
M Ha e K L Hq - U vL +c K L H f p - a r +a U wL - p He K q - f L q +q r + f L U v - q U wLL M Ha d e K L +d f L p + b K Hc f L - e qLL - Ha c K L + p qL r
,
>
In[340]:= E5 = FullSimplify @eqstates@@6DDD Out[340]= :0, L -
LU v q
, 0>
In[341]:= E6 = FullSimplify @eqstates@@7DDD Out[341]= :0,
L H-q r + r U v +d M H r - U wLL d fLM-q r
,
M H f L Hq - U vL +q H-r +U wLL d fLM-q r
>
In[342]:= E7 = FullSimplify @eqstates@@8DDD Out[342]= :0, 0, M -
MU w r
>
31
Lampiran 2 (Nilai Eigen)
Hasil ini diperoleh dengan menggunakan fasilitas Dynpac yang dapat di download di www.wolframresearch.com kemudian dijalankan pada software Mathematica 5.1. In[343]:= sysid
Mathematica 5.1.0, DynPac 10.69, 5ê 24ê 2007
In[344]:= intreset; plotreset;
In[345]:= setstate@8x, y, z
In[346]:= setparm@8 p, q, r, K, L, M, a, b, c, d, e, f, v, w
In[347]:= slopevec = : p x jj1 k
xy z
yy zy i i z - a x y - b x z , q y jj1 - zz +c x y - d y z - v U y, r z jj1 - zz +e x z - f y z - w U z>; k k L{ M{
K{
In[348]:= eqstates = findpolyeq; In[349]:= E0 = eqstates@@1DD;
In[350]:= E1 = FullSimplify @eqstates@@2DDD; In[351]:= E2 = FullSimplify @eqstates@@3DDD; In[352]:= E3 = FullSimplify @eqstates@@4DDD; In[353]:= E4 = FullSimplify @eqstates@@5DDD; In[354]:= E5 = FullSimplify @eqstates@@6DDD; In[355]:= E6 = FullSimplify @eqstates@@7DDD; In[356]:= E7 = FullSimplify @eqstates@@8DDD;
Konstruksi Matrik Jacobi In[357]:= jac = Simplify @88D@slopevec@@1DD, xD, D@slopevec@@1DD, yD, D@slopevec@@1DD, zD< , 8D@slopevec@@2DD, xD, D@slopevec@@2DD, yD, D@slopevec@@2DD, zD< , 8D@slopevec@@3DD, xD, D@slopevec@@3DD, yD, D@slopevec@@3DD, zD<
Out[358]//MatrixForm= ij p - 2 p x - a y - b z -a x jj K j jj jj c y jj jj jj e z k
yz zz zz 2q y zz - d z -d y q - U v +c x zz L zz z -f z r - U w +e x - f y - 2 r z zz M {
-b x
Pada titik tetapnya berlaku Kontruksi Matriks Jacobi untuk titik tetap E0
In[359]:= jac0 = jac ê. 8x Ø E0@@1DD, y Ø E0 @@2DD, z Ø E0@@3DD< ; In[360]:= "JHE0L=" MatrixForm @ jac0D
0 ij p 0 yz zz Out[360]= JHE0L= jjj 0 q - U v 0 z k0 0 r-U w{
32
Kontruksi Matriks Jacobi untuk titik tetap E1
In[361]:= jac1 = jac ê. 8x Ø E1@@1DD, y Ø E1 @@2DD, z Ø E1@@3DD< ; In[362]:= "JHE1L=" MatrixForm @ jac1D
-bK ij - p -a K yz zz Out[362]= JHE1L= jjj 0 c K +q - U v 0 z k0 0 eK+r-U w{
Kontruksi Matriks Jacobi untuk titik tetap E2
In[363]:= jac2 = FullSimplify @ jac ê. 8x Ø E2@@1DD, y Ø E2 @@2DD, z Ø E2@@3DD
Kontruksi Matriks Jacobi untuk titik tetap E3
In[365]:= jac3 = FullSimplify @ jac ê. 8x Ø E3@@1DD, y Ø E3 @@2DD, z Ø E3@@3DD
Kontruksi Matriks Jacobi untuk titik tetap E4
In[367]:= jac4 = FullSimplify @ jac ê. 8x Ø x4, y Ø y4 , z Ø z4
Kontruksi Matriks Jacobi untuk titik tetap E5
In[369]:= jac5 = FullSimplify @ jac ê. 8x Ø E5@@1DD, y Ø E5 @@2DD, z Ø E5@@3DD
i ij U v yz 0 yz 0 jj p +a L j-1 + zz q { k jj zz jj zz d L H- q+U vL jj c ijjL - L U v yzz zz -q +U v Out[370]= JHE5L= j zz jj k q { q zz jjj zz i y jj 0 0 r + f L j-1 + U v z - U w zz q { k k {
Kontruksi Matriks Jacobi untuk titik tetap E6 In[371]:= jac6 = FullSimplify @ jac ê. 8x Ø E6@@1DD, y Ø E6 @@2DD, z Ø E6@@3DD
Kontruksi Matriks Jacobi untuk titik tetap E7
In[373]:= jac7 = FullSimplify @ jac ê. 8x Ø E7@@1DD, y Ø E7 @@2DD, z Ø E7@@3DD
ij p + b M J -1 + U w N 0 yz 0 jj zz r jj zz U w jj 0 zz q - U v +d M J -1 + N 0 Out[374]= JHE7L= j zz jj r zz jj z f M H- r+U wL jj e J M - MU w N - r +U w zz k { r r
33
Nilai Eigen di titik tetap E0
In[375]:= l0 = FullSimplify @Eigenvalues@ jac0DD;
In[376]:= 8"l1="@l0@@1DDD, "l2 ="@l0@@2DDD, "l3="@l0@@3DDD<
Out[376]= 8l1=@ pD, l2 =@q - U vD, l3=@ r - U wD<
Nilai Eigen di titik tetap E1
In[377]:= l1 = FullSimplify @Eigenvalues@ jac1DD;
In[378]:= 8"l1="@l1@@1DDD, "l2 ="@l1@@2DDD, "l3="@l1@@3DDD<
Out[378]= 8l1=@- pD, l2=@c K +q - U vD, l3 =@e K + r - U wD<
Nilai Eigen di titik tetap E2
In[379]:= l2 = FullSimplify @Eigenvalues@ jac2DD;
In[380]:= 8"l1="@l2@@1DDD, "l2 ="@l2@@2DDD, "l3="@l2@@3DDD< Out[380]= :l1=B
1
2 Ha c K L + p qL
I-c K p q +a L p Hq - U vL - p q H p +q - U vL +
, I p I p Hq Hc K - a L + p +qL +Ha L - qL U vL2 - 4 Ha c K L + p qL Hc K +q - U vL H p q +a L H-q +U vLLMMMF ,
l2=B -
1
2 Ha c K L + p qL
Ic K p q - a L p Hq - U vL + p q H p +q - U vL +
, I p I p Hq Hc K - a L + p +qL +Ha L - qL U vL2 - 4 Ha c K L + p qL Hc K +q - U vL H p q +a L H-q +U vLLMMMF ,
l3=B r +
-c f K L p +a e K L H-q +U vL + p He K q - f L q + f L U vL
acKL+pq
- U wF>
Nilai Eigen di titik tetap E3
In[381]:= l3 = FullSimplify @Eigenvalues@ jac3DD;
In[382]:= 8"l1="@l3@@1DDD, "l2 ="@l3@@2DDD, "l3="@l3@@3DDD< Out[382]= :l1=B q + l2=B
p r Hc K - U vL - d M p He K + r - U wL - b K M Hc r +e U v - c U wL beKM+p r
1
2 H b e K M + p rL
F,
I-e K p r + b M p H r - U wL - p r H p + r - U wL +
, I p I p H r He K - b M + p + rL +H b M - rL U wL2 - 4 H b e K M + p rL He K + r - U wL H p r + b M H- r +U wLLMMMF ,
l3=B -
1
2 H b e K M + p rL
Ie K p r - b M p H r - U wL + p r H p + r - U wL +
, I p I p H r He K - b M + p + rL +H b M - rL U wL2 - 4 H b e K M + p rL He K + r - U wL H p r + b M H- r +U wLLMMMF>
Nilai Eigen di titik tetap E4
In[383]:= l4 = FullSimplify @Eigenvalues@ jac4DD;
In[384]:= 8"l1="@l4@@1DDD, "l2 ="@l4@@2DDD, "l3="@l4@@3DDD< ;
34
Nilai Eigen di titik tetap E5
In[385]:= l5 = FullSimplify @Eigenvalues@ jac5DD;
In[386]:= 8"l1="@l5@@1DDD, "l2 ="@l5@@2DDD, "l3="@l5@@3DDD< i
Out[386]= :l1=@-q +U vD, l2=B p +a L jjj-1 + k
U v yz
i U v yz zzF , l3=B r + f L jjj-1 + zz - U wF> q { q { k
Nilai Eigen di titik tetap E6
In[387]:= l6 = FullSimplify @Eigenvalues@ jac6DD;
In[388]:= 8"l1="@l6@@1DDD, "l2 ="@l6@@2DDD, "l3="@l6@@3DDD< Out[388]= :l1=B l2=B
r H-p q +a L Hq - U vLL +d L M H f p - a r +a U wL + b M H- f L q +q r + f L U v - q U wL d fLM-q r
1
F,
I- f L r Hq - U vL +d M q H-r +U wL +q r Hq + r - U H v + wLL +
, I4 Hd f L M - q rL H-q r + r U v +d M H r - U wLL H f L Hq - U vL +q H- r +U wLL +
2d fLM-2q r
H f L r Hq - U vL +q Hd M H r - U wL - r Hq + r - U H v + wLLLL MMF ,
2
l3=B
1
I- f L r Hq - U vL +d M q H-r +U wL +q r Hq + r - U H v + wLL -
, I4 Hd f L M - q rL H-q r + r U v +d M H r - U wLL H f L Hq - U vL +q H- r +U wLL +
2d fLM-2q r
H f L r Hq - U vL +q Hd M H r - U wL - r Hq + r - U H v + wLLLL2MMF>
Nilai Eigen di titik tetap E7
In[389]:= l7 = FullSimplify @Eigenvalues@ jac7DD;
In[390]:= 8"l1="@l7@@1DDD, "l2 ="@l7@@2DDD, "l3="@l7@@3DDD< i
Out[390]= :l1=@- r +U wD, l2=B p + b M jj-1 + k
Uwy
Uwy zzF , l =B q - U v +d M ijj-1 + zzF> 3 k r { r {
In[391]:= Unprotect@ p, q, r, K, L, M, a, b, c, d, e, f, v, wD; In[392]:= ClearAll @ p, q, r, K, L, M, a, b, c, d, e, f, v, wD
35
Lampiran 3 (Evaluasi Nilai-nilai Parameter)
Hasil ini diperoleh dengan menggunakan fasilitas Dynpac yang dapat di download di www.wolframresearch.com kemudian dijalankan pada software Mathematica 5.1. In[393]:= sysid
Mathematica 5.1.0, DynPac 10.69, 5ê 24ê 2007
In[394]:= intreset; plotreset;
In[395]:= setstate@8x, y, z
In[396]:= setparm@8 p, q, r, K, L, M, a, b, c, d, e, f, v, w
In[398]:= slopevec = : p x jj1 k
xy z
yy zy i i zz +e x z - f y z - w U z>; z - a x y - b x z , q y jj1 - zz +c x y - d y z - v U y, r z jj1 k k K{ L{ M{
In[399]:= eqstates = findpolyeq; In[400]:= U =
H-d L M p He K + rL v +L H p Hc K +qL r + b K M He q - c rLL v +
M Hc K L H-f p +a rL +q H-a e K L + p He K - f L + rLLL wLë 2
I2 L H b e K M + p rL v - 2 L M H b c K +a e K +Hd + fL pL v w +2 M Ha c K L + p qL w M;
2
In[401]:= U = evalparm@UD Out[401]= 8159.37 In[402]:= E0 = evalparm@eqstates@@1DDD
Out[402]= 80, 0, 0<
Nilai Eigen di titik tetap E0
In[403]:= l0 = FullSimplify @eigval@E0DD; 8"l1="@l0@@1DDD, "l2 ="@l0@@2DDD, "l3="@l0@@3DDD<
Out[403]= [email protected] D, l2 [email protected] D, [email protected] D< In[404]:= classify@E0D
unstable
In[405]:= E1 = evalparm@FullSimplify @eqstates@@2DDDD
Out[405]= 8350000, 0, 0<
Nilai Eigen di titik tetap E1
In[406]:= l1 = eigval@E1D; 8"l1="@l1@@1DDD, "l2 ="@l1@@2DDD, "l3="@l1@@3DDD<
Out[406]= [email protected] D, [email protected] D, l3 [email protected] D< In[407]:= classify@E1D
unstable
36
In[408]:= E2 = evalparm@FullSimplify @eqstates@@3DDDD
Out[408]= 8291154., 168130., 0<
Nilai Eigen di titik tetap E2
In[409]:= l2 = FullSimplify @eigval@E2DD; 8"l1="@l2@@1DDD, "l2 ="@l2@@2DDD, "l3="@l2@@3DDD<
Out[409]= [email protected] +0.479572 ÂD, [email protected] - 0.479572 ÂD, [email protected] D< In[410]:= classify@E2D
unstable
In[411]:= E3 = evalparm@FullSimplify @eqstates@@4DDDD
Out[411]= 8337777., 0, 72753.7 <
In[412]:= l3 = FullSimplify @eigval@E3DD; 8"l1="@l3@@1DDD, "l2 ="@l3@@2DDD, "l3="@l3@@3DDD<
Out[412]= [email protected] D, [email protected] D, [email protected] D< In[413]:= classify@E3D
unstable
In[414]:= E4 = evalparm@FullSimplify @eqstates@@5DDDD
Out[414]= 8301407., 126253., 26215.4 <
Nilai Eigen di titik tetap E4
In[415]:= l4 = FullSimplify @eigval@E4DD; 8"l1="@l4@@1DDD, "l2 ="@l4@@2DDD, "l3="@l4@@3DDD<
Out[415]= [email protected] +0.286408 ÂD, [email protected] - 0.286408 ÂD, [email protected] D< In[416]:= classify@E4D
strictly stable
In[417]:= E5 = evalparm@FullSimplify @eqstates@@6DDDD
Out[417]= 80, 88062.9, 0<
Nilai Eigen di titik tetap E5
In[418]:= l5 = FullSimplify @eigval@E5DD; 8"l1="@l5@@1DDD, "l2 ="@l5@@2DDD, "l3="@l5@@3DDD<
Out[418]= [email protected] D, l2 [email protected] D, l3 [email protected] D< In[419]:= classify@E5D
unstable
37
In[420]:= E6 = evalparm@FullSimplify @eqstates@@7DDDD
Out[420]= 80, 61748.9, 15433.5 <
Nilai Eigen di titik tetap E6
In[421]:= l6 = FullSimplify @eigval@E6DD; 8"l1="@l6@@1DDD, "l2 ="@l6@@2DDD, "l3="@l6@@3DDD<
Out[421]= [email protected] D, l2 [email protected] D, l3 [email protected] D< In[422]:= classify@E6D
unstable
In[423]:= E7 = evalparm@FullSimplify @eqstates@@8DDDD
Out[423]= 80, 0, 36273.7 <
Nilai Eigen di titik tetap E7
In[424]:= l7 = FullSimplify @eigval@E7DD; 8"l1="@l7@@1DDD, "l2 ="@l7@@2DDD, "l3="@l7@@3DDD<
Out[424]= [email protected] D, l2 [email protected] D, l3 [email protected] D< In[425]:= classify@E7D
unstable
In[426]:= Unprotect@ p, q, r, K, L, M, a, b, c, d, e, f, v, w, UD; In[427]:= ClearAll @ p, q, r, K, L, M, a, b, c, d, e, f, v, w, UD
38
Lampiran 4 (Kestabilan Global)
Hasil ini diperoleh dengan menggunakan software Mathematica 5.1.
* i p Hx - x L y i = jjHx - x* L jj- a H y - y* L - b Hz - z* Lzz + k { k dt K * i q Hy - y L y H y - y* L jj+c Hx - x* L - d Hz - z* Lzz + k { L * yy * i r Hz - z L Hz - z L jj+e Hx - x* L - f H y - y* Lzz zz; {{ k M
dV
In[328]:= X = 88x - x* < , 8 y - y* < , 8z - z* << ; In[329]:= X êê MatrixForm Out[329]//MatrixForm= * ij x - x yz jj * z jj y - y zzz * k z-z { In[330]:= A = ::
p a-c b-e a-c q d +f b -e d+f r , >, : >, : >>; , , , , , K 2 2 2 L 2 2 2 M
In[331]:= MatrixForm @AD Out[331]//MatrixForm= a- c ij p jj K 2 jj q jj a- c jj L jj 2 jj b- e d+ f k 2 2
b- e 2 d+f 2 r M
yz zz zz zz zz zz zz {
In[332]:= XT = Transpose @XD; In[333]:= MatrixForm @XT D
Out[333]//MatrixForm= H x - x* y - y*
z - z* L
39
Mensubstitusikan nilai a=c dan b=e ke dalammatriks A In[334]:= A = A ê. 8a Ø c, b Ø e< Out[334]= ::
p K
, 0, 0>, :0,
In[335]:= MatrixForm @AD Out[335]//MatrixForm= p ij 0 jj K jj q jjj 0 jj L jj d+ f j0 k 2
q L
,
d+f 2
>, :0,
d+f 2
,
r M
>>
yz zz z d+ f zz z 2 zzz r zzz M {
0
Mencari Nilai Eigen dari matriks A In[336]:= l = Eigenvalues@AD; In[337]:= "l1="@l@@1DDD
Out[337]= l1=B
p K
F
In[338]:= "l2="@l@@2DDD Out[338]= l2=B -
1
I-4 M q - 4 L r -
, IH4 M q +4 L rL2 +16 L M Id2 L M +2 d f L M + f2 L M - 4 q rMMMF
8LM
In[339]:= "l3="@l@@3DDD Out[339]= l3=B -
1
I-4 M q - 4 L r +
, IH4 M q +4 L rL2 +16 L M Id2 L M +2 d f L M + f2 L M - 4 q rMMMF
8LM
40
Lampiran 5 (Populasi Optimal)
Hasil ini diperoleh dengan menggunakan software Mathematica 5.1. In[328]:= Unprotect@x, y, z, p, ND; In[329]:=
In[330]:=
ClearAll@p, q, r, K, L, M, a, b, c, d, e, f, v, w, U4, x, y, z, x4, y4, z4, d, ans1, ans2, ans3, p, h1, h2, B, l1, l2, l3, M1, M2, M3, N, ans1, ans2, ans3D; H = Exp@- t dD Hh1 vy + h2wz - BL U + l1 Jl2 J-
p 2 x - axy - bxz + p xN +
K
q 2 r 2 y + cxy - dy z + y Hq - vULN + l3 Jz + exz - f y z + z Hr - wULN; L M i
p x4 y
In[331]:=
slopecostate1 = jjD -
In[332]:=
slopecostate2 = ax4 l1 - jjD -
In[333]:=
slopecostate3 = bx4 l1 - jjD -
k
zz l + cy4 l + ez4 l == 0; 1 2 3
K { i
qy4 y
i
r z4 y
k
k
zz l + f z4 l == U4h1 vExp@-d tD; 2 3
L {
zz l + dy4 l ã U4h2wExp@-d tD; 3 2
M {
Menyelesaikanpersamaan l1 denganmenggunakanoperator - D; In[334]:=
Simplify@Solve@Eliminate@8slopecostate1, slopecostate2, slopecostate3<, 8l2, l3
Out[334]= 99l1 Ø H‰
KU4 HeM HDh2Lw + dh1L vy4 - h2qwy4L z4 + cLy4 HDh1M v - h1r vz4 + f h2M wz4LLL ë 3 ID KLM + H bcf KLM + df LM p - beKM q - p qr + aKL HdeM - crLL x4y4 z4 - D2 HLM p x4 + KM qy4 + KLr z4L + D HacKLM x4y4 + M p qx4y4 + beKLM x4z4 + Lp r x4z4 - df KLM y4z4 + Kqr y4z4LM== -t d
Misalkan l1=
A1 B1
, dengan ;
In[335]:=
A1 = -t d H‰ KU4 HeM HDh2Lw + dh1L vy4 - h2qwy4L z4 + cLy4 HDh1M v - h1r vz4 + f h2M wz4LLL;
In[336]:=
B1 = 3 ID KLM + H bcf KLM + df LM p - beKM q - p qr + aKL HdeM - crLL x4y4z4 D2 HLM p x4 + KM qy4 + KLr z4L + D HacKLM x4y4 + M p qx4y4 + beKLM x4z4 + Lp r x4z4 df KLM y4z4 + Kqr y4z4LM;
Untuk penyederhanaan, kalikan A1 & B1 dengan Mengalikan A1 dengan
In[337]:= Out[337]=
A1 = ExpandBA1 *
1 K LM 1
KLM
1 K LM
;
;
F
cD ‰- t d h1U4 vy4 + De ‰- t d h2U4wz4 + de ‰- t d h1U4 vy4z4 c ‰- t d h1r U4 vy4z4 e ‰- t d h2qU4wy4z4 + c ‰- t d f h2U4wy4z4 M L
41
In[338]:= In[339]:= Out[339]=
"Mengevaluasi operator-D diperoleh:"; A1 = A1 ê. D Ø -d
c ‰- t d h1r U4 vy4z4
de ‰- t d h1U4 vy4z4 -
M
e ‰- t d h2qU4wy4z4 L In[340]:=
Out[340]=
C1 = M1;
In[342]:=
"Sehingga";
In[343]:=
A1 = M1 ‰- t d ;
Out[344]=
D3 -
D2 p x4 K
bDex4z4 +
-
FF
+ cf h2U4wy4z4 -
- ch1U4 vy4 d - eh2U4wz4 d
Mengalikan B1 dengan B1 = ExpandBB1 *
‰- t d
M
L In[341]:=
A1
ch1r U4 vy4z4
eh2qU4wy4z4
In[344]:=
- c ‰- t d h1U4 vy4 d - e ‰- t d h2U4wz4 d
M1 = ExpandBFullSimplifyB deh1U4 vy4z4 -
+ c ‰- t d f h2U4wy4z4 -
1 K LM 1
KLM
D2 qy4
L Dp r x4z4
;
F
+ acDx4y4 +
Dp qx4y4 KL Dqr y4z4
-
D2 r z4 M
+
- dDf y4z4 + + adex4y4z4 + KM LM df p x4y4z4 beqx4y4z4 acr x4y4z4 p qr x4y4z4 bcf x4y4z4 + K L M KLM
"Daftar koefisien persamaan B1 terhadap operator-D"; In[345]:=
Coeff = CoefficientList@B1, DD
Out[345]= :adex4y4z4 +
df p x4y4z4
In[347]:=
K
-
beqx4y4z4
L acr x4y4z4 p qr x4y4z4 p qx4y4 + bex4z4 + , acx4y4 + M KLM KL p r x4z4 qr y4z4 p x4 qy4 r z4 - df y4z4 + ,, 1> KM LM K L M
In[346]:=
bcf x4y4z4 +
"Misalkan,";
a0 = Collect@Coeff @@1DD, x4y4z4D
i Out[347]= jjade + k
bcf +
df p K
-
beq L
-
acr M
-
p qr y zz x4y4z4 KLM {
42
In[348]:= Out[348]=
In[349]:= Out[349]=
In[350]:= Out[350]= In[351]:= In[352]:=
a1 = Coeff @@2DD p qx4y4
acx4y4 +
KL
+ bex4z4 +
p r x4z4 KM
a2 = Coeff @@3DD -
p x4 K
-
qy4 L
a3 = Coeff @@4DD
-
- df y4z4 +
qr y4z4 LM
r z4 M
1
Unprotect@ ND;
N = a3 D3 + a2 D2 + a1 D + a0 ê. D Ø -d
i Out[352]= jjade +
p qr y zz x4y4z4 K L M KLM { k p qx4y4 p r x4z4 qr y4z4 y ij zz d + + bex4z4 + - df y4z4 + jacx4y4 + KL KM LM { k qy4 r z4 y 2 3 ij p x4 zz d - d jL M { k K bcf +
df p
-
beq
acr
-
-
"Persamaan diferensial nya dapat ditulis sebagai:" ; Ha3D +a2D +a1D+a0 Ll1= M1 ‰ 3
2
-t d
;
"Dengan menggunakan metode operator-D untuk kasus tersebut diperoleh:"; Auxiliary Equation adalah a3m3 +3a2m2 +3a1m+a0 =0, dengan akar-akar biHi=1,2,3L ; Complementary Function HC FL adalah A1‰b 1 t+A2 ‰b 2 t +A3‰b 3 t, dengan Ai Hi=1,2,3Ladalah konstanta sembarang;
Particular Integral HPIL adalah In[353]:=
M1 ‰- t d B1
;
"Solusi Umum l1 = C F + P I ";
l1 = A1 ‰b 1 t + A2 ‰b 2 t + A3 ‰b 3 t +
M1 ‰- t d N
Menyelesaikanpersamaan l2 denganmenggunakanoperator - D; In[354]:= In[355]:=
Unprotect@x, y, z, p, ND;
ClearAll@p, q, r, K, L, M, a, b, c, d, e, f, v, w, U, x, y, z, x4, y4, z4, d, ans1, ans2, ans3, p, h1, h2, B, l1, l2, l3, M1, M2, M3, N, ans1, ans2, ans3D;
43
In[356]:=
Simplify@Solve@Eliminate@8slopecostate1, slopecostate2, slopecostate3<, 8l1, l3
Out[356]= 99l2 Ø -I‰
LU4 ID2 h1KM v + H beh1KM v + h1p r v - h2M HaeK + f pL wL x4z4 + D Hf h2KM wz4 - h1 v HM p x4 + Kr z4LLMM ë 3 ID KLM + H bcf KLM + df LM p - beKM q - p qr + aKL HdeM - crLL x4y4z4 - D2 HLM p x4 + KM qy4 + KLr z4L + D HacKLM x4y4 + M p qx4y4 + beKLM x4z4 + Lp r x4z4 - df KLM y4z4 + Kqr y4z4LM== -t d
Misalkan l2=
A2 B2
, dengan ;
In[357]:=
A2 = -I‰- t d LU4 ID2 h1KM v + H beh1KM v + h1p r v - h2M HaeK + f pL wL x4z4 + D Hf h2KM wz4 - h1 v HM p x4 + Kr z4LLMM;
In[358]:=
B2 = 3 ID KLM + H bcf KLM + df LM p - beKM q - p qr + aKL HdeM - crLL x4y4z4 D2 HLM p x4 + KM qy4 + KLr z4L + D HacKLM x4y4 + M p qx4y4 + beKLM x4z4 + Lp r x4z4 df KLM y4z4 + Kqr y4z4LM;
Untuk penyederhanaan, kalikan A2 & B2 dengan Mengalikan A dengan
In[359]:=
Out[359]=
A2 = ExpandBA2 *
1 K LM 1
KLM
-D2 ‰- t d h1U4 v +
1 K LM
;
;
F
D ‰- t d h1p U4 vx4 K
+
D ‰- t d h1r U4 vz4
- D ‰- t d f h2U4wz4 - be ‰- t d h1U4 vx4z4 M ‰- t d h1p r U4 vx4z4 ‰- t d f h2p U4wx4z4 + ae ‰- t d h2U4wx4z4 + KM K
"Mengevaluasi operator-D diperoleh:"; In[360]:= Out[360]=
A2 = A2 ê. D Ø -d
- be ‰- t d h1U4 vx4z4 -
ae ‰- t d h2U4wx4z4 + ‰- t d h1r U4 vz4 d
M
‰- t d h1p r U4 vx4z4
KM t d ‰ f h2p U4wx4z4 K
+ -
‰- t d h1p U4 vx4 d
K
-
+ ‰- t d f h2U4wz4 d - ‰- t d h1U4 v d2
44
In[361]:=
Out[361]=
M2 = ExpandBFullSimplifyB - beh1U4 vx4z4 -
h1p U4 vx4 d
In[362]:=
h1p r U4 vx4z4
KM h1r U4 vz4 d
-
K
FF
A2 ‰- t d
M
+ aeh2U4wx4z4 +
f h2p U4wx4z4 K
-
+ f h2U4wz4 d - h1U4 v d2
C2 = M2;
"Sehingga"; In[363]:=
A2 = M2 ‰- t d ;
Mengalikan B2 dengan
In[364]:=
Out[364]=
B2 = ExpandBB2 * D3 -
D2 p x4
-
K
bDex4z4 +
1 K LM 1
KLM
D2 qy4
L Dp r x4z4
;
F
+ acDx4y4 +
Dp qx4y4 KL Dqr y4z4
-
D2 r z4 M
+
- dDf y4z4 + + adex4y4z4 + KM LM df p x4y4z4 beqx4y4z4 acr x4y4z4 p qr x4y4z4 bcf x4y4z4 + K L M KLM
"Daftar koefisien persamaan B2 terhadap operator-D"; In[365]:=
Coeff2 = CoefficientList@B2, DD
Out[365]= :adex4y4z4 +
bcf x4y4z4 +
df p x4y4z4
beqx4y4z4
L acr x4y4z4 p qr x4y4z4 p qx4y4 + bex4z4 + , acx4y4 + M KLM KL p r x4z4 qr y4z4 p x4 qy4 r z4 - df y4z4 + ,, 1> KM LM K L M
K
-
"Misalkan,"; In[366]:=
a0 = Collect@Coeff2@@1DD, x4y4z4D
i Out[366]= jjade + k
In[367]:= Out[367]=
bcf +
df p
a1 = Coeff2@@2DD acx4y4 +
K
p qx4y4 KL
-
beq L
-
acr M
+ bex4z4 +
-
p qr y z
z x4y4z4
KLM {
p r x4z4 KM
- df y4z4 +
qr y4z4 LM
45
In[368]:= Out[368]=
In[369]:= Out[369]=
a2 = Coeff2@@3DD -
p x4 K
-
qy4 L
-
a3 = Coeff2@@4DD
r z4 M
1
"Persamaan diferensial nya dapat ditulis sebagai:" ; Ha3D +a2D +a1D+a0 Ll2= M2 ‰ 3
2
-t d
;
"Dengan menggunakan metode operator-D untuk kasus tersebut diperoleh:"; Auxiliary Equation adalah a3m3 +3a2m2 +3a1m+a0 =0, dengan akar-akar biHi=1,2,3L ; Complementary Function HC FL adalah A1‰b 1 t+A2 ‰b 2 t +A3‰b 3 t, dengan Ai Hi=1,2,3Ladalah arbitrary constants ; Particular Integral HPIL adalah In[370]:=
M2 ‰- t d
Ha3 D3 + a2 D2 + a1 D + a0L
;
N = a3 D3 + a2 D2 + a1 D + a0 ê. D Ø -d;
Solusi Umum l2= C F + P I ; l2 = A1 ‰b 1 t + A2 ‰b 2 t + A3 ‰b 3 t +
M2 ‰- t d N
Menyelesaikanpersamaan l3 denganmenggunakanoperator - D; In[371]:=
In[372]:=
ClearAll@p, q, r, K, L, M, a, b, c, d, e, f, v, w, U, x, y, z, x4, y4, z4, d, ans1, ans2, ans3, p, h1, h2, B, l1, l2, l3, M1, M2, M3, N, ans1, ans2, ans3D;
Simplify@Solve@Eliminate@8slopecostate1, slopecostate2, slopecostate3<, 8l1, l2
Out[372]= 99l3 Ø
-I‰- t d M U4 ID2 h2KLw + H- bch1KL v - dh1Lp v + ach2KLw + h2p qwL x4
y4 + D Hdh1KL vy4 - h2w HLp x4 + Kqy4LLMM ë ID3 KLM + H bcf KLM + df LM p - beKM q - p qr + aKL HdeM - crLL x4y4z4 D2 HLM p x4 + KM qy4 + KLr z4L + D HacKLM x4y4 + M p qx4y4 + beKLM x4z4 + Lp r x4z4 - df KLM y4z4 + Kqr y4z4LM==
Misalkan l3= In[373]:=
A3 B3
, dengan ;
A3 = -I‰- t d M U4 2 ID h2KLw + H- bch1KL v - dh1Lp v + ach2KLw + h2p qwL x4y4 + D Hdh1KL vy4 - h2w HLp x4 + Kqy4LLMM;
46
In[374]:=
B3 = 3 ID KLM + H bcf KLM + df LM p - beKM q - p qr + aKL HdeM - crLL x4y4z4 D2 HLM p x4 + KM qy4 + KLr z4L + D HacKLM x4y4 + M p qx4y4 + beKLM x4z4 + Lp r x4z4 df KLM y4z4 + Kqr y4z4LM;
Untuk penyederhanaan, kalikan A & B dengan Mengalikan A3 dengan
In[375]:=
Out[375]=
A3 = ExpandBA3 *
1 K LM 1
KLM
-D2 ‰- t d h2U4w +
1 K LM
;
;
F
D ‰- t d h2p U4wx4
dD ‰- t d h1U4 vy4 +
K D ‰- t d h2qU4wy4
L
d ‰- t d h1p U4 vx4y4
+ bc ‰- t d h1U4 vx4y4 +
- ac ‰- t d h2U4wx4y4 -
K
‰- t d h2p qU4wx4y4
KL
"Mengevaluasi operator-D diperoleh:"; In[376]:=
A3 = A3 ê. D Ø -d;
In[377]:=
M3 = ExpandBFullSimplifyB
Out[377]=
bch1U4 vx4y4 + h2p U4wx4 d K
In[378]:=
A3 ‰- t d
FF
dh1p U4 vx4y4 K
+ dh1U4 vy4 d -
- ach2U4wx4y4 -
h2qU4wy4 d L
h2p qU4wx4y4 KL
-
- h2U4w d2
C3 = M3;
"Sehingga"; A3 = M3 ‰- t d ; Mengalikan B3 dengan
In[379]:=
Out[379]=
B3 = ExpandBB3 * D3 -
D2 p x4 K
bDex4z4 +
-
1 K LM 1
KLM
D2 qy4
L Dp r x4z4
;
F
+ acDx4y4 +
Dp qx4y4 KL Dqr y4z4
-
D2 r z4 M
+
- dDf y4z4 + + adex4y4z4 + KM LM df p x4y4z4 beqx4y4z4 acr x4y4z4 p qr x4y4z4 bcf x4y4z4 + K L M KLM
47
"Daftar koefisien persamaan B3 terhadap operator-D"; In[380]:= Out[380]=
Coeff3 = CoefficientList@B3, DD
:a d e x4 y4 z4 + b c f x4 y4 z4 +
b e q x4 y4 z4
-
d f p x4 y4 z4
a c r x4 y4 z4
-
K p q r x4 y4 z4
, KL M p q x4 y4 p r x4 z4 q r y4 z4 + b e x4 z4 + - d f y4 z4 + , a c x4 y4 + KL KM LM p x4 q y4 r z4 , 1> K L M L
M
"Misalkan,"; In[381]:= Out[381]=
a0 = Coeff3@@1DD a d e x4 y4 z4 + b c f x4 y4 z4 + b e q x4 y4 z4 L
In[382]:= In[383]:= Out[383]=
In[384]:= In[385]:= Out[385]=
In[386]:= In[387]:= Out[387]= In[388]:=
-
d f p x4 y4 z4
a c r x4 y4 z4 M
-
K p q r x4 y4 z4
KL M
A0 = a0;
a1 = Coeff3@@2DD a c x4 y4 +
p q x4 y4 KL
+ b e x4 z4 +
p r x4 z4 KM
- d f y4 z4 +
q r y4 z4 LM
A1 = a1;
a2 = Coeff3@@3DD -
p x4 K
-
q y4 L
-
r z4 M
A2 = a2;
a3 = Coeff3@@4DD 1 A3 = a3;
"Persamaan diferensial nya dapat ditulis sebagai:" ; Ha3D +a2D +a1D+a0 Ll3= M3 ‰ 3
2
-t d
;
"Dengan menggunakan metode operator-D untuk kasus tersebut diperoleh:"; Auxiliary Equation adalah a3m3 +3 a2m2 +3a1m+a0 =0, dengan akar-akar biHi=1,2,3L ;
48
Complementary Function HC FL adalah A1‰ b1 t+A2 b t b t ‰ 2 +A3‰ 3 , dengan Ai Hi=1,2,3Ladalah arbitrary constants ;
Particular Integral HPIL adalah
Ha3 D3+a2 D2+a1 D+a0L
M2 −t δ
;
Solusi Umum l3= C F + P I ; l3 = A1 ‰ b1 t + A2 ‰ b2 t + A3 ‰ b3 t + In[389]:= In[390]:=
M3 ‰-t d
ClearAll@M1, M2, M3, ND
N
H = Exp@- t d D Hh1 v y + h2 w z - BL U + l1 Jl2 Jl3 J-
p K
y2 + c x y - d y z + y Hq - v ULN +
q L r M
‰-t d M1
z2 + e x z - f y z + z Hr - w ULN;
In[391]:=
l1 =
In[392]:=
l2 =
In[393]:=
l3 =
In[394]:=
ans2 = Flatten@Solve@v y l2 + w z l3 == 0, yDD
Out[394]=
In[395]:=
N ‰-t d M2
N ‰-t d M3
N
:y Ø -
; ; ;
M3 w z M2 v
>
L21 = FullSimplify@-∑x HD ê. ans2 ‰-t d I2 M1 p x - K Ie M3 z -
Out[395]=
In[396]:=
In[397]:=
Out[398]=
+ M1 Ip - b z +
a M3 w z MMM M2 v
L22 = FullSimplify@-∑ y HD ê. ans2 LN
J‰-t d J-
2 M3 q w z v
+
L H-h1 N U v + a M1 x + f M3 z + M2 H- q + U v - c x + d zLLNN
L23 = FullSimplify@-∑z HD ê. ans2
‰- t d J2M3r z + M J- h2 NUw + bM1x - d M3 w z + M3 J-r + Uw - ex - f M3 w z NNN v
Out[397]=
In[398]:=
c M3 w z v
KN
1 Out[396]=
x2 - a x y - b x z + p xN +
M2 v
MN Tl1 = ∑t l1 -
‰- t d M1 d
N
49
In[399]:= Out[399]=
In[400]:= Out[400]=
Tl2 = ∑t l2 -
‰- t d M2 d
N
Tl3 = ∑t l3 -
‰- t d M3 d
N
In[401]:=
ansU = Flatten@FullSimplify@Solve@Eliminate@8Tl1 == L21, Tl2 == L22, Tl3 == L23<, 8x, z
In[402]:=
ans1 = Flatten@FullSimplify@Solve@Eliminate@8Tl1 == L21, Tl2 == L22, Tl3 == L23<, 8U, z
In[403]:=
ans3 = Flatten@FullSimplify@Solve@Eliminate@8Tl1 == L21, Tl2 == L22, Tl3 == L23<, 8U, x
In[404]:=
p = 2.5;
In[405]:=
q = 2;
In[406]:=
r = 2.5;
In[407]:=
K = 350000;
In[408]:=
L = 220000;
In[409]:=
M = 225000;
In[410]:=
a = 2.5 * 10-6 ;
In[411]:=
b = 1.2 * 10-6 ;
In[412]:=
c = 2.5 * 10-6 ;
In[413]:=
d = 1.55 * 10-5 ;
50
In[414]:=
e = 1.2 * 10-6 ;
In[415]:=
f = 3.75 * 10-6 ;
In[416]:=
v = 0.000147;
In[417]:=
w = 0.000257;
In[418]:=
B = 2500;
In[419]:=
h1 = 700;
In[420]:=
h2 = 500; 4r q
In[421]:= Out[421]=
LM True cr
In[422]:= Out[422]=
Out[423]= In[424]:=
Out[424]= In[425]:=
Out[425]= In[426]:=
Out[426]= In[427]:=
Out[427]=
>M
de True eq
In[423]:=
> Hd + f L2
cf
>L
True U4 = H-dLM p HeK + rL v + L Hp HcK + qL r + bKM Heq - crLL v + M HcKL H-f p + arL + q H- aeKL + p HeK - f L + rLLL wL ë 2 2 I2L H beKM + p rL v - 2LM H bcK + aeK + Hd + f L pL vw + 2M HacKL + p qL w M 8159.37 x4 = HK Hr H-aLq + p q + aLU4 vL - dLM Hf p - ar + aU4wL + bM H-qr + f L Hq - U4 vL + qU4wLLL ê H-LM HadeK + bcf K + df pL + beKM q + HacKL + p qL rL 301407. y4 = HL Hp r HcK + q - U4 vL - dM p HeK + r - U4wL + bKM Heq - cr - eU4 v + cU4wLLL ê H-LM HadeK + bcf K + df pL + beKM q + HacKL + p qL rL 126253. z4 = HM HaeKL Hq - U4 vL + cKL Hf p - ar + aU4wL p HeKq - f Lq + qr + f LU4 v - qU4wLLL ê HM HadeKL + df Lp + bK Hcf L - eqLL - HacKL + p qL rL 26215.4
51
In[428]:=
d = 0.1;
In[429]:=
M1 = C1
Out[429]=
-60.6264
In[430]:=
M2 = C2
Out[430]=
-492.611
In[431]:=
M3 = C3
Out[431]=
600.723
In[432]:=
Unprotect@ ND;
In[433]:=
a0 = A0;
In[434]:=
a1 = A1;
In[435]:=
a2 = A2;
In[436]:=
a3 = A3;
In[437]:= Out[437]= In[438]:= Out[438]= In[439]:= Out[439]= In[440]:= Out[440]= In[441]:= Out[441]= In[442]:= Out[442]= In[443]:= Out[443]= In[444]:= Out[444]=
N = a3 D3 + a2 D2 + a1 D + a0 ê. D Ø -d
-0.716155
FullSimplify@Tl1 == L21D
‰-0.1 t H- 168000. + 1. x - 1.74214zL ã 0
FullSimplify@Tl2 == L22D
‰-0.1 t H732063. + 1. U + 0.844693x- 19.1456zL ã 0
FullSimplify@Tl3 == L23D
‰-0.1 t H- 5850.92 + 1. U - 0.0032207x + 0.00760269zL ã 0
x = x ê. ans1 254768. x4 301407.
z = z ê. ans3 49805.6 z4 26215.4
52
z4 In[445]:= Out[445]= In[446]:= Out[446]= In[447]:= Out[447]=
M
11.6513
y = y ê. ans2 106185. y4 126253. y4
In[448]:= Out[448]= In[449]:= Out[449]= In[450]:= Out[450]= In[451]:= Out[451]= In[452]:= Out[452]= In[453]:= Out[453]= In[454]:= Out[454]= In[455]:=
100
L
100
57.3876 T = vy4U4 + wz4U4 206404.
U = U ê. ansU 6292.8 T1 = vy U 98225.7 T2 = wz U 80548.1 Td = T1 + T2 178774. U4 8159.37 U
Out[455]=
6292.8
In[456]:=
U4 - U
Out[456]=
1866.58
In[457]:= In[458]:= Out[458]= In[459]:=
Unprotect@pD;
p = Hh1 vy + h2wz - BL U
9.33 μ 107
ClearAll@p, q, r, K, L, M, a, b, c, d, e, f, v, w, U, x, y, z, x4, y4, z4, d, ans1, ans2, ans3, p, h1, h2, BD
53
Lampiran 6 (Kurva MSY)
Hasil ini diperoleh dengan menggunakan software Mathematica 5.1. In[328]:=
p = 2.5;
In[329]:=
q = 2;
In[330]:=
r = 2.5;
In[331]:=
K = 350000;
In[332]:=
L = 220000;
In[333]:=
M = 225000;
In[334]:=
a = 2.5 * 10-6 ;
In[335]:=
b = 1.2 * 10-6 ;
In[336]:=
c = 2.5 * 10-6 ;
In[337]:=
d = 1.55 * 10-5 ;
In[338]:=
e = 1.2 * 10-6 ;
In[339]:=
f = 3.75 * 10-6 ;
In[340]:=
v = 0.000147;
In[341]:=
w = 0.000257;
In[342]:=
In[343]:=
In[344]:=
y4 = HL Hp r Hc K + q - U vL - d M p He K + r - U w L + b K M He q - c r - e U v + c U w LLL ê H-L M Ha d e K + b c f K + d f pL + b e K M q + Ha c K L + p qL rL; z4 = HM Ha e K L Hq - U vL + c K L Hf p - a r + a U w L p He K q - f L q + q r + f L U v - q U w LLL ê HM Ha d e K L + d f L p + b K Hc f L - e qLL - Ha c K L + p qL rL; T = v U y4 + w U z4;
54
Solusi Pemanenan Maksimum In[345]:=
sol1 = FullSimplify@Solve@∑U T ã 0, UDD
In[346]:=
T = T ê. sol1
In[347]:=
ClearAll@TD
Out[345]= 88U →
8159.37<<
Out[346]= 8206404.<
In[348]:=
Plot@v U y4 + w U z4, 8U, 0 , 16500<, AxesLabel Ø 8U, T<, TextStyle Ø 8FontFamily Ø "Times", FontSize Ø 10
200000 150000 100000 50000
2500 Out[348]=
5000
7500
10000 12500 15000
U
Graphics
55
Lampiran 7 (Dinamika Sistem)
Hasil ini diperoleh dengan menggunakan software Mathematica 5.1. In[349]:= Unprotect@x, y, z, p , N, p, q, r, K, L, M, a, b, c, d, e, f, v, w D ; In[350]:=
ClearAll@parmval, x, y, z, p , N, p, q, r, K, L, M, a, b, c, d, e, f, v, w D;
In[351]:=
sysid
In[352]:=
intreset;
In[353]:=
plotreset;
In[354]:=
setstate[{x,y,z}];
In[355]:=
setparm[{p,q,r,K,L,M,a,b,c,d,e,f,v,w}];
In[356]:=
Mathematica 5.1.0, DynPac 10.69, 5ê24ê2007
slopevec = :p x J1 q y J1 -
y L
x K
N- ax y -b x z,
N + c x y - d y z - v U y, r z J1 -
z M
N + e x z - f y z - w U z>;
In[357]:=
eqstates = findpolyeq;
In[358]:=
parmval = 82.5, 2, 2.5, 350000, 220000, 225000, 2.5 * 10-6 , 1.2 * 10-6, 2.5 * 10-6, 1.55 * 10-5 , 1.2 * 10-6, 3.75 * 10-6 , 0.000147, 0.000257<;
In[359]:=
In[360]:=
E4 = FullSimplify@eqstates@@5DDD;
U = H- d L M p He K + rL v + L Hp Hc K + qL r + b K M He q - c rLL v + M Hc K L H-f p + a rL + q H- a e K L + p He K - f L + rLLL w L ê 2 H2 L Hb e K M + p rL v - 2 L M Hb c K + a e K + Hd + f L pL v w + 2 M Ha c K L + p qL w 2L;
In[361]:=
E4 = eqstateval@E4 D
In[362]:=
eigval@E4D
In[363]:=
classify@E4D;
Out[361]= 8301407., 126253., 26215.4< Out[362]= 8−1.74304 + 0.286408
, − 1.74304 − 0.286408 , −0.105861<
strictly stable
56
In[364]:=
evalparm@slopevec D
In[365]:=
<< Graphics`ParametricPlot3D`
In[366]:=
solution = NDSolveB:
Out[364]= :2.5 J1 −
x N x − 2.5 × 10−6 x y − 1.2 × 10−6 x z, 350000 y −6 −1.19943 y + 2.5 × 10 x y + 2 J1 − N y − 0.0000155 y z, 220000 z −6 −6 −2.09696 z + 1.2 × 10 x z − 3.75 × 10 y z + 2.5 J1 − N z> 225000
y£@tD y@tD z£@tD z@tD
x£@tD x@tD
== 2.5` J1 -
x@tD 350000
N - 2.5*^-6 y@tD - 1.2`*^-6 z@tD,
== -1.20648 + 2.5*^-6 x@tD + 2 J1 -
y@tD
220000
N - 0.0000155` z@tD,
== -2.10313 + 1.2`*^-6 x@tD - 3.75*^-6 y@tD + 2.5` J1 -
x@0D == 255297, y@0D == 107165, z@0D == 48695>,
z@tD
225000
N ,
8x, y, z< , 8t, 0, 150< , Method Ø Automatic F;
ParametricPlot3D@Evaluate@8 x@tD, y@tD, z@tD< ê. solutionD, 8 t, 0, 150<, PlotPoints Ø 1000, BoxRatios Ø 81, 1, 1<, AxesLabel Ø 8 "x", "y", "z"<, PlotRange Ø AllD; y 120000 110000 100000 00000
45000 40000 z 35000 30000 25000 260000 280000 x
In[367]:=
300000
intreset;
In[368]:= In[369]:=
plotreset;
57
In[370]:=
8X, Y, Z< = 8 x@tD, y@tD, z@tD< ê. Flatten@solutionD
Out[370]= 8InterpolatingFunction@880., 150.<<, <>D@tD,
InterpolatingFunction@880., 150.<<, <>D@tD, InterpolatingFunction@880., 150.<<, <>D@tD<
In[371]:=
Plot@8 X<, 8 t, 0, 25< , PlotRange Ø 8 250000, 320000< , AxesLabel Ø 8 "t", "x"< , PlotStyle Ø 8 RGBColor@0, 0, 1D
Out[371]= In[372]:=
10
15
20
25
Graphics
Plot@8 Y<, 8 t, 0, 50< , PlotRange Ø 8 80000, 130000<, AxesLabel Ø 8 "t", "y"< , PlotStyle Ø 8 RGBColor@0, 0, 1D
Out[372]= In[373]:=
20
30
40
50
Graphics
Plot@8 Z<, 8 t, 0, 50< , PlotRange Ø 8 20000, 50000< , AxesLabel Ø 8 "t", "z"<, PlotStyle Ø 8RGBColor@0, 0, 1D
Out[373]=
20
30
40
50
Graphics
58
Lampiran 8 (Syarat Teorema 7)
à Misalkan f0 = ‰-t d U Hv y h1 + w z h2 L, f1 = In[386]:=
f0 = Exp@-d tD Hh1 v y + h2 w zL U;
In[387]:=
f1 = p x J1 -
In[388]:=
f2 = q y J1 -
In[389]:=
f3 = r z J1 -
dx dy dz . , f2 = , f3 = dt dt dt
x N - a x y - b x z; K
y N + c x y - d y z - v U y; L
z N + e x z - f y z - w U z; M
dengan (x,y,z) adalah variabel state (autonomous) , t adalah waktu, dan sisanya adalah parameter. In[390]:= ∑t f0
Out[390]= - I‰-t d H- d L M p He K + rL v + L Hp Hc K + qL r + b K M He q - c rLL v +
M Hc K L H- f p + a rL + q H- a e K L + p He K - f L + rLLL wL d Hv y h1 + w z h2 LM ë I2 L Hb e K M + p rL v2 2 L M Hb c K + a e K + Hd + f L pL v w + 2 M Ha c K L + p qL w2M
Turunan parsial f1 terhadap x In[391]:= ∑x f1 Out[391]= -
px x + p J1 - N - a y - b z K K
Turunan parsial f2 terhadap y In[392]:= ∑y f2
Out[392]= - Hv H- d L M p He K + rL v + L Hp Hc K + qL r + b K M He q - c rLL v +
M Hc K L H- f p + a rL + q H- a e K L + p He K - f L + rLLL wLL ë
I2 L Hb e K M + p rL v2 - 2 L M Hb c K + a e K + Hd + f L pL v w + qy y 2 M Ha c K L + p qL w2M + c x + q J1 - N - d z L L
59
Turunan parsial f3 terhadap z In[393]:= ∑z f3
Out[393]= - Hw H- d L M p He K + rL v + L Hp Hc K + qL r + b K M He q - c rLL v +
M Hc K L H- f p + a rL + q H- a e K L + p He K - f L + rLLL wLL ë
I2 L Hb e K M + p rL v2 - 2 L M Hb c K + a e K + Hd + f L pL v w + rz z 2 M Ha c K L + p qL w2M + e x - f y + r J1 N M M
Terlihat bahwa turunan parsial dari f0, f1 , f2 , f3 kontinu In[382]:=
sehingga U dapat kontinu bagian demi bagian. H2.6 .3L
à ‡ f0 „ t ¶
0
Jika d > 0 maka Int = ‡ Exp@-d tD Hh1 v y + h2 w zL U „ t T
In[394]:=
0
Out[394]=
II1 - ‰-T dM H- d L M p He K + rL v + L Hp Hc K + qL r + b K M He q - c rLL v +
M Hc K L H- f p + a rL + q H- a e K L + p He K - f L + rLLL wL Hv y h1 + w z h2 LM ë II2 L Hb e K M + p rL v2 -
2 L M Hb c K + a e K + Hd + f L pL v w + 2 M Ha c K L + p qL w2M dM
In[395]:= Out[395]=
Limit@Int, T Ø ¶, Assumptions Ø d > 0D
HHc K L p r v + L p q r v - d L M p He K + rL v + b K L M He q - c rL v c f K L Mpw - a e K L Mqw +e K Mpqwf L M p q w + a c K L M r w + M p q r wL Hv y h1 + w z h2 LL ë I2 IM p q w2 + b K L M v He v - c wL +
L Ha K M w H- e v + c wL + p v Hr v - Hd + f L M wLLM dM
dari Definisi 18 akibatnya
integral pada (2.6.3) konvergen.
60
In[384]:=
Karena ‡ Exp@-d tD Hh1 v y + h2 w zL U „ t ada untuk d > 0, maka
In[384]:=
‡ Abs @Exp@-d Ht + sLD Hh1 v y + h2 w zL UD „ HtL
T
0
T
0
-Re@d H t +sLD ‡ ‰
T
Out[384]=
0
AbsAHH-d L M p He K + rL v + L Hp Hc K + qL r + b K M He q - c rLL v +
M Hc K L H-f p + a rL + q H- a e K L + p He K - f L + rLLL wL Hv y h1 + w z h2 LL ë I2 L Hb e K M + p rL v2 -
2 L M Hb c K + a e K + Hd + f L pL v w + 2 M Ha c K L + p qL w2 ME „ t
juga ada untuk d>0, s(-e,e),
In[385]:= In[396]:= Out[396]=
akibatnya untuk suatu fungsi yang kontinu bagian demi akibatnya bagian demi fungsi kontinu suatu untuk yang
In[386]:=
bagian a HtL, Abs @∑t Exp@-d Ht +sLD Hh1 v y + h2 w zL U D < a HtL
In[386]:=
untuk setiap s H-e, eL, t œ @0, ¶L , dan ‡ a HtL „ t < ¶, ¶
0
In[386]:=
sehingga
Teorema 7 dapat digunakan.
61
