Abstrak —Reaksi  Belousov  –  Zhabotinsky  merupakan  reaksi yang  berosilasi,  dimana  osilasi  yang  terjadi  pada  reaksi  ini memiliki  kemiripan  dengan  osilasi  yang  terjadi  pada  sistem biologi  seperti  gelombang  elektromagnetik  pada  jaringan jantung,  saraf  maupun  otot.  Reaksi  ini  memiliki  perilaku dinamik yang menarik untuk dikaji melalui analisis kestabilan dan bifurkasi. Selain untuk mengetahui perilaku dinamik suatu sistem,  analisis  mengenai  kestabilan  dan  bifurkasi  juga  dapat digunakan  untuk  mengetahui  faktor  apa  yang  mempengaruhi terjadinya  osilasi  pada  sistem  tersebut.  Analisis  kestabilan dapat dilakukan dengan menggunakan kriteria kestabilan nilai eigen dan Routh – Hurwitz, sementara untuk analisis bifurkasi dilakukan dengan menggunakan teorema bifurkasi Hopf. Dari hasil  analisis  didapatkan  bahwa  osilasi  terjadi  ketika  nilai koefisien stokiometri Cerium(IV) berada pada rentang 𝟏 𝟒 < 𝒇 < 𝟏 𝟐 .  Hasil  yang  didapatkan  ini  selanjutnya  disimulasikan menggunakan MATLAB.
Kata Kunci —Bifurkasi, Kestabilan, Oregonator.

## I. PENDAHULUAN
ELAMA  150  tahun  terakhir  keterkaitan  antara  ilmu biologi  dan  ilmu  kimia  menjadi  semakin  kuat,  dimulai dengan  munculnya  kimia  organik  dan  berlanjut  dengan biokimia.  Ilmu  kimia  tidak  hanya  sekedar  menyediakan pemahaman tentang dasar – dasar molekuler pada fenomena biologi, tetapi juga dapat menyediakan model dinamik untuk sistem biologi meskipun dasar yang ada dalam kedua sistem sangat berbeda. Hal ini didukung dengan teori kompleksitas yang  menyatakan bahwa  meskipun sistem  – sistem tersebut mempunyai  dasar  yang  berbeda,  seringkali  sistem  –  sistem tersebut  memiliki  karakteristik  dinamik  yang  sama  yang dapat  dianalisis  dan  dimodelkan  [1].  Salah  satu  contohnya adalah  reaksi  Belousov  –  Zhabotinsky  (B  –  Z)  yang  dapat digunakan untuk memodelkan sistem biologi.
Reaksi Belousov  – Zhabotinsky  merupakan  suatu reaksi osilasi  yang  dianggap  sebagai  purwarupa  osilator.  Osilasi yang  terjadi  pada  reaksi  Belousov  –  Zhabotinsky  memiliki kemiripan  dengan  osilasi  yang  terjadi  pada  sistem  biologis lain,  seperti  pada  gelombang  gerak  elektromagnetik  pada jaringan otot, jaringan saraf, dan agregasi jamur lendir [2].
Penelitian  mengenai  reaksi  osilator  telah  membuka kesempatan  bagi  para  peneliti  untuk  mengembangkan penerapan dari teori persamaan diferensial. Dalam penelitian yang  dilakukan  oleh  Field, Körös dan  Noyes,  didapatkan suatu model matematika berupa sistem persamaan diferensial yang dapat menggambarkan perilaku osilasi yang terjadi pada reaksi  Belousov  –  Zhabotinsky  [2].  Reaksi  kimia  dapat dianggap  sebagai  suatu  sistem  yang  perilaku  dinamiknya dapat  dikaji.  Dalam  reaksi  Belousov  –  Zhabotinsky  ini perilaku  dinamik  yang  menarik  untuk  dikaji  adalah bagaimana  perubahan  osilasi  yang  terjadi  seiring  dengan
perubahan  waktu  serta  bagaimana  perilaku  sistem  apabila nilai  parameter  tertentu  berubah.  Osilasi  yang  terjadi mungkin  akan  menghilang  atau  membesar  seiring  dengan berjalannya  waktu,  dimana  dalam  teori  matematika  sistem perilaku  ini  disebut  sebagai  kestabilan  sistem.  Perubahan perilaku  sistem  akibat  perubahan  nilai  parameter  disebut sebagai bifurkasi, dimana masalah ini memiliki kaitan yang erat  dengan  kestabilan  sistem.  Selain  untuk  mengetahui perilaku dinamik suatu sistem, analisis mengenai kestabilan dan bifurkasi juga dapat digunakan untuk mengetahui faktor apa  yang  mempengaruhi  terjadinya  osilasi  pada  sistem tersebut.
Penelitian  ini  akan  dilakukan  dalam  beberapa  tahap. Pertama  akan  dilakukan  studi  literatur,  dilanjutkan  dengan penentuan  titik  setimbang,  pelinieran  sistem,  uji  kestabilan menggunakan  kriteria  nilai  eigen  dan  Routh  –  Hurwitz. Setelah  kestabilan  sistem  didapatkan,  akan  dilanjutkan dengan  analisis  bifurkasi  yang  terdiri  dari  penentuan parameter  bifurkasi,  pembuktian  bahwa  terjadi  bifurkasi Hopf, serta penentuan rentang nilai parameter bifurkasi yang menyebabkan  sistem  tidak  stabil.  Hasil  yang  didapat selanjutnya  akan  disimulasikan  menggunakan  MATLAB sehingga  bisa  diambil  kesimpulan  dari  hasil  analisis  serta simulasinya.

## II. URAIAN PENELITIAN

