PENERAPAN RANTAI-AR (ALGORITHMIC REGULARIZATIONS) DALAM MASALAH N-BENDA

158

Prosiding Pertemuan Ilmiah XXIV HFI Jateng & DIY, Semarang 10 April 2010 hal. 158-167

PENERAPAN RANTAI-AR (ALGORITHMIC REGULARIZATIONS) DALAM MASALAH N-BENDA Nugroho Adi P Kelompok Riset Kosmologi, Astrofsika, dan Fisika Matematik (KAM) Jurusan Fisika FMIPA UGM Jurusan Fisika Universitas Negeri Malang e-mail:[email protected] M. F. Rosyid Kelompok Riset Kosmologi, Astrofisika, dan Fisika Matematik (KAM) Jurusan Fisika FMIPA UGM, Yogyakarta e-mail:[email protected]

INTISARI Dipelajari simulasi n-benda akibat interaksi gravitasi. Dikaji kemungkinan yang terjadi jika menggunakan metode numerik biasa. Ditunjukkan sebagian besar integrasi kehilangan ketelitian pada beberapa situasi dikarenakan faktor 1⁄r2 pada gaya gravitasi. Dibicarakan berbagai metode untuk mengatasi hal tersebut. Dalam makalah ini hanya dibicarakan metode rantai-AR yang merupakan algoritma yang memiliki hasil tanpa menggunakan transformasi koordinat. Kata kunci: rantai-AR, n-benda

I. PENDAHULUAN Pada simulasi N-benda, umumnya interaksi yang paling kuat disebabkan karena mendekatnya dua benda. Kebanyakan integrasi numerik kehilangan ketelitian untuk keadaan tertentu karena singularitas 1⁄r2 dari gaya bersama pada dua benda [Mikkola, 2008]. Akibatnya, hal yang paling penting dari regularisasi algoritme adalah bahwa algoritma tersebut menangani masalah dua benda yang mengalami gangguan. Pada dasarnya ada dua jenis metode yang tersedia yaitu transformasi waktu dan koordinat dan algoritma yang menghasilkan hasil biasa tanpa transformasi koordinat. Artikel ini membahas tentang Algorithm Regularization dan secara khusus membahas tentang rantai-AR yang memberi penekanan pada pentingnya struktur rantai untuk mencegah penjumlahan dengan bilangan yang besar. Hal ini penting karena penjumlahan dengan bilangan yang besar akan megakibatkan peniadaan bilangan-bilangan kecil yang penting dan mengakibatkan galat yang tidak dapat diperbaiki. II. METODE PENELITIAN II.1. Algorithmic Regularization Algorithmic regularization tidak menggunakan transformasi koordinat namun hanya menggunakan transformasi waktu dan algoritme yang sesuai yang menghasilkan hasil biasa dengan memperhatikan singularitas gaya. II.1.1 Logarithmic Hamiltonian (LogH) untuk Dua Benda Hamiltonian di ruang-fase yang diperluas adalah

