sekali magnet yang di sini. Karena

Eko Sulisstya/ Simulasi Model M Ising 2 Dim mensi dengan Algoritma A Metrop polis pada Lemb barkerja Excel

35

Simu ulasi Moodel Ising 2 Diimensi dengan d Algorittma Meetropoliis padaa Lembaarkerjaa Excel Eko Sullistya Jurusan Fiisika, FMIPA UGM, U Sekip Uttara Kotak Pos Bls B 21 Yogyakaarta 55281 [email protected]

Abstrak – Pengujiann dan peneraapan algoritm ma Metropoliss telah dilakkukan pada llembarkerja Excel. E Pengujjian dilakukann dengan carra mensimulaasikan pembalikan spin da an menghitungg magnetisassi bahan ferom magnetik berrupa model Ising 2D. Ukurran kisi bisa dipilih d mulai dari 10×10 sampai dengaan 100×100. Magnetisasi dihitung sebaagai fungsi suuhu, mulai daari 0,1 sampaai dengan 5 satuan s J/kB untuk ukuran kisi 20×20. Hasil yang diperoleh d adaalah semakin tinggi suhu, maka magnnetisasi semakkin kecil, meenurun tajam m dan mendeekati 0 pada T = 2,26 J/k J B. Keuntunggan menggunnakan Excel adalah a bahwaa proses pem mbalikan spin dapat divisuualisasikan tanpa memerluukan modul grrafik, dan setiaap langkah koomputasi dapaat ditampilkan n. Kata kun nci: simulasi, feromagnetikk, model Ising,, algoritma Metropolis M Abstractt – Testing annd application of the Metroppolis algorithm has been done do on the Exxcel worksheet. Testing is done d by simulaating the spin reversal andd calculating the t magnetizattion of ferrom magnetic mateerials in the fo orm of a 2D Issing model. The T lattice sizze can be seelected from 10×10 to 10 00 ×100. Maggnetization iss calculated as a a functionn of temperatture, ranging from f 0.1 to 5 units of J/kB for lattice siize 20×20. Thhe result is thhe higher the temperature, the magnetizzation is gettinng smaller, decreasing d shaarply and app proached zeroo at T = 2.266 J/kB. The ad dvantage of ussing Excel is that the spinn reversal proocess can be visualized without w the neeed for graphics module, and a each stepp of computattion can be dissplayed. Key worrds: simulationn, ferromagneetics, Ising moodel, Metropolis algorithm I. PEND DAHULUAN Bahann feromagnetiik tersusun oleh banyakk sekali atom, yaang menyumbbangkan mom men magnet berasal dari gerakkan elektron bebas. b Magneetisasi bahann paramagnettik ditentukaan oleh keadaan//arah dari momen m dipol magnet atoom-atom penyusunnnya. Jika baanyak momenn dipol magnnet yang searah, maka m magneetisasinya besar, sebaliknnya jika momen dipolnya d sebaanding antara yang searah dengan yang berllawanan, makka magnetisasiinya kecil. Secaraa statistik, maagnetisasi bahhan adalah nillai ratarata darii keadaan miikroskopis peenyusunnya. Karena jumlahnyya sangat besar, maka nilaai rata-rata tiddak bisa (tidak mungkin, m atauu sukar) dipperoleh dengaan cara menghituung masing-m masing nilai, sehingga caranya adalah melakukan m sam mpling secaraa acak beberaapa nilai untuk mewakili m keseeluruhan, kem mudian diambbil nilai rata-ratannya, yaitu denngan metode Monte M Carlo. Nilai harap h yang diicari dengan metode Montte Carlo adalah juumlah dari selluruh hasil kalli nilai peluanng untuk memperooleh suatu nilaai dengan nilaai fungsi dari bilangan b acak yang dibangkitkaan dibagi denggan ba-nyaknyya cacah samplingg. Cara ini bissa dilakukan jika fungsi distribusi peluang dan d bentuk fuungsi yang diccari berupa fuungsi sederhana.