## A. Kinetika Kimia
Laju dari reaksi 𝑎𝐴 + 𝑏𝐵 → 𝑐𝐶 + 𝑑𝐷 diberikan oleh
𝑙𝑎𝑗𝑢 = − 1 𝑎 ∆[𝐴] ∆𝑡 = − 1 𝑏 ∆[𝐵] ∆𝑡 = 1 𝑐 ∆[𝐶] ∆𝑡 = 1 𝑑 ∆[𝐷] ∆𝑡
Dengan hukum laju, didapatkan
𝑙𝑎𝑗𝑢 = 𝑘[𝐴] 𝑥 [𝐵] 𝑦
Sehingga laju perubahan konsentrasi menjadi ∆[𝐴] ∆𝑡 = −𝑎𝑘[𝐴] 𝑥 [𝐵] 𝑦 ∆[𝐵] ∆𝑡 = −𝑏𝑘[𝐴] 𝑥 [𝐵] 𝑦 ∆[𝐶] ∆𝑡 = 𝑐𝑘[𝐴] 𝑥 [𝐵] 𝑦 ∆[𝐷] ∆𝑡 = 𝑑𝑘[𝐴] 𝑥 [𝐵] 𝑦
Laju  perubahan  konsentrasi  tersebut  dapat  ditulis  sebagai suatu sistem persamaan diferensial
𝐴̇ = −𝑎𝑘[𝐴] 𝑥 [𝐵] 𝑦 𝐵̇ = −𝑏𝑘[𝐴] 𝑥 [𝐵] 𝑦
𝐶̇ = 𝑐𝑘[𝐴] 𝑥 [𝐵] 𝑦
𝐷̇ = 𝑑𝑘[𝐴] 𝑥 [𝐵] 𝑦 Dengan konsentrasi awal 𝐴(0) =   𝐴 0 , 𝐵(0) =   𝐵 0 , 𝐶(0) =   𝐶 0 , 𝐷(0) =   𝐷 0 [3].

## Analisis Kestabilan dan Bifurkasi Model Oregonator dari Reaksi Belousov – Zhabotinsky
Amelia M. S., Erna Apriliani, dan I Gst Ngr Rai Usadha Departemen Matematika, Fakultas Matematika, Komputasi dan Sains Data,
Institut Teknologi Sepuluh Nopember (ITS)
e-mail : april@matematika.its.ac.id
S

## B. Model Oregonator
Menurut  Leigh  seperti  dikutip  [4],  pada  tahun  1970 beberapa peneliti dari Universitas Oregon mulai mengembangkan  suatu  model  yang  dapat  menjelaskan perilaku eksperimental dari reaksi Belousov  – Zhabotinsky. Model  yang  dikenal  sebagai  model  FKN  disederhanakan menjadi model Oregonator.
Reaksi kimia yang digunakan untuk mendapatkan model Oregonator adalah sebagai berikut.
(𝑅3) 𝐵𝑟𝑂 3 − + 𝐵𝑟 − + 2𝐻 + → 𝐻𝐵𝑟𝑂 2 + 𝐻𝑂𝐵𝑟 (𝑅2) 𝐻𝐵𝑟𝑂 2 + 𝐵𝑟 − + 𝐻 + → 2𝐻𝑂𝐵𝑟 (𝐺) 2𝐶𝑒 3+ + 𝐵𝑟𝑂 3 − + 𝐻𝐵𝑟𝑂 2 + 3𝐻 + → 2𝐶𝑒 4+ + 2𝐻𝐵𝑟𝑂 2 + 𝐻 2 𝑂 (𝑅4) 2𝐻𝐵𝑟𝑂 2 → 𝐵𝑟𝑂 3 − + 𝐻𝑂𝐵𝑟 + 𝐻 + (𝑅10) 4𝐶𝑒 4+ + 𝐵𝑟𝑀𝐴 + 2𝐻 2 𝑂 → 4𝐶𝑒 3+ + 𝐻𝐶𝑂𝑂𝐻 + 2𝐶𝑂 2 + 5𝐻 +
konsentrasi 𝐵𝑟𝑂 3 − , 𝐵𝑟𝑀𝐴 dan 𝐻𝑂𝐵𝑟 adalah  sama  karena ketiga zat ini  masuk dan  keluar dari medium reaksi dengan laju konstan. Sebaliknya, konsentrasi 𝐻𝐵𝑟𝑂 2 , 𝐵𝑟 − dan 𝐶𝑒 4+ berubah seiring berjalannya waktu [5]. Selanjutnya, misalkan X = [ 𝐻𝐵𝑟𝑂 2 ] Y = [ 𝐵𝑟 − ] Z = [ 𝐶𝑒 4+ ] A = [ 𝐵𝑟𝑂 3 − ] B = [𝐵𝑟𝑀𝐴]
P = [ 𝐻𝑂𝐵𝑟 ]
Maka reaksinya menjadi (𝑀1) 𝐴 + 𝑌 → 𝑋 + 𝑃 (𝑀2) 𝑋 + 𝑌 → 2𝑃 (𝑀3) 𝐴 + 𝑋 → 2𝑋 + 2𝑍 (𝑀4) 2𝑋 → 𝐴 + 𝑃 (𝑀5)𝐵 +  𝑍 → 𝑓𝑌 Dengan laju 𝑟 (𝑀1) = 𝑘 𝑅3 [𝐻 + ] 2 𝐴𝑌 𝑟 (𝑀2) = 𝑘 𝑅2 [𝐻 + ]𝑋𝑌 𝑟 (𝑀3) = 𝑘 𝑅5 [𝐻 + ]𝐴𝑋 𝑟 (𝑀4) = 𝑘 𝑅4 𝑋 2 𝑟 (𝑀5) = 𝑘 0 𝐵𝑍
Dengan  menggunakan  hukum  laju,  didapat  suatu  sistem persamaan diferensial tak linier sebagai berikut.
𝑑𝑋 𝑑𝑡 = 𝑘 3 𝐴𝑌 − 𝑘 2 𝑋𝑌 + 𝑘 5 𝐴𝑋 − 2𝑘 4 𝑋 2 𝑑𝑌 𝑑𝑡 = −𝑘 3 𝐴𝑌 − 𝑘 2 𝑋𝑌 + 𝑓𝑘 0 𝐵𝑍 𝑑𝑍 𝑑𝑡 = 2𝑘 5 𝐴𝑋 − 𝑘 0 𝐵𝑍
Dimana: X = konsentrasi asam Bromat / [𝐻𝐵𝑟𝑂 2 ] Y = konsentrasi ion Bromida / [𝐵𝑟 − ] Z = konsentrasi ion Ce(IV) / [𝐶𝑒 4+ ] A = konsentrasi ion Bromat / [𝐵𝑟𝑂 3 − ] B = konsentrasi asam Bromo Malonat / [𝐵𝑟𝑀𝐴] P = konsentrasi asam Hipobromat / [𝐻𝑂𝐵𝑟] 𝑓 = konstanta stoikiometri
Berdasarkan [6], apabila digunakan suatu larutan dengan [𝐻 + ]~1 𝑀 pada 0.8 𝑀 𝐻 2 𝑆𝑂 4 dengan  mengasumsikan konsentrasi awal 𝐴 0 = [𝐵𝑟𝑂 3− ] 0 = 0.06 𝑀 𝐵 0 = [𝐵𝑟𝑀𝐴] 0 = 10 −3 𝑀 maka didapat nilai parameter yang ditunjukkan pada Tabel 1:
Tabel 1. Parameter Model Oregonator [6] Parameter Deskripsi Nilai 𝑘 0 Konstanta laju dari reaksi 5 0.4 𝑀 −1 𝑆 −1 𝑘 2 konstanta laju dari reaksi 2 2 × 10 9 𝑀 −1 𝑆 −1 𝑘 3 konstanta laju dari reaksi 1 2 𝑀 −1 𝑆 −1 𝑘 4 konstanta laju dari reaksi 4 4 × 10 7 𝑀 −1 𝑆 −1 𝑘 5 konstanta laju dari reaksi 3 10 4 𝑀 −1 𝑆 −1

