JMP : Vol. 8 No. 2, Des. 2016, hal. 19-32                ISSN 2085-1456
MODEL MATEMATIKA GERAK PENDULUM DENGAN MEMPERTIMBANGKAN GAYA GESEK UDARA
Rukmono Budi Utomo
Universitas Muhammadiyah Tangerang Email: rukmono.budi.u@mail.ugm.ac.id
ABSTRACT. This  paper  aims  to  construct  a  mathematical  model  of  the pendulum  motion  by  considering  the  frictional  force  of  air.  The  study  was conducted by first making a pendulum motion illustration with a length of rope l and  pendulum  mass m   vertically  that  given  deviation  and  speed.  Based  on  that created a mathematical model of motion of the pendulum and its determination to find out  the time t it takes  the pendulum passing through the zero point and the current  time  there  is  no  change  in  the  speed  deviation  of  the  pendulum.  In  this study  will  also  be  given  a  simulation  of  mathematical  models  that  have  been made .
Keywords : Pendulum Motion Mathematical Model, Air Friction
ABSTRAK. Tulisan  ini  bertujuan  untuk  mengkonstruksikan  model  matematika dari  gerak  pendulum  dengan  mempertimbangkan  adanya  gaya  gesek  udara. Penelitian  dilakukan  dengan  terlebih  dahulu  membuat  ilustrasi  gerak  pendulum sederhana  dengan  panjang  tali l dan  beban  massa  pendulum m yang  diikatkan secara  vertikal  dan  diberikan  simpangan  serta  kecepatan.    Berdasarkan  hal tersebut  dibuat  model  matematis  gerak  pendulum  beserta  penentuan t untuk mengetahui  waktu  yang  diperlukan  pendulum  melewati  titik  nol  dan  waktu  saat tidak  ada  perubahan  kecepatan  simpangan  pendulum.  Dalam  penelitian  ini  juga akan diberikan simulasi perhitungan dari model matematis yang telah dibuat.
Kata Kunci : Model Matematis Gerak Pendulum, Gesekan Udara

## 1. PENDAHULUAN
Pendulum atau bandul merupakan benda yang terikat pada sebuah tali dan dapat berayun secara bebas dan periodik. Suatu pendulum sederhana tersusun atas sebuah  bandul  bermassa m yang  diikat  dengan  tali  dengan  panjang l dan ditempatkan menggantung vertikal ke bawah.
Saat tidak diberikan perlakukan apapun, dapat dikatakan bahwa pendulum berada  pada  posisi  diam  atau  posisi  nol.  Pada  posisi  ini  dikatakan  bahwa simpangan  pendulum  adalah  nol.  Lebih  lanjut  ketika  pendulum  ditarik,  maka
pendulum  memiliki  nilai  simpangan  tertentu.  Nilai  simpangan  ini  dapat  berupa positif 0   ataupun negatif 0   .
Lebih  lanjut  apabila  pendulum  diberikan  suatu  kecepatan  tertentu  baik dengan  cara  diayunkan  (arah  kiri) 0   ataupun  didorong  terlebih  dahulu  ke belakang  (arah  kanan) 0   ,  maka  yang  menjadi  pertanyaan  selanjutnya  adalah bagaimana  menjelaskan  perilaku  gerakan  pendulum  tersebut  dalam  suatu  model matematika?  Apabila  pada  gerak  pendulum  berlaku  gaya  gesek  udara  seperti keadaan  sebenarnya,  bagaimana  model  matematikanya?  Selanjutnya  bagaimana cara  menentukan  waktu  agar  pendulum  melewati  titik  nol?  Serta  kapan  waktu yang  menunjukkan  tidak  ada  perubahan  kecepatan  pendulum?  Serangkaian perrtanyaan di atas akan coba penulis ulas dalam penelitian ini. Dalam penelitian ini  juga,  peneliti  tidak  hanya  mengkonstruksikan  model  matematis  gerak pendulum,  namun  juga  dilakukan  simulasi  terhadap  model  yang  telah  terbentuk dengan memasukan beberapa nilai parameter yang terkait.
Penelitian  mengenai  gerakan  pendulum  yang  dilakukan  oleh  peneliti bukanlah  yang  pertama  kalinya.  Awal  mula  studi  tentang  pendulum  ini  diawali dengan  teori  hukum    Hooke  tentang  gaya  pemulih  pada  pegas.  Hukum  Hook sendiri  telah  dikaji  dan  dibahas  secara  lengkap  oleh  (Rychlewski,  1984)  pada artikelnya yang berjudul On Hook’s Law .
Penelitian  tentang  gerakan  pendulum  juga  dilakukan  oleh  (Kecik  et.al, 2006)  yakni  tentang  Chaos  pada  gerakan  pendulum  sebagai  sistem  dengan parametrik  resonansi.  Selain  itu,  (Kecik  et.al,  2012)  juga  meneliti  gerakan autoparametrik dari sistem nonlinear pendulum. Penelitian serupa tentang gerakan pendulum  juga  telah  dilakukan  oleh  (Leeet.al,  1994)  dengan  topik  penelitian tentang  analisis  global  dari  sistem  gerakan  harmonik  pendulum  dengan  internal resonansi