SAN TEORI II. LANDAS A. Feromagnetik Dalam feromagnetism f me, medan magnit yang y ditimbulkan oleh bahan disebabkan oleh spin dari d elektron yanng tidak berppasangan. Tiaap spin terseebut “senang” meenunjuk pada arah yang saama dengan arah a spin tetanggaanya. Alasan mengapa beg gitu, sepenuhhnya adalah fenom mena mekanikaa kuantum, yaang tidak dibaahas di sini. ya terjadi pada p Penyearahaan spin-spin itu biasany daerah-daerahh kecil (sekiitar 10−3 mm m3) yang diseebut domain. Tiapp domain bissa terdiri darri milyaran sppin, semuanya searah, tetapi tiiap-tiap domaain arahnya acak a sehingga seccara makroskopik medan n magnet yang y ditimbulkan bahan itu 00. Itulah seb babnya menggapa magnit perman nen. Gambarr 1 sepotong beesi bukan m menampilkann contoh domaain feromagneetik.

Gamb bar 1. Domain Feromagnetik [Griffiths]. [

Prosiding Peertemuan Ilmiah h XXVI HFI Jatteng & DIY, Purrworejo 14 Aprill 2012 ISSN : 0853 3-0823

36

Eko Sulistya/ Simulasi Model Ising 2 Dimensi dengan Algoritma Metropolis pada Lembarkerja Excel

Salah satu hal yang dapat menghilangkan arah spin yang seragam adalah gerak acak akibat perubahan suhu. Ini terjadi tepat pada suhu Curie, saat sifat feromagnetik tiba-tiba berubah menjadi paramagnetik. B. Magnetisasi Menurut distribusi Boltzmann, peluang bahwa bahan berada pada suatu keadaan S dinyatakan dengan [2]

P( S ) ∝ e − E ( S )

k BT

(1)

dengan E(S) adalah tenaga pada keadaan S, kB adalah konstanta Boltzmann, dan T adalah temperatur mutlak pada keadaan itu. Tidak mudah untuk melakukan sampling dengan bentuk distribusi peluang seperti di atas, sehingga dilakukan sampling dengan cara yang lain. Ada beberapa algoritma yang bisa digunakan untuk menguji model Ising, antara lain: heat bath algorithm, Metropolis algorithm, Single-bond algorithm, Swendsen-Wang algorithm, Wolff’s algorithm, dan Propp-Wilson algorithm. C. Algoritma Metropolis Dalam paper aslinya dikatakan bahwa, metode Metropolis[3] merupakan varian dari metode Monte Carlo, yaitu tidak memilih distribusi secara acak kemudian dibobot dengan exp(−E/kBT), namun memilih distribusi dengan pe-luang exp(−E/kBT) dan membobotnya secara merata. Caranya adalah, sejumlah N partikel diletakkan pada suatu distribusi sembarang, misalnya pada tiap kisi yang beraturan. Kemudian tiap partikel secara berurutan dipindahkan. Pada setiap pemindahan satu partikel, perubahan energi dE dari distribusi itu dihitung. Jika dE bernilai negatif, berarti terjadi penurunan energi akibat pemindahan satu partikel itu, maka pemindahan itu dibolehkan dan partikel diletakkan pada posisi yang baru. Jika dE bernilai positif, pemindahan diperbolehkan dengan peluang exp(−E/kBT). D. Model Ising 2 dimensi Model Ising diambil dari nama Ising (1925), yang menyelesaikan model yang diusulkan oleh Lenz (1920) untuk mempelajari fase transisi feromagnetik pada temperatur Curie. Ising menyelesaian kasus untuk 1 dimensi, kemudian Onsager menyelesaikan untuk kasus 2 dimensi 20 tahun kemudian. Gambar 2 merupakan bagian dari lattice ukuran N×N, dengan keadaan tiap sel dinyatakan dengan spin-up atau spin-down. Setiap perubahan spin pada suatu sel, akan menyebabkan perubahan energi sistem

. Gambar 2. Model sederhana bahan feromagnetik 2 dimensi.

Perbandingan distribusi peluang antara dua keadaan A dan B, diperoleh berdasarkan persamaan (1) : P ( B ) exp(− E ( B ) / k BT ) ⎞ (2) = = exp ⎛⎜ − dE k BT ⎟⎠ ⎝ P ( A) exp(− E ( A) / k BT ) dengan dE = E(B) − E(A). A adalah keadaan sebelum diubah dan B keadaan setelah diubah. Jika dE negatif, maka dari persamaan (2), P(B) lebih besar daripada P(A), sehingga keadaan B bisa diterima. Sedangkan jika dE positif, maka terjadi sebaliknya, P(B) lebih kecil daripada P(A), namun tidak berarti perubahan langsung ditolak. Ditolak atau diterima bergantung pada nilai exp(−E/kT), dibandingkan dengan suatu nilai acak antara 0 dan 1. Dengan kata lain, tidak apa-apa melakukan kesalahan langkah, dengan peluang sebesar exp(−E/kT). Dalam model Ising 2 dimensi seperti pada Gambar 2, tanpa medan magnet luar, dE ditentukan oleh interaksi antara sel yang diubah spin-nya dengan nilai-nilai spin sel tetangganya (atas, bawah, kiri, dan kanan) menurut persamaan

dE = −J ∑ Si S j

(3)

i≠ j

dengan J suatu konstanta, disebut exchange constant, nilainya positif untuk interaksi feromagnetik, negatif untuk interaksi antiferomagnetik, dan 0 jika tidak ada interaksi. Si adalah spin yang diubah, sedangkan Sj adalah spin di sel atas, bawah, kiri, dan kanan. Magnetisasi dihitung berdasarkan persamaan [4]

M = ∑ Si

(4)

i

III. METODE PENELITIAN/IMPLEMENTASI Dibuat satu file Excel yang terdiri dari 10 lembarkerja yang berturut-turut diberi nama 10×10, 20×20, dan seterusnya sampai 100×100. Masing-masing lembarkerja mewakili ukuran kisi sesuai dengan namanya. Pilihan berapa ukuran kisi yang akan divisualisasikan dan untuk mengisikan variabel masukan dilakukan pada userform seperti ditunjukkan pada Gambar 3.

Prosiding Pertemuan Ilmiah XXVI HFI Jateng & DIY, Purworejo 14 April 2012 ISSN : 0853-0823

Eko Sulisstya/ Simulasi Model M Ising 2 Dim mensi dengan Algoritma A Metrop polis pada Lemb barkerja Excel

37

sistem termaagnetisasi. Paada keadaan ini, jika iteerasi diteruskan, salah satu spinn bisa saja membalik, m nam mun hanya sesaat, kemudian keembali sejajar dengan spin-sspin tetangganya. Hasil-hasill simulasi yanng ditampilkan n pada Gambaar 4 dan Gambarr 5 sesuai deengan simulaasi-simulasi yang y dihasilkan oleh Jim Ma [55] yang meng ggunakan Pythhon, dan Velasco [6] yang mengggunakan Java. Gambarr 3. Form untukk pemilihan ukuuran kisi, dan masukan. m

Keadaaan awal spinn setiap sel dittentukan secaara acak. Keadaan acak diperooleh dengan cara c membanngkitkan satu bilaangan acak. Jiika bilangan acak yang diiperoleh bernilai lebih besar dari 0,5 maaka spin di sel itu ditentukaan bernilai +11 dan sel dibeeri warna puttih. Jika bilangan acak bernilaii kurang dari 0,5 maka spiin di sel itu ditenttukan bernilai –1 dan sel dibberi warna meerah. Langkah--langkah/algooritmanya adallah sebagai beerikut: 1. Pilih P suatu sel secara acak. 2. Hitung H perubaahan energi dE dengan perrsamaan (33) 3. Jika dE < 0 , terima t perubaahan (spin dibaalik), ke laangkah 5. 4. Jika dE > 0 , bangkitkaan bilangann acak jika terima perubahan 0 <η <1, e − dE / T ) > η , exp( 5. Diulang D sebanyyak pengulanggan yang diinnginkan 6. Hitung H magnettisasi dengan persamaan p (4))