## C. Kriteria Nilai Eigen
Teorema 2.5.3.1.1 [7] Diberikan  persamaan  diferensial ẋ = Ax dengan  matriks  A berukuran n × n mempunyai nilai karakteristik yang berbeda λ 1 , … , λ k (k ≤ n) . 1. Titik  asal 𝑥̅ = 0 adalah  stabil  asimtotik  bila  dan  hanya bila 𝑅𝑒(𝜆 𝑖 ) < 0 untuk semua 𝑖 = 1, 2, … , 𝑘 . 2. Titik  asal 𝑥̅ = 0 adalah  stabil  bila  dan  hanya  bila 𝑅𝑒(𝜆 𝑖 ) ≤ 0 untuk semua 𝑖 = 1, 2, … , 𝑘 , dan untuk semua 𝜆 𝑖 dengan 𝑅𝑒(𝜆 𝑖 ) = 0 multiplisitas aljabar sama dengan multisiplitas geometrinya.
3. Titik asal 𝑥̅ = 0 adalah tidak stabil bila dan hanya bila
𝑅𝑒(𝜆 𝑖 ) > 0 untuk  beberapa 𝑖 = 1, 2, … , 𝑘 ,  atau  ada 𝜆 𝑖 dengan 𝑅𝑒(𝜆 𝑖 ) = 0 dengan  multiplisitas  aljabar  lebih besar dari multiplisitas geometrinya.
D. Aturan Tanda Descartes
Diberikan suatu polinomial karakteristik 𝑝(𝜆) = 𝑎 0 𝜆 𝑛 + 𝑎 1 𝜆 𝑛−1 + ⋯ + 𝑎 𝑛−1 𝜆 + 𝑎 𝑛 = 0 Tanpa  mengurangi  keumuman,  diasumsikan 𝑎 𝑛 > 0 . Misalkan 𝑁 adalah banyaknya perubahan tanda pada barisan koefisien {𝑎 𝑛 , 𝑎 𝑛−1 , … , 𝑎 0 } tanpa  menghiraukan  koefisien yang bernilai 0.  Aturan tanda Descartes  menyatakan bahwa 𝑝(𝜆) mempunyai paling banyak 𝑁 akar – akar yang real dan positif, dan lagi, ada 𝑁, 𝑁 − 2 atau 𝑁 − 4, … akar real positif [8].