## 2. METODE PENELITIAN
Metode yang dilakukan dalam penelitian ini adalah dengan studi pustaka, yakni  dengan  memperlajari  terlebih  dahulu  teori  tentang  gerakan  pendulum melalui  buku  Agarwal  (2008),  Boyce  (1992),  Haberman  (1971),  Perko  (2000),
Model Matematika Gerak Pendulum                                                                             21
Wiggins  (2003)  serta  buku-buku  lain  yang  dapat  dibaca  dalam  daftar  pustaka . Selain buku, beberapa artikel atau jurnal juga digunakan sebagai  referensi dalam penelitian ini, antara lain artikel J. Rychlewski (1984), Kecik, K., dan Warminski, J (2012), Lee, W.K. dan Hsu, C.S. (1994) dan jurnal lain yang dapat dibaca dalam daftar pustaka .
3. HASIL DAN PEMBAHASAN

## 3.1 Model dan Penyelesaian
Diberikan  pendulum  sederhana  dengan  panjang  tali l dan  massa  bandul m yang diikatkan secara vertikal menggantung ke  bawah. Mula-mula pendulum tersebut  berada  pada  posisi  diam  atau  posisi  nol  sampai  diberikan  simpangan  kecepatan simpangan  sehingga berosilasi seperti pada gambar berikut ini.
Gambar 1. Pendulum berosilasi
Ketika  pendulum  diikat  secara  vertikal  menggantung  ke  bawah  seperti pada  gambar  1  di  atas,  maka  secara  alamiah  massa  pendulum m akan  mendapat pengaruh  gaya  gravitasi g .  Lebih  lanjut,  ketika  pendulum  mulai  diayunkan dengan cara ditarik, maka  terdapat  besarnya sudut  pergerakan pendulum  yang dikenal  dengan  simpangan  pendulum.  Nilai  dianggap  bernilai  positif  apabila pendulum ditarik ke arah kanan, dan sebaliknya bernilai negatif apabila pendulum ditarik  ke  arah  kiri.  Sebelum  diayunkan,  didefinsikan  yang  merupakan kecepatan  pendulum.  Nilai  ini  juga  dapat  berupa  nilai  positif  ataupun  negatif tergantung  cara  memberikan  kecepatan  pada  pendulum  tersebut  (didorong  atau
ditarik). Saat pendulum bergerak, maka pergerakan pendulum akan mendapatkan gaya gesek udara yang merupakan gaya redam dan disimbolkan dengan d a dt  .
Sebagaimana  diketahui  bahwa  gaya  gesek  ataupun  gaya  redam  selalu berlawanan dengan arah gerakan pendulum itu sendiri sebagaimana hukum Hooke menjelaskan tantang gaya gesek lantai yang berlawanan terhadap tarikan maupun dorongan pada pegas secara horizontal.
Dalam penelitian ini diasumsikan bahwa mula-mula pegas akan ditarik ke arah  kanan  sehingga  nilai  simpangan  pendulum  adalah  positif.  Lebih  lanjut, kecepatan  yang  diberikan  pada  pendulum  juga  dalam  hal  ini  dimisalkan  positif yakni 0 d dt   yang  memberi  arti  bahwa  pendulum  diperlakukan  tarikan  atau gerakan  ke  arah  kanan  dan  bukan  berupa  dorongan.  Berdasarkan  hal  tersebut, model  matematis  gerak  pendulum  sebagaimana  diilustrasikan  pada  gambar  1  di atas dapat dituliskan sebagai berikut :
sin 0 ml a mg         1
Karena d dt    yang  merupakan  kecepatan  dari  gerak  pendulum  dan
2 d dt    yang  merupakan  percepatan  gerak  pendulum,  maka  berdasarkan  hal tersebut persamaan   1 dapat dituliskan kembali sebagai berikut :
2 sin 0 d d ml a mg dt dt  
      2