B adalah momentum dari waktu (yang merupakan koordinat: fungsi

dapat digunakan sebagai Hamiltonian. Misal pada kasus gerak dua-benda H = p2/2 -M/r memberikan

ISSN 0853 - 0823

. Jika B(0) = - H(0), maka

Nugroho Adi P, dkk/ Penerapan Rantai -AR (Algorithmic Regularizations) dalam Masalah -N Benda

159

Transformasi waktunya adalah

B tetap konstan, B = -(p2/2 -M/r). Variabel baru s adalah

merupakan kuantitas yang sebanding dengan peningkatan anomali eksentris. Dapat dilihat bahwa sigularitas saat r→0 tidak hilang II.1.2. Time-Transformed Leapfrog (TTL) Tinjau sistem umum

Diperkenalkan transformasi waktu di mana Ω(r) > 0 sebarang. jika W = Ω, maka dapat ditulis

Jika W didapat dari persamaan differensial

alih-alih W = Ω , kita akan mendapatkan

Leapfrog sendiri, dalam banyak kasus, tidak begitu akurat [Mikkola, 2008]. Keakuratan dapat diperoleh dengan menggunakan algoritme leapfrog orde tinggi. Alternatif lain adalah, dengan menggunakan metode ekstrapolasi II.2. Rantai-AR Perlu ditonjolkan keunggulan struktur rantai di algorithmic regularization. Alasannya adalah galat pembulatan. Jika menggunakan koordinat pusat massa, koordinat relatif dari jarak pasangan terdekat, sebagai penyebut, akan memiliki nilai yang besar dan akan meniadakan/mengabaikan beberapa nilai penting, yang mengakibatkan galat yang tidak bisa diperbaiki [Mikkola, 2008]. II.2.1. Menemukan dan Memperbaharui Rantai Kita mulai dengan menemukan vektor interpartikel terpendek untuk bagian pertama rantai. Berikutnya kita mencari partikel yang dekat dengan satu ujung atau ujung yang lainnya yang sudah merupakan bagian dari rantai. Partikel ini ditambahkan ke ujung yang terdekat dari rantai yang telah ada. Ini diulang hingga semua partikel termasuk di dalam rantai. Setelah setiap langkah integrasi, kita cek apakah rantai perlu diperbaharui. Gambar 1 menggambarkan kasus dari sebuah rantai 10-partikel. Untuk menghindari beberapa masalah pembulatan, akan sangat menguntungkan untuk membawa transformasi dari vektor-vektor rantai yang lama Xk ke ke yang baru secara langsung dengan menyatakan vektor-vektor rantai yang baru sebagai jumlahan dari yang lama. ISSN 0853 - 0823

160


Misal nama-nama fisis dari partikel-partikel 1,...., N adalah I1, I2, ..., IN dan misal kita gunakan notasi Ikold dan Iknew untuk nama di rantai lama dan baru. Maka kita dapat menuliskan

Jadi kita membutuhkan hubungan antara indeks lama dan baru untuk menyatakan vektor-vektor rantai baru X dalam suku lama. Telah ditemukan bahwa jika k0 dan k1 adalah dua indeks sedemikian sehingga dan maka

Gambar 1. Penamaan Rantai. dimana Bµυ = +1 jika (k1 > υ dan k0 ≤ υ ) dan Bµυ = -1 jika (k1≤ υ dan k0 > υ ), selain nilai tersebut Bµυ = 0. II.3. Transformasi Setelah memilih rantai, dan menamai partikel sebagai ' 1,2, ... ,N sepanjang rantai, kita dapat mengevaluasi nilai-nilai awal untuk vektor-vektor rantai dan kecepatan sebagai

dimana vk =

. Di waktu yang sama kita dapat memeriksa pusat massa

Transformasi balik ke r, v dapat dilakukan dengan jumlahan sederhana

diikuti dengan reduksi pusat masa

ISSN 0853 - 0823


161

Bagaimanapun, kita tidak harus mereduksi koordinat ke sistem pusat massa karena percepatan hanya bergantung pada perbedaan. II.4. Persamaan Gerak dan Leapfrog Persamaan gerak membaca

dimana percepatan Ak, dengan gaya luar yang mungkin adalah fk, adalah

dan, untuk j < k

Untuk k > j kita gunakan fakta bahwa rjk = -rkj . Penggunaan Xj dan Xj + Xj+1 mereduksi efek pembulatan secara signifikan. Secara lebih umum, kita dapat juga menggunakan

tetapi untuk banyak benda, lebih cepat dengan menggunakan (30) dan alternatif yang disebut terakhir terlihat tidak meningkatkan hasilnya Energi kinetik adalah

dan energi potensial

yang diperiksa bersama dengan percepatan menurut (30). Diperkenalkan sebuah fungsi transformasi lebih jauh

di mana Ωij adalah koefisien (dibahas di subbab di bawah). Didefinisikan dua transformasi waktu

di mana α, β, γ dan adalah konstanta yang dapat disesuaikan, B = U − T adalah energi ikat N-benda dan ω didefinisikan oleh persamaan differensial

dan nilai awal ω(0) = Ω (0). Energi ikat B berubah menurut

ISSN 0853 - 0823

162


Persamaan gerak yang dapat digunakan untuk membangun leapfrog yang menyediakan algorihmic regularization, untuk waktu dan koordinat yang diberikan, adalah

dan untuk kecepatan B dan ω

Leapfrog untuk vektor rantai Xk dan Vk dapat ditulis paling mudah dalam bentuk dua pemetaan

di mana < vk > adalah rata-rata v awal dan akhir. Perhatikan juga bahwa perlu untuk memeriksa kecepatan individual vk, karena, jika tidak, B' dan ω' akan sulit dihitung. Satu langkah leapfrog dapat ditulis sebagai

dan sekuens yang lebih panjang dari n-langkah

Rumus ini akan digunakan dengan metode ekstrapolasi saat memproses sebuah interval waktu total dari panjang nh. II.5. Transformasi Waktu Alternatif

ISSN 0853 - 0823


163

II.6. Algoritma Dasar untuk Metode Extrapolasi II.6.1. Leapfrog

III. HASIL DAN PEMBAHASAN Hasil yang diharapkan dari rantai-AR adalah pengurangan nilai galat pada penghitungan numerik. Gambar 2 menunjukkan galat relatif dari energi yang berayun di sekitar nilai 6×10-13 dan memiliki kecenderungan untuk membesar. Hal tersebut kemungkinan disebabkan karena solusi pada langkah ganjil hasil metode leapfrog memiliki kecenderungan menjauh dari solusi pada langkah genap, sehingga sudah wajar untuk menghentikan integrasi tiap 20 langkah dan mengawali kembali dengan metode Euler maju [Shampine, ]. Namun cara tersebut mengakibatkan metode ini secara keseluruhan, dengan penghentian tiap 20 langkah, tidak stabil. Diharapkan dengan menggunakan metode Runge-Kutta, galat tersebut dapat diperkecil baik nilai galat dan ayunan galat tersebut.

Gambar 2. Galat di masalah 10-benda yang diintegralkan menggunakan kode rantai-AR. Kurva paling atas adalah galat relatif dari energi. IV. KESIMPULAN DAN SARAN Rantai-AR merupakan alternatif yang baik dalam simulasi N-benda. Bagaimanapun Rantai-AR adalah metode untuk menentukan vektor jarak antar partikel pada simulasi N-benda, sedangkan pada simulasi N-benda juga melibatkan persamaan gerak yang pada artikel ini diselesaikan dengan metode leapfrog yang dalam banyak kasus kurang akurat. Untuk itu diperlukan metode leapfrog dengan orde yang lebih tinggi atau menggunakan metode lain yang akurat, misal Runge-Kutta.

ISSN 0853 - 0823

164


Program Untuk Rantai-AR program rantaiAR implicit none real m(10) ,f(10) ,x(10) ,y(10) ,z(10) , rx (10) , ry (10) , rz (10) ,& &vx (10) , vy (10) , vz (10) , rvx (10) , rvy (10) , rvz (10) , ax (10) , ay (10) , az (10) ,& &e(10) ,u(10) , rpmx ,rpmy ,vpmx ,vpmy ,& &t,tmax ,h,dummy ,mtot ,rjk ,ek ,ep ,ep0 ,ek0 ,ohm , ohmjk integer n,i,j,k t=0 h =0.01 tmax =2 mtot =0 ek =0 ep =0 ohm =0 rpmx =0 rpmy =0 vpmx =0 vpmy =0 read (* ,*)n do i=1,n read (* ,*)m(i),x(i),y(i),vx(i),vy(i) ax(i)=0 f(i)=0 end do ! nilai awal print *,' Posisi Awal ' do i=1,n print *,i,m(i),x(i),y(i),vx(i),vy(i) end do ! buat rantai 1 continue ! catat posisi dan kecepatan do i=1,n -1 rx(i)=x(i+1) -x(i) ry(i)=y(i+1) -y(i) rvx(i)= vx(i+1) - vx(i) rvy(i)= vy(i+1) - vy(i) end do ! periksa pusat massa do i=1,n mtot = mtot +m(i) end do do i=1,n rpmx = rpmx +m(i)*x(i)/ mtot rpmy = rpmy +m(i)*y(i)/ mtot vpmx = vpmx +m(i)* vx(i)/ mtot vpmy = vpmy +m(i)* vy(i)/ mtot end do ! Xdot_ {k}= V_{k} do k=1,n -1 rx(k)=h*rvx (k) end do ! Vdot_k =A_{k+1} - A_{k} do k=1,n -1 rvx(k)=h*( ax(k+1) - ax(k)) end do ISSN 0853 - 0823


165

! tentukan r do k=1,n -1 do j=1,n -1 if(j.ne.k) then if(j .lt. k) then if (k .gt.j +2) then rjk=x(k)-x(j) !end if else if(k==j+1) then rjk=rx(j) !end if else if(k==j+2) then rjk=rx(j)+ rx(j+1) end if ! hitung energi potensial else if (j .gt.k+2) then rjk=x(j)-x(k) !end if else if(j==k+1) then rjk=rx(k) !end if else if(j==k+2) then rjk=rx(k)+ rx(k+1) end if ! hitung energi potensial end if ep=ep+m(j)*m(k)/ abs (rjk ) ohm=ohm+ ohmjk /abs(rjk) !A_{k}=\ sum_ {j\ne k} ax(k)= -( ax(k)+m(j)* rjk /( abs (rjk ))**3+ f(k)) end if end do end do ! hitung energi kinetik do k=1,n ek=ek +0.5* m(k)* vx(k )**2 end do !tt =1/( alfa *( ek+b)+ beta * omega + gama ) ! Transformasi balik x (1)=0 vx (1)=0 do k=1,n -1 x(k +1)= x(k)+ rx(k) vx(k +1)= vx(k)+ rvx (k) end do t=t+h print *,x(2) , vx (2) , ep+ek if(t.lt. tmax ) goto 1 print *,' Posisi akhir ' do i=1,n print *,i,m(i),x(i),y(i),vx(i),vy(i) end do print *,' Selesai ' end program rantaiAR Untuk Time-Transformed Leapfrog implicit real *8 (a-h,m,o-z) ISSN 0853 - 0823

166


read (5 ,*)h,tincr ,tmx , mass ! read step ,tincr , maxtime , mass read (5 ,*)x,y,z,vx ,vy ,vz ! read initial coords / vels

!

! 1

!

2

tnext =0 initializations t=0 r= sqrt (x*x+y*y+z*z) ! distance vv=vx*vx+vy*vy+vz*vz !v- square E0=vv /2- mass /r W= mass /r Integration of two - body motion continue dt=h/W/2 ! time increment t=t+dt x=x+dt*vx y=y+dt*vy z=z+dt*vz dtc=h/(x*x+y*y+z*z) dw= -(x*vx+y*vy+z*vz )* dtc /2 vx=vx -x*dtc vy=vy -y*dtc vz=vz -z*dtc W=W+dw -(x*vx+y*vy+z*vz )* dtc /2 dt=h/W/2 ! new time increment t=t+dt ! this has an O(h^3) error x=x+dt*vx y=y+dt*vy z=z+dt*vz diagnostics if(t.lt. tnext ) goto 1 tnext = tnext + tincr r= sqrt (x*x+y*y+z*z) err=-E0 +( vx*vx+vy*vy+vz*vz )/2 - mass /r write (6 ,2)t,x,y,z,err*r ,W*r- mass ! time , coords & error if(t.lt.Tmx ) goto 1 format (1x ,1p ,10 g12 .4) end

Untuk leapfrog ! Leapfrog code for a harmonic oscillator: ! ----------------------------------------------implicit real *8 (a-h,o-z) x=1 p=0 h =0.01 d0 E0 =(p*p+x*x)/2 t=0 1 continue x=x+h/2*p ! this is p=p-h*x ! a leapfrog x=x+h/2*p ! step t=t+h ! diagnostics E=(p*p+x*x)/2 write (6 ,*)t,x,p,E-E0 if(t.lt .100.) goto 1 ! max time =100. End

ISSN 0853 - 0823


167

Untuk turunan leapfrog ! Differentiated leapfrog for harmonic oscillator !---------------------------------------------implicit real *8 (a-h,o-z) x=1 dx =1 ! var p=0 dp =0 ! var E0 =(p*p+x*x)/2 dE0=p*dp+x*dx ! var t=0 h =0.01 d0 ! stepsize 1 continue x=x+h/2*p ! this is dx=dx+h/2* dp ! var p=p-h*x ! a leapfrog dp=dp -h*dx ! var x=x+h/2*p ! step dx=dx+h/2* dp ! var t=t+h ! diagnostics E=(p*p+x*x)/2 dE=p*dp+x*dx ! var ( this should be constant !) write (6 ,*)t,x,p,E-E0 ,dE -dE0 if(t.lt .100.) goto 1 ! max time =100. end [Kosasih, 2006], [Chapra and Canale, 1994], [Koonin and Meredith, 1990] V. DAFTAR PUSTAKA Chapra, S. C. and Canale, R. P. (1994). Metode Numerik. 1st. Penerbit Erlangga, Jakarta, 2nd edition. Koonin, S. E. and Meredith, D. C. (1990). COMPUTATIONAL PHYSICS, Fortran Edition. Westview. Kosasih, P. B. (2006). Komputasi Numerik Teori dan Aplikasi.Penerbit Andi, Yogyakarta. Mikkola, S. (2008). Regular algorithm for the few-body problem. In The Cambridge N-Body Lectures, pages 31-58. Springer. Shampine, L. Stability of the Leapfrog/Midpoint Method. Mathematics Departement Shouthern Methodist University.

ISSN 0853 - 0823

PENERAPAN RANTAI-AR (ALGORITHMIC REGULARIZATIONS) DALAM MASALAH N-BENDA

Recommend Documents