## E. Kriteria Routh – Hurwitz
Diberikan suatu polinomial 𝑞(𝑠) = 𝑎 𝑛 𝑠 𝑛 + 𝑎 𝑛−1 𝑠 𝑛−1 + ⋯ + 𝑎 1 𝑠 + 𝑎 0 , 𝑎 𝑛 ≠ 0 Dari  polinomial  tersebut,  dapat  disusun  suatu  Tabel  seperti yang dapat dilihat pada Tabel 2:
Tabel 2.
Tabel Routh Hurwitz 𝑠 𝑛 𝑎 𝑛 𝑎 𝑛−2 𝑎 𝑛−4 … 𝑠 𝑛−1 𝑎 𝑛−1 𝑎 𝑛−3 𝑎 𝑛−5 … 𝑠 𝑛−2 𝑏 1 𝑏 2 𝑏 3 𝑠 𝑛−3 𝑐 1 𝑐 2 𝑐 3 ⋮ ⋮ 𝑠 0
𝑞
dimana 𝑏 1 , 𝑏 2 , … , 𝑐 1 , 𝑐 2 , … dan 𝑞 secara  rekursif  didapatkan dari:
𝑏 1 = 𝑎 𝑛−1 𝑎 𝑛−2 − 𝑎 𝑛 𝑎 𝑛−3 𝑎 𝑛−1 , 𝑏 1 = 𝑎 𝑛−1 𝑎 𝑛−4 − 𝑎 𝑛 𝑎 𝑛−5 𝑎 𝑛−1 , … 𝑐 1 = 𝑏 1 𝑎 𝑛−3 − 𝑏 2 𝑎 𝑛−1 𝑏 1 , 𝑐 2 = 𝑏 1 𝑎 𝑛−5 − 𝑏 3 𝑎 𝑛−1 𝑏 1 , …
Kriteria  Routh  –  Hurwitz  menyimpulkan  bahwa  banyaknya perubahan tanda dalam kolom pertama pada tabel diatas sama dengan banyaknya akar – akar polinomial 𝑞(𝑠) yang bagian realnya positif [8].
F. Orbit Periodik
Untuk suatu medan vector 𝑥̇ = 𝑓(𝑥, 𝑦)
𝑦̇ = 𝑔(𝑥, 𝑦)
Dimana 𝑓, 𝑔  ∈ 𝐶 1 dengan 𝐶 1 merupakan  himpunan fungsi kontinu yang turunan pertamanya juga kontinu, berlaku teorema berikut.
Teorema 2.5.4.1 [9]
Jika  pada  suatu  daerah  terhubung  sederhana 𝐷 ⊂ ℝ 2
(misalnya D tidak mempunyai lubang) 𝜕𝑓 𝜕𝑥 + 𝜕𝑔 𝜕𝑦 tidak nol dan tidak  berubah  tanda,  maka  system  tidak  memiliki  orbit tertutup yang terletak seluruhnya di D.
G. Bifurkasi
Bifurkasi  adalah  perubahan  kualitatif  dalam  struktur asimtotik suatu sistem dinamik.
Definisi 2.7.1 [9] Suatu  titik  tetap (𝑥, 𝜇) = (0, 0) dari  suatu  keluarga  vector fields  satu  dimensi  dengan  satu  parameter  dikatakan mengalami bifurkasi pada 𝜇 = 0 jika flow untuk 𝜇 disekitar 0 dan 𝑥 disekitar 0 tidak sama secara kualitatif dengan flow disekitar 𝑥 = 0 pada 𝜇 = 0 .
Bifurkasi Hopf adalah bifurkasi periodik dimana suatu limit cycle baru  muncul  dari  suatu  solusi  stasioner  ketika  suatu parameter  melewati  nilai  kritis.  Ketika limit  cycle muncul, kestabilan dari solusi stasioner turut berubah. Bifurkasi Hopf hanya bisa terjadi pada dimensi 2 atau lebih.
Definisi 2.7.2 [9] Misalkan 𝑥 0 adalah  titik  kesetimbangan  regular  untuk 𝜇 = 𝜇 0 . 𝑥 0 disebut sebagai titik bifurkasi Hopf atau disingkat titik Hopf untuk 𝜇 =   𝜇 0 jika dan hanya jika dua kondisi berikut terpenuhi:
1. Kondisi nilai eigen
J mempunyai sepasang nilai eigen ±𝑖𝜔 0 dimana 𝜔 ≠ 0 dan  untuk 𝜇 = 0 ,  tidak  ada  nilai  eigen  dari J yang merupakan nilai eigen regular.
2. Kondisi transversal
𝛼 ′ (𝜇 0 ) ≠ 0 (Dua  kurva  nilai  eigen 𝛼(𝜇) ± 𝑖𝛽(𝜇) melintasi  sumbu imajiner dengan kecepatan tak nol)
Pada bifurkasi Hopf, nilai – nilai eigen melewati sumbu imajiner.  Titik  merah  pada  gambar  2.4  mengindikasikan bifurkasi  Hopf.  Jika limit  cycle yang  muncul  stabil,  maka bifurkasi  Hopf  disebut supercritical ,  sedangkan  jika limit cycle yang muncul tidak stabil, maka bifurkasi Hopf disebut subcritical .  Periode  dari  limit  cycle  yang  muncul  diberikan oleh:
𝑇 = 2𝜋 |𝐼𝑚(𝜆(𝜇 0 ))|

## III. ANALISIS DAN PEMBAHASAN

## A. Model Matematika Oregonator
Dalam [6] diambil parameter tanpa dimensional: 𝑥 = 2𝑘4 𝑘5𝐴 𝑋
𝑦 = 𝑘2 𝑘5𝐴 𝑌 𝑧 = 𝑘0𝑘4𝐵 (𝑘5𝐴)2 𝑍 𝜏 = 𝑘 0 𝐵𝑡 𝜀 = 𝑘0 𝑘5𝐴 𝐵 𝛿 = 2𝑘0𝑘4 𝑘2𝑘5𝐴 𝐵
𝑞 = 2𝑘3𝑘4 𝑘2𝑘5
Maka didapat suatu system persamaan baru sebagai berikut:
𝑑𝑥 𝑑𝜏 = 1 𝜀 (𝑞𝑦 − 𝑥𝑦 + 𝑥(1 − 𝑥))
𝑑𝑦 𝑑𝜏 = 1 𝛿 (−𝑞𝑦 − 𝑥𝑦 + 2𝑓𝑧) (3.1) 𝑑𝑧 𝑑𝜏 = 𝑥 − 𝑧
Sistem  (3.1)  memiliki  dua  titik  setimbang,  yaitu 𝐸 1 ∗ = (0,0,0) dan 𝐸 2 ∗ = (𝑥 ∗ , 𝑦 ∗ , 𝑧 ∗ ) dimana
𝑥 ∗ = 1−𝑓−𝑞+√1+𝑞2+4𝑓2+2𝑞+12𝑞𝑓−4𝑓 2
𝑦 ∗ = 2𝑓(1−𝑓−𝑞+√1+𝑞2+4𝑓2+2𝑞+12𝑞𝑓−4𝑓)
1+𝑞−𝑓+√1+𝑞2+4𝑓2+2𝑞+12𝑞𝑓−4𝑓
𝑧 ∗ = 1−𝑓−𝑞+√1+𝑞2+4𝑓2+2𝑞+12𝑞𝑓−4𝑓 2 Dengan 1 − 2𝑓 − 𝑞 > 0 dan 1 + 𝑞 2 + 4𝑓 2 + 2𝑞 + 12𝑞𝑓 − 4𝑓 > 0 .