Apabila  didefinisikan a ml   dan g l   ,  maka  persamaan   2 di  atas dapat dituliskan sebagai berikut
2 sin 0
d d dt dt           3
Misalkan 1 x   dan 2 x   ,  maka  model  persamaan   3   dapat  dituliskan kembali sebagai sistem model simutan sebagai berikut
1 2 2 1 2 sin x x x x x         4 Titik kesetimbangan
  1 2 , x x dari sistem   4 diperoleh dengan menyamadengankan  nol 1 x dan 2 x ,  dan  berdasarkan  hal  tersebut  diperoleh  titik kesetimbangan sistem   4 yakni   0, 0  dan   , 0 n  dengan n N  .
Karena  sistem   4 merupakan  model  sistem  nonlinear,  maka  akan dilakukan  linearisasi  sistem  disekitar  titik  kesetimbangan (equilibrium) . Linearisasi  dilakukan  dengan  Jakobian,  dan  berdasarkan  hal  tersebut  diperoleh sistem sebagai berikut
          1 1 2 1 1 2 1 2 1 1 1 2 1 2 2 2 1 2 2 1 2 1 2 , , 0 1 , cos
, , f x x f x x x x x x Df x x x x x f x x f x x x x      
                                    
 

  5
dengan   1 1 2 2 , f x x x  dan   2 1 2 1 2 , sin f x x x x      .
Selanjutnya  akan  dilakukan  analisis  kestabilan  terhadap  sistem  model terlinearisasi  disekitar titik kesetimbangan   0, 0 dan untuk analisis disekitar titik kesetimbangan   , 0 n  ,  dapat  dilakukan  dengan  cara  dan  langkah  yang  sama.
Berdasarkan  sistem  linearisasi  pendulum  sesuai  persamaan   5 ,  maka  analisis kestabilan  sistem  model  terlinerisasi  disekitar  titik  kesetimbangan   0, 0   dapat diuraikan sebagai berikut:
  1 2 0 1 0, 0 x Df x                    6
Selanjutnya nilai eigen dari sistem   6  diperoleh nilai-nilai eigen sebagai berikut :
2 1,2 4 2 2           7
Solusi   1 x t dan   2 x t dari sistem   6 berdasarkan nilai eigen   7 masing-masing dituliskan sebagai berikut:
  2 2 4 4 2 2 2 2 1 1 2 t t x t C e C e                                 8
  2 2 4 4 2 2 2 2 2 2 2 1 2 4
4
2 2 2 2 t t x t C e
C e              
                                              9