a. Keaddaan awal

b. Setelah h 500 kali iteraasi

c. Setelah 10000 kali iterasi

d. Setelah h 1500 kali iteraasi

e. Setelah 20000 kali iterasi

f. Setelah h 2500 kali iteraasi

g. Setelah 30000 kali iterasi

h. Setelah h 3500 kali iteraasi

i. Setelah 4000 kali iterasi

j. Setellah 4500 iterasi

Untuk keperluan koomputasi, diggunakan J = 2 , dan satuan unntuk T adalah J/kB. SIL DAN PEMBAHASAN N IV. HAS A. Visuallisasi Spesifi fikasi kompuuter yang digunakan d addalah : Prosesor AMD Athlonn(tm) II Neo K235 K Dual-Coore 1,30 GHz, dann RAM 5 GB.. Gambaar 3 menampiilkan urutan keadaan k untukk ukuran kisi 20× ×20 sel. Padda keadaan awal, spin tiap t sel dibangkittkan secara acak. Suhu dipilih pada T = 1 dengan satuan joule/kB. Dari Gambar G 4 terrlihat bahwa terjadi penyyearahan spin dippol magnit yang berdekkatan, menunnjukkan pembentuukan domain magnit. Setelaah 10.000 kali iterasi, diperolehh keadaan yanng ditunjukkann pada Gambaar 5a. Gambaar 5b dan seterusnya s menampilkan keadaan k akhir settelah 10.000 kali k iterasi unntuk suhu-suhhu yang berbeda. Setiap pemillihan suhu yaang berbeda, dimulai dari keaddaan awal yang acak. Gambaar 5 menunjuukkan bahwa semakin tingggi suhu, ukuran domain semakiin kecil, dan pada p suhu yanng tinggi keadaan spin menjadi semakin accak. Secara numerik, n dengan semakin s tinggginya suhu, algoritma a Meetropolis menerim ma pembalikann spin yang anntiparalel. Pada suhu rendah terjadi t pengellompokan spiin-spin yang searah, dan bisa terjadi penyeearahan keseluuruhan spin sehingga s