## B. Analisa Kestabilan
Sistem  dinamik  model  Oregonator  pada  (3.1)  merupakan sistem yang tidak linier, sehingga untuk melakukan analisis kestabilan  perlu  dilakukan  pelinieran  di  sekitar  titik setimbangnya  terlebih  dahulu  dengan  membentuk  matriks Jacobian.  Matriks  Jacobian  dari  sistem  (3.1)  di  sekitar  titik setimbang 𝐸 1 ∗ = (0,0,0) adalah:
𝑨 𝟏 = ( 1 𝜀 𝑞 𝜀 0 0 − 𝑞 𝛿 2𝑓 𝛿 1 0 −1)
Persamaan karakteristik dari matriks 𝐴 1 diberikan oleh 𝜆 3 + ( 𝜀(𝛿+𝑞)−𝛿 𝜀𝛿 ) 𝜆 2 + ( 𝑞(𝜀+1)−𝛿 𝜀𝛿 ) 𝜆 − ( 𝑞(2𝑓+1) 𝜀𝛿 ) = 0
Berdasarkan aturan tanda Descartes dan kriteria nilai eigen, titik setimbang 𝐸 1 ∗ = (0,0,0) tidak stabil.
Matriks  Jacobian  dari  sistem  (3.1)  di  sekitar  titik setimbang 𝐸 2 ∗ = (𝑥 ∗ , 𝑦 ∗ , 𝑧 ∗ ) adalah:
𝑨 𝟐 = ( −𝑦 ∗ + 1 − 2𝑥 ∗ 𝜀 𝑞 − 𝑥 ∗ 𝜀 0 −𝑦 ∗ 𝛿 − 𝑞 + 𝑥 ∗ 𝛿 2𝑓 𝛿 1 0 −1)
Persamaan karakteristik dari matriks 𝐴 2 diberikan oleh 𝜆 3 + 𝐴𝜆 2 + 𝐵𝜆 + 𝐶 = 0 Dimana:
𝐴 = 𝑥∗ 2 +𝑞(𝑥∗+2𝑓) 𝑞+𝑥∗ 𝜀 + 𝑞+𝑥∗ 𝛿 + 1 > 0 𝐵 = 𝑥∗ 2 +𝑞(𝑥∗+2𝑓) 𝑞+𝑥∗ 𝜀 + 𝑞+𝑥∗ 𝛿 + ( 𝑥∗ 2 +𝑞(𝑥∗+2𝑓)
𝑞+𝑥∗ )(𝑞+𝑥∗)+ 2𝑓𝑥∗ 𝑞+𝑥∗ (𝑞−𝑥∗)
𝜀𝛿
𝐶 = 𝑥∗ 2 +𝑞(1+2𝑓) 𝜀𝛿 > 0
Berdasarkan  kriteria  Routh  –  Hurwitz,  kestabilan  sistem dipengaruhi oleh 𝐵 . Apabila 𝐵 > 0 , sistem stabil dan supaya sistem  stabil  haruslah 𝐴𝐵 − 𝐶 > 0 .  Apabila 𝐵 < 0 ,  sistem tidak stabil dan supaya sistem stabil haruslah 𝐴𝐵 − 𝐶 < 0 . Titik  setimbang 𝐸 2 ∗ = (𝑥 ∗ , 𝑦 ∗ , 𝑧 ∗ ) mengalami  bifurkasi ketika 𝐴𝐵 − 𝐶 = 0 dengan  nilai  eigen 𝜆 1 = −𝐴, 𝜆 2,3 =
±𝑖√ 𝐶 𝐴 .

## C. Orbit Periodik
Asumsikan  perubahan  konsentrasi 𝑧 adalah  konstan, dimisalkan 𝑀 , dengan 𝑀 > 0 , sehingga didapat: 𝑓(𝑥, 𝑦, 𝜏) = 𝑑𝑥 𝑑𝜏 = 1 𝜀 (𝑞𝑦 − 𝑥𝑦 + 𝑥(1 − 𝑥))
𝑔(𝑥, 𝑦, 𝜏) = 𝑑𝑦 𝑑𝜏 = 1 𝛿 (−𝑞𝑦 − 𝑥𝑦 + 𝑓𝑥 − 𝑓𝑀)
(3.2)
Dari sistem (3.2) didapat;
𝜕𝑓 𝜕𝑥 + 𝜕𝑔 𝜕𝑦 = ( 1 𝜀 (−𝑦 + 1 − 2𝑥)) + ( 1 𝛿 (−𝑞 − 𝑥)) 𝜕𝑓 𝜕𝑥 + 𝜕𝑔 𝜕𝑦 = −𝑦𝛿−𝑥(2𝛿+𝜀)+𝛿−𝑞𝜀 𝜀𝛿 Karena 𝑥, 𝑦, 𝜀, 𝛿, 𝑞 > 0 dan 𝛿 ≪ 𝜀 , maka 𝜕𝑓
𝜕𝑥 + 𝜕𝑔 𝜕𝑦 = −𝑦𝛿−𝑥(2𝛿+𝜀)+𝛿−𝑞𝜀 𝜀𝛿 < 0
Karena 𝜕𝑓 𝜕𝑥 + 𝜕𝑔 𝜕𝑦 ≠ 0 dan tidak berubah tanda, maka menurut Kriteria  Bendixson  system  tidak  mempunyai  orbit  tertutup yang berada seluruhnya di suatu daerah 𝐷 ⊂ ℝ 2 . Selanjutnya dengan mengasumsikan bahwa 𝜕𝑓 𝜕𝑥 + 𝜕𝑔 𝜕𝑦 = 0
𝜕𝑓 𝜕𝑥 + 𝜕𝑔 𝜕𝑦 = −𝑦𝛿−𝑥(2𝛿+𝜀)+𝛿−𝑞𝜀 𝜀𝛿
= 0

## Didapat
𝑥 = 𝛿(1−𝑦)−𝑞𝜀 2𝛿+𝜀 Kedua  garis  tersebut  membagi  bidang  menjadi  tiga  bidang terpisah, yaitu 𝑅 1 : 𝑥 < 𝑥 1 , 𝑅 2 :  𝑥 1 < 𝑥 < 𝑥 2 dan 𝑅 3 : 𝑥 > 𝑥 2 . Berdasarkan  kriteria  Bendixson,  sistem  (3.2)  tidak mempunyai  orbit  periodik  yang  terletak  sepenuhnya  pada 𝑅 1 , 𝑅 2 atau 𝑅 3 ,  tapi  tidak  menutup  kemungkinan  bahwa sistem  (3.2)  mempunyai  orbit  periodic  yang  terletak  di gabungan dari dua atau lebih bidang – bidang tersebut ( 𝑅 1 , 𝑅 2 dan 𝑅 3 ).