Berdasarkan  nilai-nilai  eigen  pada  persamaan   7 di  atas,  dapat  dikaji  tiga kondisi berdasarkan
2 4    yang selanjutnya dijelaskan sebagai berikut:
1) Kondisi 2 4 0    
Apabila  nilai 2 4 0     ,  maka  diperoleh    nilai  eigen  sistem   6
1 2 2       .  Dengan  mengingat  bahwa 0 a ml    ,  maka 0 2    .
Berdasarkan  hal  tersebut  solusi  matematis   1 x t dan   2 x t berdasarkan kondisi 1 ini dapat dituliskan sebagai berikut
  2 2 1 1 2 t t x t C e C te         10
  2 2 2 2 1 2 2 2 2 t t t t t x t
C e C e tC e         
      11 Dengan  syarat  awal   1 10 0 x x  dan   2 20 0 x x  ,  maka  solusi   10 dan
  11 masing-masing dapat dituliskan kembali sebagai berikut
  2 2 2 1 10 20 2 t t t x t x e te x te                
  12
  2 2 2 2 2 10 20 4 2 t t t x t te x e te x       
             13
Lebih lanjut berdasarkan solusi   12 dan   13  akan ditentukan waktu agar pendulum berada pada titik nol dan waktu ketika tidak terdapat perubahan kecepetan  pendulum.  Masing-masing  waktu  tersebut  dicari  dengan
menyamadengankan  nol   1 x t dan   2 x t sehingga  berdasarkan  hal tersebut  diperoleh  waktu
10 1
10 20 2 x t x x     dan
20 2 2 10 20 4 2 x t x x     .
Masing-masing  nilai 1 t dan 2 t akan  ada  apabila
10 20 2 x x    dan 2 10 20 4 2 x x     .
2) Kondisi 2 4 0    
Apabila  nilai 2 4 0     ,  maka  diperoleh  nilai  eigen  sistem   6 2 1 4 2 2         dan 2 2 4 2 2         dengan 2 1 0     .
Berdasarkan  hal  tersebut  solusi  matematis   1 x t dan   2 x t berdasarkan kondisi 2 ditunjukkan dengan persamaan   8 dan   9 . Dengan syarat awal   1 10 0 x x  dan   2 20 0 x x  ,  maka  solusi   8 dan   9   masing-masing dapat dituliskan kembali sebagai berikut
      1 2 1 2 10 20 20 1 10 2 1 1 t t x t x x e x x e                  14
      1 2 2 1 2 1 2 1 2 10 2 1 20 2 1 1 t t t t x t e e x e e x                        15 .
Lebih  lanjut  berdasarkan  solusi   14 dan   15   waktu  agar  pendulum berada pada titik nol dan waktu ketika tidak terdapat perubahan kecepetan
pendulum adalah 2 10 20 1 2 1 1 10 20 1 ln x x t x x               dan
1 2 10 1 20 2 2 1 1 2 10 1 20 1 ln x x
t x x   
    
          .  Kedua  waktu 1 t dan 2 t akan  ada  apabila
2 10 20 1 10 20 1 x x x x      dan 1 2 10 1 20 1 2 10 1 20 1 x x x x          . 3) Kondisi 2 4 0    
Apabila  nilai 2 4 0     ,  maka  nilai  eigen  sistem   6 adalah  imajiner
yakni 1 i      dan 2 i      dengan
2     dan 2 4 2      .
Berdasarkan hal tersebut solusi matematis  
1 x t dan   2 x t dengan syarat
awal   1 10 0 x x  dan   2 20 0 x x  berdasarkan kondisi 3 dituliskan kembali sebagai berikut
  1 10 20 sin sin e e x t x t x t          
  16       2 2 2 10 20 sin sin cos e e x t x t x t t                   17
Lebih  lanjut  berdasarkan  solusi   16 dan   17   waktu  agar  pendulum berada pada titik nol dan waktu ketika tidak terdapat perubahan kecepetan
pendulum. adalah 10 1 20 10 1 arctan x t x x 
           dan
  20 2 2 2 10 1 arctan x t x         
        . Kedua  waktu 1 t dan 2 t akan  ada  apabila  periode 10 20 10 x x x     dan
  20 2 2 10 x x        bukan 2  atau kelipatannya.