Gam mbar 4. Urutan kkeadaan spin kiisi 20×20.

Prosiding Peertemuan Ilmiah h XXVI HFI Jatteng & DIY, Purrworejo 14 Aprill 2012 ISSN : 0853 3-0823

38


T = 2, 26 ×

5, 0 ×10−21 J = 409 K 2 × kB

kB = 1,38 ×10−23 J K , dengan Boltzmann.

a. T = 1

b. T = 1,5

c. T = 2

d. T = 2,5

adalah

konstanta

Gambar 6. Grafik magnetisasi ternormalisasi sebagai fungsi suhu.

e. T = 3

g. T = 5

Nilai suhu yang diperoleh pada Gambar 6 dekat dengan suhu Curie untuk Nikel, yaitu 358K. Di atas suhu ini, bahan feromagnetik berubah menjadi paramagnetik dengan spin dipol magnet yang arahnya acak, ditunjukkan dengan mengecilnya magnetisasi.

f. T = 4

h. T = 10

Gambar 5. Keadaan akhir untuk beberapa suhu dan 10.000 iterasi.

B. Magnetisasi Magnetisasi dihitung dengan persamaan (4) dan dinormalisasi. Jika seluruh spin searah maka magnetisasinya 1, atau 100%. Grafik nilai mutlak magnetisasi ternormalisasi sebagai fungsi suhu, mulai dari T = 0,1 J/kB sampai dengan T = 5 J/kB dengan interval 0,1 ditunjukkan pada Gambar 6, cocok dengan yang dihasilkan oleh Ishita Argawal [7], bahwa magnetisasi turun tajam pada T = 2,26 J/kB. Dari persamaan (3), orde besar dE sama dengan orde besarnya konstanta interaksi J. Nilai J untuk bahan Ni telah diperoleh oleh P. M. Hemenger dan H. Weik [8], yaitu J = 5,0×10-21 J. Dalam perhitungan untuk komputasi ini telah digunakan nilai J = 2, sehingga bisa diperkirakan besar suhu saat magnetisasi mulai turun, yaitu

C. Waktu komputasi Tabel 1 menampilkan waktu yang diperlukan untuk menjalankan algoritma Metropolis, mulai dari ukuran kisi 10x10 sampai 100x100, suhu T = 1 J/kB, dan 5.000 iterasi. Microsoft Excel dirancang untuk multitasking, sehingga jika dijalankan bersama dengan aplikasi lain, maka akan berpengaruh pada waktu eksekusinya. Hasil yang disajikan pada Tabel 1 bisa saja berbeda-beda untuk tiap komputer, selain karena prosesornya, juga karena aplikasi lain yang dijalankan bersama Excel. Karena program Ising yang dibuat dengan bahasa pemrograman lain tidak membatasi simulasi pada jumlah iterasi tertentu, maka waktu komputasi tidak bisa dibandingkan dengan waktu komputasi dari simulasi Ising yang dibuat dengan Excel ini. Tabel 1. Data waktu komputasi untuk berbagai ukuran kisi. Ukuran kisi Waktu (detik) 10×10 1,02 20×20 1,14 30×20 1,46 40×40 1,76 50×50 2,27 60×60 2,38 70×70 2,44 80×80 2,47 90×90 2,55 100×100 2,98