## D. Analisis Bifurkasi
Ketika 𝐴𝐵 − 𝐶 = 0 , didapat nilai eigen
𝜆 1 = − 𝑥∗ 2 +𝑞(𝑥∗+2𝑓)
𝑞+𝑥∗ 𝜀 − 𝑞+𝑥∗ 𝛿 − 1 𝜆 2,3 = ±𝑖√ (𝑞+𝑥∗)(𝑞+𝑥∗2)+2𝑓(𝑞+𝑥∗)
2𝑞𝛿𝑓+𝛿𝑥(𝑞+𝑥∗)+𝜀(𝑞+𝑥∗)(𝑞+𝑥∗+1)
Dapat dilihat bahwa 𝐴 dan 𝐶 adalah fungsi dari 𝑞, 𝜀, 𝛿, 𝑓 dan 𝑥 ∗ dimana 𝑞, 𝜀, 𝑑 adalah  parameter  tanpa  dimensi  dan 𝑥 ∗ adalah  titik  setimbang  dan  merupakan  fungsi  dari 𝑓 dan 𝑞 . Nilai parameter tanpa dimensi 𝑞 , 𝜀, 𝛿 bergantung pada reaksi yang ada pada tabel 2.1 sedangkan nilai 𝑓 tidak bergantung pada laju reaksi ataupun reaktan, sehingga 𝑓 dapat dianggap sebagai parameter bebas [2]. Oleh karena itu, akan diambil 𝑓 sebagai parameter bifurkasi.
Karena 𝑨 𝟐 mempunyai  sepasang  nilai  eigen  dengan bagian  real  0,  maka  kondisi  nilai  eigen  dari  definisi  Hopf terpenuhi.  Setelah  dilakukan  analisa  terhadap  turunan  total dari persamaan karakteristik 𝑨 𝟐 didapatkan 𝑅𝑒 ( 𝑑 𝑑𝑓 𝜆(𝑓 0 )) > 0 sehingga kondisi  transversalitas terpenuhi dan  limit cycle yang muncul memiliki periode
𝑇 = 2𝜋
√ (𝑞+𝑥∗)(𝑞+𝑥∗ 2 )+2𝑓(𝑞+𝑥∗) 2𝑞𝛿𝑓+𝛿𝑥(𝑞+𝑥∗)+𝜀(𝑞+𝑥∗)(𝑞+𝑥∗+1) .
Dari hasil analisa kestabilan, diketahui bahwa sistem mengalami bifurkasi ketika 𝐴𝐵 − 𝐶 = 0 𝐵 = 𝐶 𝐴 > 0
Sesuai dengan hasil analisa tabel Routh – Hurwitz, ketika 𝐵 > 0 kondisi yang harus dipenuhi supaya sistem tidak stabil diberikan oleh:
𝐴𝐵 − 𝐶 < 0 𝑁𝛿2+𝑀𝛿+𝐿 𝛿2 < 0 Dengan
𝑁 = (𝑥∗ 2 +𝑞(𝑥∗+2𝑓))(1+ 𝑥∗ 2 +𝑞(𝑥∗+2𝑓) 𝜀(𝑞+𝑥∗) )
𝜀(𝑞+𝑥∗) ,
𝑀 =
𝑥∗(2𝑥∗+3𝑞)+𝑞(4𝑓−1) 𝜀 + (𝑥∗ 2 +𝑞(𝑥∗+2𝑓)) 2
𝜀2(𝑞+𝑥∗) + 2𝑓𝑥∗(𝑞−𝑥∗)(𝑥∗ 2 +𝑞(𝑥∗+2𝑓)) 𝜀2(𝑞+𝑥∗)2 + 2𝑓𝑥∗(𝑞−𝑥∗) 𝜀(𝑞+𝑥∗) +(𝑞+𝑥∗) 𝛿 ,
𝐿 = (𝑞 + 𝑥 ∗ ) ((𝑞 + 𝑥 ∗ ) + 𝑥∗(1−𝑞−4𝑓)+2𝑞(1+3𝑓)
𝜀 ) . sehingga
𝑁𝛿 2 + 𝑀𝛿 + 𝐿 < 0
(𝛿 − ( −𝑀+√𝑀2−4𝑁𝐿 2𝑁 )) (𝛿 − ( −𝑀−√𝑀2−4𝑁𝐿 2𝑁 )) < 0
Supaya 𝛿 bernilai real, maka kondisi berikut harus terpenuhi:
𝑀 2 − 4𝑁𝐿 > 0 Karena 𝑀 2 > 0 dan 𝑁 > 0 , maka haruslah: 𝐿 < 0
(𝑞 + 𝑥 ∗ ) ((𝑞 + 𝑥 ∗ ) + 𝑥∗(1−𝑞−4𝑓)+2𝑞(1+3𝑓) 𝜀 ) < 0
(1 − 2𝑓 − 𝑞 + √(1 − 𝑞 − 2𝑓) 2 + 4𝑞(2𝑓 + 1))(1 − 𝑞 − 4𝑓) + 4𝑞(1 + 3𝑓) < 0
(3.3) Untuk suatu nilai 𝑞 yang sangat kecil, didapat
(1 − 2𝑓)(1 − 4𝑓) < 0
Penyelesaian  dari  (3.3)  juga  merupakan  rentang  nilai 𝑓 supaya sistem tidak stabil diberikan oleh:
1 4 < 𝑓 < 1 2