## 3.2 Simulasi
Simulasi  dilakukan  untuk  tiap  kondisi  yang  dibahas  pada  penelitian  ini,  dan masing-masing dijelaskan pada uraian di bawah ini
1) Kondisi 2 4 0
   
Ada  banyak  sekali  nilai  parameter  dan  agar  nilai 2 4 0     ,
namun  pada  intinya  kondisi  ini  terjadi  apabila  nilai 2 4    .  Sebagai
pengingat  bahwa  nilai a ml   dan g l   ,  dan  misalkan  dalam  hal  ini  di
ambil parameter 20, 1, 10 a m l    dan 2 10 / g m s  atau 2   dan 1   dan dengan  kondisi  awal     1 2 0 0 10 x x   .  Berdasarkan  hal  tersebut diperoleh simulasi sebagai berikut.
Gambar 2.   Plot   1 x t dan   2 x t
Gambar 3. Potret fase Sistem   6
Sistem   6 kondisi 1 kondisi 1
Gambar  2  diatas  menunjukkan  grafik  plot  sistem  pendulum terlinearisasi di sekitar titik kesetimbangan   0, 0 . Grafik berwarna merah menunjukkan    simpangan  atau  posisi  pendulum   1 x t   dan  garfik berwarna hijau menunjukkan kecepatan pendulum   2 x t   .
Dari gambar 2 terlihat bahwa pendulum yang digerakan suatu saat akan berhenti dengan kata lain simpangan  dan kecepatan  akan habis atau  nol.  Hal  demikian  dapat  terjadi  karena  adanya  gaya  gesekan  udara d a dt  dengan  nilai  konstanta  gaya  sesek  udara  dalam  simulasi  ini
ditetapkan 20 a  .
Lebih lanjut gambar 3 menunjukkan potret fase sistem   6 dengan nilai parameter yang diberikan pada kondisi 1 ini. Potret fase pada gambar 3 menunjukkan bahwa titik kesetimbangan   0, 0 stabil asimptotik dengan
28                                                                                                      Rukmono Budi Utomo
tipe sink . Kestabilan asimptotik ini diperkuat dengan nilai eigen kondisi 1 yakni 1 2 0 2        .
2) Kondisi 2 4 0    
Agar  terbentuk  kondisi 2 4 0     ,  misalkan  dalam  hal  ini  di  ambil parameter 20, 0.5, 10 a m l    dan 2 10 / g m s  atau 4   dan 1   dan dengan  kondisi  awal     1 2 0 0 10 x x   .  Berdasarkan  hal  tersebut diperoleh simulasi sebagai berikut
Gambar 4.  Plot   1 x t dan   2 x t Gambar 5. Potret fase Sistem   6
Sistem   6 kondisi 2 kondisi 2
Gambar  4  diatas  menunjukkan  grafik  plot  sistem  pendulum terlinearisasi di sekitar titik kesetimbangan   0, 0 . Grafik berwarna merah menunjukkan    simpangan  atau  posisi  pendulum   1 x t   dan  garfik berwarna hijau menunjukkan kecepatan pendulum   2 x t   . Dari gambar 4 terlihat bahwa pendulum yang digerakan suatu saat akan berhenti dengan kata  lain  simpangan  dan  kecepatan  akan  habis  atau  nol.  Hal demikian  dapat  terjadi  karena  adanya  gaya  gesekan  udara d a dt  dengan nilai  konstanta  gaya  sesek  udara  dalam  simulasi  ini  ditetapkan 20 a  . Meski  pendulum  akan  berhenti  sebagimana  pada  kondisi  1,  namun
perilaku  pendulum  menuju  berhenti  sangatlah  berbeda.  Perbedaan  ini dapat dilihat pada potret fase gambar 5.
Potret  fase  sistem   6 sebegaimana  ditunjukkan  gambar  6 menunjukkan  titik  kesetimbangan   0, 0 stabil  asimptotik  dengan  tipe node .    Tipe  ini  berbeda  dengan  potret  fase  kondisi  1  yang  bertipe  node. Kestabilan  asimptotik  ini  diperkuat  dengan  nilai  eigen  kondisi  2  yakni
2
1 4 2 2         dan 2 2 4 2 2         dengan 2 1 0    
3) Kondisi 2 4 0    
Agar terbentuk kondisi 2 4 0     , misalkan dalam hal ini di ambil parameter 20, 2, 10 a m l    dan 2 10 / g m s  atau 1   dan 1   dan
dengan  kondisi  awal     1 2 0 0 10 x x   .  Berdasarkan  hal  tersebut diperoleh simulasi sebagai berikut:
Gambar 6.  Plot   1 x t dan   2 x t Gambar 7. Potret fase Sistem   6
Sistem   6 kondisi 3                               Kondisi 3
Gambar  6  diatas  menunjukkan  grafik  plot  sistem  pendulum terlinearisasi di sekitar titik kesetimbangan   0, 0 . Grafik berwarna merah menunjukkan    simpangan  pendulum   1 x t   dan  garfik  berwarna  hijau menunjukkan  kecepatan  pendulum   2 x t   .  Dari  gambar  6  terlihat
bahwa pendulum yang digerakan suatu saat akan berhenti atau dengan kata lain  simpangan  dan  kecepatan  akan  habis  atau  nol.  Hal  demikian dapat  terjadi  karena  adanya  gaya  gesekan  udara d a dt  dengan  nilai konstanta  gaya  sesek  udara  dalam  simulasi  ini  ditetapkan 20 a  .  Cara pendulum  ini  berhenti  berbeda  dengan  cara  pendulum  berhenti  pada kondisi  1  dan  2.  Pada  kondisi  3  ini  sebelum  pendulum  berhenti,  terlebih dahulu  pendulum  akan  berosilasi.  Hal  ini  dikarenakan  nilai  eigen  yang bernilai  imajner.  Berdasarkan  hal  ini  pula,  potret  fase  yang  ditunjukkan pada gambar 7 menjelaskan bahwa titik   0, 0 stabil asimptotik dengan tipe potert fase adalah spiral