Tabel 2 menampilkan waktu komputasi untuk ukuran kisi 50×50 untuk beberapa jumlah iterasi. Tabel 2. Waktu komputasi untuk beberapa jumlah iterasi. Jumlah iterasi Waktu (detik) 1000 0,59 2000 0,89 4000 1,66 6000 2,26 10.000 3,06 20.000 5,12 50.000 9,81 75.000 13,37 100.000 16,39 200.000 29,86

[4]

[5] [6]

Nampak bahwa untuk jumlah iterasi yang cukup besar, 200.000 kali pengulangan, waktu yang diperlukan masih di bawah 1 menit. Ada hubungan antara kesetimbangan sistem dengan jumlah iterasi. Jika jumlah iterasi yang ditentukan tidak cukup besar, maka ada kemungkinan sistem belum mencapai kesetimbangan pada saat iterasi berakhir. Berhubungan juga dengan ukuran kisi. Sebagai contoh, untuk ukuran kisi 10×10, dengan jumlah iterasi 5000 kali, kesetimbangan sistem sudah tercapai, namun belum tercapai untuk ukuran kisi yang lebih besar. V. KESIMPULAN Dari hasil-hasil yang diperoleh dapat ditunjukkan bahwa lembarkerja Microsoft Excel dapat digunakan sebagai perangkat untuk komputasi numerik, simulasi dan visualisasi. Hasil simulasi cocok dengan hasil yang diperoleh dengan program-program lain, yaitu Python dan Java. Kekurangan pada lembarkerja adalah waktu komputasi yang diperlukan cukup lama jika visualisasi ditampilkan serta jumlah data/pengulangan yang besar. PUSTAKA [1] [2]

[3]

D. J. Griffiths, Introduction to Electrodynamics, third edition, Prentice Hall K. Binder and D.W. Heermann, Monte Carlo Simulation in Statistical Physics- An Introducion, 2002, fourth edition, Springer

[7] [8]

39

N. Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenbluth, Augusta H. Teller, and Edward Teller, 1953, Equation of state calculations by fast computing machines., 1953, The Journal of Chemical Physics, 21(6):1087–1092 D. Marchand, Classical Monte Carlo and the Metropolis Algorithm: Revisiting the 2D Ising Model, 2005, Department of Physics and Astronomi, University of British Columbia, Vancouver. J. Ma, 2D Ising Model Simulation, 2007, Department of Physics, UCDAVIS, [email protected] E.S. Velasco, A Java Applet to simuate a 2-D Ising system, 1998, Website: http://www2.truman.edu/~velasco/ising.html, diakses pada tanggal 8 Februari 2012 I. Argawal, Numerical Analysis of 2-D Ising Model, March 2011, physics report, University of Bonn. P.M. Hemenger and H. Weik, Determination of the Ferromagnetic Exchange Energy Constant in Ferromagnetic Ni and NixCu Films by Means of Domain Wall Width Measurements, 1965, Journal Applied Physics 36.

TANYA JAWAB Ari Setiawan, UGM ? Dalam ratlab bisa memvisualisasi dalam 3D. Bagaimana dengan Excell apakah dimungkinkan untuk mengembangkan dalam tampilan 3D. Eko S, UGM √ Excel tidak memiliki kemampuan untuk menampilkan visualisasi 3 dimensi. Namun perhitungan magnetisasi yang melibatkan spin-spin tetangga (atas, bawah, kiri, kanan, depan, dan belakang, untuk 3 dimensi) bisa dilakukan komputasi numeriknya. Jika 2 dimensi perlu 2 indeks ( i dan j, maka untuk 3 dimensi, perlu 3 indeks I, j, k array 3 dimensi. Visualisasi 3 dimensinya bisa dilakukan dengan program lain (mis pythom) jika angka? Sudah diperoleh dengan excel.


sekali magnet yang di sini. Karena

Recommend Documents