## IV. SIMULASI DAN ANALISIS
Simulasi  ini  akan  menunjukkan  pengaruh  perubahan  nilai parameter 𝑓 terhadap osilasi yang terjadi pada reaksi B-Z. Nilai parameter yang digunakan adalah
𝜀 = 5 × 10 −5 𝛿 = 2 × 10 −4 𝑞 = 8 × 10 −4
Dari [3] didapat: 𝑥 0 ≈ 1.27 × 10 −7 𝑀 𝑦 0 ≈ 6 × 10 −7 𝑀 𝑧 0 ≈ 5 × 10 −3 𝑀 1. 𝑓 = 0.1
Gambar 1. Simulasi untuk 𝑓 = 0.1 .
Karena 𝑓 < 0.25 ,  maka  reaksi  stabil.  Titik  setimbang sistem adalah 𝐸 2 ∗ = (0.798, 0.2, 0.798) , yang berarti reaksi
mencapai  kesetimbangan  ketika  konsentrasi 𝐻𝐵𝑟𝑂 2 dan 𝐶𝑒(𝐼𝑉) sebesar 0.798 𝑀 ,  dan  konsentrasi 𝐵𝑟 − sebesar 0.2 𝑀 . Pada awal reaksi, konsentrasi 𝐻𝐵𝑟𝑂 2 mengalami kenaikan yang  paling  tinggi,  menandakan  bahwa 𝐻𝐵𝑟𝑂 2 banyak dihasilkan.  Setelah  itu,  konsentrasi 𝐻𝐵𝑟𝑂 2 menurun  secara drastis  sedangkan  konsentrasi 𝐵𝑟 − naik  secara  drastis  dan 𝐶𝑒(𝐼𝑉) naik sedikit, menandakan 𝐻𝐵𝑟𝑂 2 banyak digunakan sebagai  reaktan  sementara 𝐵𝑟 − banyak  dihasilkan  dan 𝐶𝑒(𝐼𝑉) dihasilkan sedikit. Setelah itu, konsentrasi 𝐵𝑟 − tetap naik sedangkan konsentrasi 𝐶𝑒(𝐼𝑉) mulai turun, menandakan 𝐵𝑟 − dihasilkan  sedikit  dan 𝐶𝑒(𝐼𝑉) digunakan  sebagai reaktan.  Konsentrasi 𝐵𝑟 − mulai  turun  segera  setelah  satuan waktu ke 0.9. Seiring berjalannya waktu, konsentrasi 𝐵𝑟 − dan 𝐶𝑒(𝐼𝑉) terus turun, sedangkan konsentrasi 𝐻𝐵𝑟𝑂 2 mulai naik setelah satuan waktu keempat. Saat konsentrasi 𝐻𝐵𝑟𝑂 2 dan 𝐶𝑒(𝐼𝑉) sama  pada  satuan  waktu  keenam,  konsentrasi ketiganya  menjadi  konstan,  menandakan  reaksi  menjadi setimbang.
2. 𝑓 = 0.2 Karena 𝑓 < 0.25 ,  maka  reaksi  stabil.  Titik  setimbang sistem  adalah 𝐸 2 ∗ = (0.598, 0.399, 0.598) yang  berarti reaksi  mencapai  kesetimbangan  ketika  konsentrasi 𝐻𝐵𝑟𝑂 2 dan 𝐶𝑒(𝐼𝑉) sebesar 0.598 𝑀 ,  dan  konsentrasi 𝐵𝑟 − sebesar 0.399 𝑀 .
Gambar 2. Simulasi untuk 𝑓 = 0.2 .
Gambar  2  hampir  sama  seperti  gambar  4.1,  namun penurunan konsentrasi 𝐻𝐵𝑟𝑂 2 dan kenaikan konsentrasi 𝐵𝑟 − dan 𝐶𝑒(𝐼𝑉) pada  saat 𝑓 = 0.2 terjadi  lebih  lambat.  Setelah itu, konsentrasi 𝐵𝑟 − tetap naik sedangkan konsentrasi 𝐶𝑒(𝐼𝑉) mulai turun, menandakan 𝐵𝑟 − dihasilkan sedikit dan 𝐶𝑒(𝐼𝑉) digunakan  sebagai  reaktan.  Konsentrasi 𝐵𝑟 − mulai  turun sebelum  satuan  waktu  ke  0.9.  Seiring  berjalannya  waktu, konsentrasi 𝐵𝑟 − dan 𝐶𝑒(𝐼𝑉) terus  turun,  sedangkan konsentrasi 𝐻𝐵𝑟𝑂 2 mulai  naik  sebelum  satuan  waktu keempat. Saat konsentrasi 𝐻𝐵𝑟𝑂 2 dan 𝐶𝑒(𝐼𝑉) sama, sebelum satuan  waktu  keenam,  konsentrasi  ketiganya  menjadi konstan,  menandakan  reaksi  menjadi  setimbang.  Jika dibandingkan  dengan  simulasi  sebelumnya,  ketika 𝑓 = 0.2 reaksi menjadi setimbang sedikit lebih cepat. 3. 𝑓 = 0.3
Karena 𝑓 < 0.25 ,  maka  reaksi  stabil.  Titik  setimbang sistem  adalah 𝐸 2 ∗ = (0.396, 0.599, 0.396) yang  berarti reaksi  mencapai  kesetimbangan  ketika  konsentrasi
𝐻𝐵𝑟𝑂 2 dan 𝐶𝑒(𝐼𝑉) sebesar 0.396 𝑀 ,  dan  konsentrasi 𝐵𝑟 − sebesar 0.599 𝑀 .
Gambar 2. Simulasi untuk 𝑓 = 0.3 .
Gambar  3  juga  hampir  sama  dengan  dengan  gambar  4.2, namun  penurunan  konsentrasi 𝐻𝐵𝑟𝑂 2 dan  kenaikan konsentrasi 𝐵𝑟 − dan 𝐶𝑒(𝐼𝑉) pada saat 𝑓 = 0.3 terjadi lebih lambat.  Setelah  itu,  konsentrasi 𝐵𝑟 − tetap  naik  sedangkan konsentrasi 𝐶𝑒(𝐼𝑉) mulai turun, menandakan 𝐵𝑟 − dihasilkan sedikit  dan 𝐶𝑒(𝐼𝑉) digunakan  sebagai  reaktan.  Konsentrasi 𝐵𝑟 − mulai  turun  sebelum  satuan  waktu  ke  0.9.  Jika dibandingkan  dengan  simulasi  sebelumnya,  ketika 𝑓 = 0.3 reaksi menjadi setimbang sedikit lebih cepat.
4. 𝑓 = 0.4
Karena 0.25 < 𝑓 < 1+√2 2 ,  maka  reaksi  tidak  stabil.  Titik setimbang  sistem  adalah 𝐸 2 ∗ = (0.192, 0.797, 0.192) yang berarti  reaksi  mencapai  kesetimbangan  ketika  konsentrasi 𝐻𝐵𝑟𝑂 2 dan 𝐶𝑒(𝐼𝑉) sebesar 0.192 𝑀 ,  dan  konsentrasi 𝐵𝑟 − sebesar 0.797 𝑀 . Limit cycle yang muncul memiliki periode 𝑇 = 6822.44 .
Gambar 4. Simulasi untuk 𝑓 = 0.4 .
Ketika 𝑓 = 0.4 , konsentrasi 𝐻𝐵𝑟𝑂 2 , 𝐵𝑟 − dan 𝐶𝑒(𝐼𝑉) naik dan turun secara berulang – ulang dan sistem tidak kunjung setimbang.  Naik  turunnya  konsentrasi 𝐻𝐵𝑟𝑂 2 , 𝐵𝑟 − dan 𝐶𝑒(𝐼𝑉) ini  dapat  diamati  sebagai  munculnya  perubahan warna  yang  membentuk  suatu  pola  pada  reaksi  saat eksperimen dilakukan. Warna yang muncul pada eksperimen berasal  dari  ion 𝐶𝑒(𝐼𝐼𝐼) dan 𝐶𝑒(𝐼𝑉) .  Oleh  karena 𝐶𝑒(𝐼𝐼𝐼) diubah menjadi 𝐶𝑒(𝐼𝑉) dan 𝐶𝑒(𝐼𝑉) diubah menjadi 𝐶𝑒(𝐼𝐼𝐼)
secara  terus  –  menerus,  muncullah  perubahan  warna  yang akhirnya membentuk suatu pola tertentu. Perubahan 𝐶𝑒(𝐼𝐼𝐼) menjadi 𝐶𝑒(𝐼𝑉) dan  sebaliknya  yang  terus  –  menerus menandakan bahwa kelima reaksi terjadi secara berkesinambungan  dan  terus  –  menerus.  Bidang  fase  untuk 𝑓 = 0.4 ditunjukkan oleh Gambar 5, 6 dan 7.
Gambar 3. Bidang fase untuk 𝐻𝐵𝑟𝑂 2 𝑣𝑠. 𝐵𝑟 − .
Gambar  4. Bidang fase untuk 𝐵𝑟 − 𝑣𝑠. 𝐶𝑒(𝐼𝑉) .
Gambar 5. Bidang fase untuk 𝐻𝐵𝑟𝑂 2 𝑣𝑠. 𝐶𝑒(𝐼𝑉) .
V. KESIMPULAN
Kesimpulan yang didapat dari penelitian ini adalah:
1) Model  oregonator  memiliki  dua  titik  kesetimbangan, yaitu:
𝐸 1 ∗ = (0,0,0) dan
𝐸 2 ∗ = (𝑥 ∗ , 𝑦 ∗ , 𝑧 ∗ ) dimana:
𝑥 ∗ = 1−2𝑓−𝑞+√1+𝑞2+4𝑓2+2𝑞+12𝑞𝑓−4𝑓 2 , 𝑦 ∗ = 2𝑓(1−𝑓−𝑞+√1+𝑞2+4𝑓2+2𝑞+12𝑞𝑓−4𝑓)
1+𝑞−𝑓+√1+𝑞2+4𝑓2+2𝑞+12𝑞𝑓−4𝑓 ,
𝑧 ∗ = 1−2𝑓−𝑞+√1+𝑞2+4𝑓2+2𝑞+12𝑞𝑓−4𝑓 2
Dengan 𝑓 < 1−𝑞 2 dan 1 + 𝑞 2 + 4𝑓 2 + 2𝑞 + 12𝑞𝑓 − 4𝑓 > 0 . 2) Titik setimbang 𝐸 1 ∗ = (0,0,0) tidak stabil sedangkan pada titik setimbang 𝐸 2 ∗ = (𝑥 ∗ , 𝑦 ∗ , 𝑧 ∗ ) terjadi bifurkasi.
3) Bifurkasi  yang  terjadi  pada  titik  setimbang 𝐸 2 ∗ = (𝑥 ∗ , 𝑦 ∗ , 𝑧 ∗ ) adalah bifurkasi Hopf dan sistem mengalami osilasi apabila nilai 𝑓 berada pada rentang 1 4 < 𝑓 < 1 2 .
4) Dari  hasil  simulasi,  untuk  sistem  yang  stabil  semakin kecil  nilai 𝑓 ,  semakin  sedikit  pula  waktu  yang dibutuhkan untuk sistem menjadi setimbang.