## 4. KESIMPULAN DAN SARAN
Terdapat beberapa kesimpulan dan saran dalam penelitian ini ynag masing-masing diuraikan sebagai berikut.

## 4.1 Kesimpulan
Kesimpulan yang diperoleh adalah :
1) Model  matematika  yang  menjelaskan  gerakan  pendulum  dengan  adanya gaya gesek udara ditunjukkan pada persamaan   1 atau   4 dengan sistem yang  telah  terlinierisasi  disekitar  titik  kesetimbangan   0, 0 ditunjukkan pada persamaan   6  dengan solusi   1 x t dan   2 x t sistem   6  ditunjukkan pada persamaan   8 dan   9 .
2) Nilai eigen dari  sistem   6  memberikan nilai 2 1,2 4 2 2        
3) Berdasarkan nilai eigen pada poin 2 diatas, dibagi tiga kondisi berdasarkan nilai 2 4    ,  yakni  kondisi  saat
2 4 0     ,  kondisi  saat
2 4 0     dan 2 4 0    
4) Untuk  kondisi 2 4 0     ,  solusi   1 x t dan   2 x t sistem   6 ditunjukkan  pada  persamaan   12 dan   13 dengan  waktu  pendulum mencapai  titik  nol 1 t dan  waktu  saat  tidak  terdapat  perubahan  kecepatan pendulum 2 t yang dapat dilihatpada pembahasan kondisi 1
5) Untuk  kondisi 2 4 0     ,  solusi   1 x t dan   2 x t sistem   6 ditunjukkan  pada  persamaan   14 dan   15   dengan  waktu  pendulum mencapai  titik  nol 1 t dan  waktu  saat  tidak  terdapat  perubahan  kecepatan pendulum 2 t yang dapat dilihat pada pembahasan kondisi 2
6) Untuk  kondisi
2 4 0     ,  solusi   1 x t dan   2 x t sistem   6 ditunjukkan  pada  persamaan   16 dan   17 dengan  waktu  pendulum mencapai  titik  nol 1 t dan  waktu  saat  tidak  terdapat  perubahan  kecepatan pendulum 2 t yang dapat dilihat pada pembahasan kondisi 3
7) Berdasarkan simulasi yang dilakukan untuk ketiga kondisi di atas, terlihat bahwa  pendulum  yang  digerakan  dengan  diberikan  simpangan  dan kecepatan suatu saat pendulum tersebut akan berhenti. Perilaku pendulum menuju  berhentinya  untuk  masing-masing  kondisi  berbeda-beda tergantung nilai eigenya dan dapat dilihat pada pembahasan simulasi.

## 4.2 Saran
Saran yang dapat diberikan adalah :
1) Perlu  dilakukan  simulasi  terhadap  kondisi  3  saat  nilai  eigen  sistem   6 adalah imajiner murni.
2) Perlu  diselidiki  linearisasi  sistem   4   disekitar  titik  kesetimbangan
 
, 0 n  dengan 1, 2, . n 
3) Perlu  dikembangkan  sistem  model  matematis  gerak  pendulum  dengan mempertimbangan faktor lain yang terjadi pada kehidupan  sebenarnya
4) Perlu dikembangkan model matematis gerak pendulum 3 Dimensi.
32                                                                                                      Rukmono Budi Utomo

## UCAPAN TERIMAKASIH
Ucapan  terimakasih  ditujukan  kepada  kampus  Universitas  Muhammadiyah Tangerang dan kampus Institut Teknologi Bandung atas support sehingga tulisan ini  dapat  selesai  dan  penulis  dapat  mengikuti  seminar  nasional  matematika  dan terapannya di kampus Universitas Jendral Soedirman Tanggal 3 Desember 2016

## DAFTAR PUSTAKA
Agarwal, R. D, dan Regan, D. O., An Introduction to Differential Equation ,
Springer, New York, 2008.
Boyce, W.E. dan  Richard, Elementary Differential Equation and Boundary Value Problem , fifth edition, John Wiley and Sons, Inc., 1992.
Haberman,  R .,  Mathematical  Models  :  Mechanical  Vibration,  Population Dynamics and Traffic Flow , Prentice Hall, Inc., 1971.
Kecik,  K.,  dan  Warminski,  J. Chaos  in  Mechanical  Pendulum-Like  System  near Main Parametric Resonance, Procedia IUTAM, 5 (2012), 249-258. Lee, W.K. dan Hsu, C.S., A Global Analysis of an Harmonically Excited Spring- Pendulum  System  with  Internal  Resonance , Journal  of  Sound  and Vibration, 171 (3) (1994) , 335-359.
Perko,  Lawrence., Differensial  Equation  and  Dynamical  System, Third  Edition, Springer-Verlag, New York, 2000.
Rychlewski,  J., On  Hook’s  Law ,  Pergamon  Press,  Ltd.,  PMM  U.S.S.R., 48 (3) 1984, 303-314.
Warminski, J. dan Kecik, K. Autoparametric Vibration of Nonlinear System with Pendulum, Mathematical  Problems  in  Engineering  Volume  2006,  2006, Article  ID  80705,  19  pages, http://dx.doi.org/10.1155/MPE/2006/80705 , diakses pada 26 November 2016.
Wiggins,  S., Introduction  to  Applied  Nonlinear  Dynamical  Systems  and  Chaos , Springer-Verlag, New York, 2003
Widowati  dan  Sutimin, Pemodelan  Matematika:  Analisis  dan  aplikasinya , Penerbit Undip Press, Semarang, 2013 .