## DAFTAR PUSTAKA
[1] N.  Shanks,  “Modeling  Biological  Systems:  The  Belousov– Zhabotinsky Reaction,” Found.  Chem. ,  vol.  3,  no.  1,  pp.  33–53, 2001. [2]
C.  Gray,  “An  Analysis  of  the  Belousov-Zhabotinskii  Reaction,”
2002. [3] R. Chang, Chemistry . New York: Mc Graw – Hill Inc, 2010.
[4] X.  Wang,  Y.  Chang,  and  D.  Xu,  “Bifurcation  Analysis  of  a Belousov–Zhabotinsky  Reaction  Model,” Int.  J.  Bifurc.  Chaos , vol. 25, no. 6, p. 1550093, Jun. 2015.
[5] M. Dodd, “The Belousov-Zhabotinsky Oscillator: An Overview,” 2010.
[6] J.  Tyson, Lecture  Notes  in  Biomethematics ,  1st  ed.  Berlin: Springer, 1976.
[7] Subiono,  “Sistem  Linear  dan  Kontrol  Optimal,”  2013.  [Online]. Available: https://docplayer.info/30430283-Sistem-linear-dan- kontrol-optimal.html. [Accessed: 28-Mar-2019].
[8] J.  Murray, Mathematical  Biology ,  3rd  ed.  New  York:  Springer,
2002.
[9] S. Wiggins, “Introduction,” in
Introduction to Applied Nonlinear
Dynamical Systems and Chaos , New York: Springer-Verlag, 2003, pp. 1–4.