View source for N-Body Problems

'''N-body problem''' adalah masalah memprediksi gerakan individual sekumpulan objek bermassa yang berinteraksi satu sama lain melalui [[gaya gravitasi]]. Untuk dua body (''two-body problem''), terdapat solusi analitik tertutup yang diberikan oleh [[hukum Kepler]]. Namun untuk <math>N \geq 3</math>, sistem umumnya '''tidak memiliki solusi analitik''' dan harus diselesaikan secara numerik melalui simulasi.

Artikel ini menjelaskan formulasi matematika N-body problem berbasis [[hukum gravitasi Newton]] dan memberikan contoh perhitungan manual untuk kasus tiga body.

== Hukum Gravitasi Newton ==

Antara dua body bermassa <math>m_1</math> dan <math>m_2</math> yang berjarak <math>r</math>, terdapat gaya tarik gravitasi yang dirumuskan oleh [[Isaac Newton]] sebagai:

<math display="block">F = G \frac{m_1 m_2}{r^2}</math>

di mana:
* <math>G = 6{,}674 \times 10^{-11} \; \mathrm{N \cdot m^2 / kg^2}</math> adalah [[konstanta gravitasi universal]]
* <math>r</math> adalah jarak antara kedua body (meter)
* <math>F</math> adalah besar gaya (Newton)

Gaya gravitasi memiliki dua sifat penting:
# Selalu '''menarik''' (''attractive''), tidak pernah menolak.
# Bersifat '''mutual''' — body 1 menarik body 2 dengan gaya yang sama besar dan berlawanan arah dengan gaya body 2 ke body 1 ([[hukum III Newton]]).

== Formulasi Vektor ==

Formula skalar <math>F = G m_1 m_2 / r^2</math> hanya memberikan ''besar'' gaya. Untuk simulasi numerik, kita memerlukan informasi '''arah''' gaya. 

Misalkan body <math>i</math> berada pada posisi <math>\vec{r}_i = (x_i, y_i, z_i)</math> dan body <math>j</math> pada posisi <math>\vec{r}_j</math>. Maka:

'''Vektor jarak''' dari <math>i</math> ke <math>j</math>:

<math display="block">\vec{r}_{ij} = \vec{r}_j - \vec{r}_i = (x_j - x_i,\; y_j - y_i,\; z_j - z_i)</math>

'''Magnitudo jarak''':

<math display="block">r_{ij} = \| \vec{r}_{ij} \| = \sqrt{(x_j - x_i)^2 + (y_j - y_i)^2 + (z_j - z_i)^2}</math>

'''Gaya gravitasi pada body <math>i</math> yang disebabkan oleh body <math>j</math>''':

<math display="block">\vec{F}_{i \leftarrow j} = G \frac{m_i m_j}{r_{ij}^2} \cdot \hat{r}_{ij} = G \frac{m_i m_j}{r_{ij}^3} \, \vec{r}_{ij}</math>

Pembagi <math>r_{ij}^3</math> (bukan <math>r_{ij}^2</math>) muncul karena vektor satuan <math>\hat{r}_{ij} = \vec{r}_{ij} / r_{ij}</math> memerlukan pembagian tambahan dengan <math>r_{ij}</math>.

== Generalisasi ke N Body ==

Untuk sistem dengan <math>N</math> body, '''total gaya''' pada body <math>i</math> adalah jumlah vektor dari semua gaya yang dikerjakan oleh body-body lain:

<math display="block">\vec{F}_i = \sum_{\begin{smallmatrix} j=1 \\ j \neq i \end{smallmatrix}}^{N} G \frac{m_i m_j}{r_{ij}^3} \, \vec{r}_{ij}</math>

Persamaan ini adalah inti dari setiap implementasi N-body. Kompleksitas perhitungannya adalah <math>O(N^2)</math>, karena untuk setiap body harus dilakukan iterasi terhadap <math>N-1</math> body lain.

== Hukum II Newton dan Integrasi Numerik ==

Setelah total gaya pada body <math>i</math> diketahui, '''percepatan''' diperoleh dari [[hukum II Newton]]:

<math display="block">\vec{F}_i = m_i \vec{a}_i \quad \Longrightarrow \quad \vec{a}_i = \frac{\vec{F}_i}{m_i}</math>

Untuk memperbarui kecepatan dan posisi seiring waktu, digunakan metode [[integrasi numerik]]. Metode paling sederhana adalah [[metode Euler|Euler semi-implisit]]:

<math display="block">\vec{v}_i(t + \Delta t) = \vec{v}_i(t) + \vec{a}_i(t) \cdot \Delta t</math>

<math display="block">\vec{r}_i(t + \Delta t) = \vec{r}_i(t) + \vec{v}_i(t + \Delta t) \cdot \Delta t</math>

Metode yang lebih akurat seperti [[Leapfrog integration|Leapfrog]] atau [[metode Runge-Kutta|Runge-Kutta]] umum digunakan dalam simulasi production karena memiliki sifat konservasi energi yang lebih baik (lihat bagian {{Section link||Drift Energi}}).

== Contoh Perhitungan: 3-Body ==

Untuk memberikan ilustrasi konkret, berikut adalah perhitungan satu ''time step'' untuk sistem 3 body di bidang dua dimensi (sumbu <math>z = 0</math>).

=== Kondisi Awal ===

{| class="wikitable" style="text-align: center;"
! Body !! Massa (kg) !! Posisi (m) !! Kecepatan (m/s)
|-
| A || <math>1 \times 10^{10}</math> || <math>(0,\; 0)</math> || <math>(0,\; 0)</math>
|-
| B || <math>1 \times 10^{10}</math> || <math>(1,\; 0)</math> || <math>(0,\; 0)</math>
|-
| C || <math>1 \times 10^{10}</math> || <math>(0,\; 1)</math> || <math>(0,\; 0)</math>
|}

Ketiga body berada dalam keadaan diam, membentuk segitiga siku-siku. Gunakan <math>\Delta t = 1</math> detik dan abaikan [[#Catatan Implementasi|softening]] (<math>\epsilon = 0</math>) untuk kebersihan perhitungan.

=== Langkah 1: Vektor Jarak Antar Pasangan ===

'''Pasangan A–B:'''

<math display="block">\vec{r}_{AB} = (1 - 0,\; 0 - 0) = (1,\; 0), \quad r_{AB} = 1, \quad r_{AB}^3 = 1</math>

'''Pasangan A–C:'''

<math display="block">\vec{r}_{AC} = (0 - 0,\; 1 - 0) = (0,\; 1), \quad r_{AC} = 1, \quad r_{AC}^3 = 1</math>

'''Pasangan B–C:'''

<math display="block">\vec{r}_{BC} = (0 - 1,\; 1 - 0) = (-1,\; 1), \quad r_{BC} = \sqrt{2} \approx 1{,}414, \quad r_{BC}^3 = 2\sqrt{2} \approx 2{,}828</math>

=== Langkah 2: Perhitungan Gaya ===

Karena semua body bermassa sama, faktor <math>G \cdot m_i \cdot m_j</math> bernilai konstan untuk semua pasangan:

<math display="block">K \equiv G \cdot m_A \cdot m_B = 6{,}674 \times 10^{-11} \times 10^{10} \times 10^{10} = 6{,}674 \times 10^{9} \; \mathrm{N \cdot m^2}</math>

'''Gaya pada body A''' (kontribusi dari B dan C):

<math display="block">\vec{F}_A = \frac{K}{r_{AB}^3} \vec{r}_{AB} + \frac{K}{r_{AC}^3} \vec{r}_{AC} = \frac{K}{1}(1,\; 0) + \frac{K}{1}(0,\; 1)</math>

<math display="block">\vec{F}_A = (6{,}674 \times 10^9,\; 6{,}674 \times 10^9) \; \mathrm{N}</math>

''Interpretasi: A ditarik ke kanan (oleh B) dan ke atas (oleh C).''

'''Gaya pada body B''' (kontribusi dari A dan C):

Vektor dari B ke A adalah <math>(-1,\; 0)</math>, dan vektor dari B ke C adalah <math>(-1,\; 1)</math>:

<math display="block">\vec{F}_B = \frac{K}{1}(-1,\; 0) + \frac{K}{2\sqrt{2}}(-1,\; 1)</math>

<math display="block">\vec{F}_B = (-6{,}674 \times 10^9,\; 0) + (-2{,}360 \times 10^9,\; 2{,}360 \times 10^9)</math>

<math display="block">\vec{F}_B = (-9{,}034 \times 10^9,\; 2{,}360 \times 10^9) \; \mathrm{N}</math>

''Interpretasi: B ditarik kuat ke kiri (oleh A yang dekat) dan sedikit ke atas (oleh C yang lebih jauh).''

'''Gaya pada body C''' (kontribusi dari A dan B):

Vektor dari C ke A adalah <math>(0,\; -1)</math>, dan vektor dari C ke B adalah <math>(1,\; -1)</math>:

<math display="block">\vec{F}_C = \frac{K}{1}(0,\; -1) + \frac{K}{2\sqrt{2}}(1,\; -1)</math>

<math display="block">\vec{F}_C = (2{,}360 \times 10^9,\; -9{,}034 \times 10^9) \; \mathrm{N}</math>

''Interpretasi: C ditarik kuat ke bawah (oleh A) dan sedikit ke kanan (oleh B).''

=== Langkah 3: Verifikasi Hukum III Newton ===

Untuk sistem terisolasi, jumlah total semua gaya internal harus nol:

<math display="block">\vec{F}_A + \vec{F}_B + \vec{F}_C \stackrel{?}{=} \vec{0}</math>

Komponen <math>x</math>:

<math display="block">6{,}674 + (-9{,}034) + 2{,}360 = 0 \quad \checkmark</math>

Komponen <math>y</math>:

<math display="block">6{,}674 + 2{,}360 + (-9{,}034) = 0 \quad \checkmark</math>

Sanity check ini sangat penting dalam implementasi: jika jumlah total gaya tidak nol (di luar ''floating-point error''), kemungkinan besar ada kesalahan dalam perhitungan.

=== Langkah 4: Percepatan ===

Bagi masing-masing gaya dengan massa body (semua <math>10^{10}</math> kg):

{| class="wikitable" style="text-align: center;"
! Body !! <math>a_x \; (\mathrm{m/s^2})</math> !! <math>a_y \; (\mathrm{m/s^2})</math>
|-
| A || <math>0{,}6674</math> || <math>0{,}6674</math>
|-
| B || <math>-0{,}9034</math> || <math>0{,}2360</math>
|-
| C || <math>0{,}2360</math> || <math>-0{,}9034</math>
|}

=== Langkah 5: Update Kecepatan dan Posisi ===

Dengan <math>\Delta t = 1</math> detik dan kecepatan awal nol, metode Euler semi-implisit menghasilkan:

{| class="wikitable" style="text-align: center;"
! Body !! Kecepatan baru (m/s) !! Posisi baru (m)
|-
| A || <math>(0{,}6674,\; 0{,}6674)</math> || <math>(0{,}6674,\; 0{,}6674)</math>
|-
| B || <math>(-0{,}9034,\; 0{,}2360)</math> || <math>(0{,}0966,\; 0{,}2360)</math>
|-
| C || <math>(0{,}2360,\; -0{,}9034)</math> || <math>(0{,}2360,\; 0{,}0966)</math>
|}

Satu ''time step'' telah selesai. Untuk ''time step'' berikutnya, ulangi prosedur dari [[#Langkah 1: Vektor Jarak Antar Pasangan|Langkah 1]] dengan posisi dan kecepatan baru ini.

== Implementasi dalam C++ ==

Berikut implementasi minimal dari satu ''time step'' untuk verifikasi perhitungan di atas:

<syntaxhighlight lang="cpp">
#include <iostream>
#include <vector>
#include <cmath>

constexpr double G = 6.674e-11;

int main() {
    std::vector<double> x = {0.0, 1.0, 0.0};
    std::vector<double> y = {0.0, 0.0, 1.0};
    std::vector<double> vx(3, 0), vy(3, 0);
    std::vector<double> mass = {1e10, 1e10, 1e10};
    int n = 3;
    double dt = 1.0;

    std::vector<double> fx(n, 0), fy(n, 0);
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            if (i == j) continue;
            double dx = x[j] - x[i];
            double dy = y[j] - y[i];
            double r2 = dx*dx + dy*dy;
            double r3 = r2 * std::sqrt(r2);
            double s = G * mass[i] * mass[j] / r3;
            fx[i] += s * dx;
            fy[i] += s * dy;
        }
    }

    // Semi-implicit Euler integration
    for (int i = 0; i < n; ++i) {
        vx[i] += fx[i] / mass[i] * dt;
        vy[i] += fy[i] / mass[i] * dt;
        x[i]  += vx[i] * dt;
        y[i]  += vy[i] * dt;
    }
}
</syntaxhighlight>

Eksekusi program ini menghasilkan nilai posisi dan kecepatan yang konsisten dengan perhitungan manual hingga presisi ''floating-point''.

== Ringkasan Algoritma ==

Untuk setiap ''time step'', siklus perhitungan adalah:

# '''Reset gaya''': <math>\vec{F}_i = \vec{0}</math> untuk semua <math>i</math>.
# '''Loop atas pasangan body''': untuk setiap <math>i, j</math> dengan <math>j \neq i</math>, akumulasi:
#: <math>\vec{F}_i \mathrel{+}= G \frac{m_i m_j}{r_{ij}^3} \vec{r}_{ij}</math>
# '''Hitung percepatan''': <math>\vec{a}_i = \vec{F}_i / m_i</math>.
# '''Update kecepatan''': <math>\vec{v}_i \mathrel{+}= \vec{a}_i \cdot \Delta t</math>.
# '''Update posisi''': <math>\vec{r}_i \mathrel{+}= \vec{v}_i \cdot \Delta t</math>.

Siklus ini diulang untuk ribuan hingga jutaan ''time step'', tergantung pada durasi simulasi yang diinginkan.

== Catatan Implementasi ==

=== Kompleksitas dan Paralelisasi ===

[[Kompleksitas algoritma|Kompleksitas]] perhitungan langsung (''direct sum'') adalah <math>O(N^2)</math> per ''time step''. Untuk <math>N = 10^3</math>, terdapat <math>10^6</math> perhitungan pasangan per ''step''; untuk <math>N = 10^6</math>, terdapat <math>10^{12}</math> perhitungan per ''step''. Kebutuhan komputasi yang masif ini memotivasi penggunaan [[komputasi paralel]] menggunakan teknologi seperti [[CUDA]], [[OpenCL]], dan [[Message Passing Interface|MPI]].

Algoritma alternatif dengan kompleksitas lebih rendah meliputi:
* [[Algoritma Barnes-Hut]] — kompleksitas <math>O(N \log N)</math>, menggunakan struktur [[octree]] untuk meng-aproksimasi gaya dari kelompok body yang jauh.
* [[Fast Multipole Method]] — kompleksitas <math>O(N)</math>, menggunakan ekspansi multipolar.

=== Softening Parameter ===

Jika dua body mendekati jarak nol, faktor <math>1/r^3</math> menyebabkan gaya menuju tak hingga (''singularitas''), yang menyebabkan simulasi tidak stabil. Untuk mengatasinya, digunakan [[Plummer softening|softening parameter]] <math>\epsilon</math>:

<math display="block">r_{ij}^{\,\text{soft}} = \sqrt{r_{ij}^2 + \epsilon^2}</math>

Secara fisik, ini setara dengan memodelkan body sebagai distribusi massa terdistribusi (bukan titik massa), sehingga interaksi pada jarak sangat dekat tetap finit.

=== Drift Energi ===

Untuk sistem gravitasi tertutup, energi total <math>E = E_\text{kinetik} + E_\text{potensial}</math> harus kekal:

<math display="block">E_\text{tot} = \sum_i \frac{1}{2} m_i v_i^2 - \sum_{i<j} \frac{G m_i m_j}{r_{ij}} = \text{konstan}</math>

Namun integrasi numerik selalu memperkenalkan kesalahan ''truncation error'' kecil di setiap ''step''. Akumulasi kesalahan ini menyebabkan '''drift energi''' — perubahan bertahap dalam energi total simulasi.

Metode Euler memiliki ''drift'' yang tumbuh linear terhadap waktu. Metode [[Leapfrog integration|Leapfrog]] yang bersifat [[Symplectic integrator|simplektik]] memiliki ''drift'' yang ''bounded'' (berosilasi tetapi tidak tumbuh), sehingga lebih disukai untuk simulasi jangka panjang.

=== Optimisasi Hukum III Newton ===

Karena <math>\vec{F}_{i \leftarrow j} = -\vec{F}_{j \leftarrow i}</math>, secara teoritis kita dapat menghemat 50% perhitungan dengan iterasi <math>j > i</math> saja dan mengakumulasi gaya ke kedua body sekaligus. Namun optimisasi ini menyulitkan paralelisasi pada [[GPU]] karena menyebabkan ''write conflict'' antar [[thread (komputasi)|thread]]. Pada implementasi paralel, perhitungan <math>N^2</math> penuh seringkali lebih efisien.

== Lihat Pula ==

* [[Hukum gravitasi Newton]]
* [[Two-body problem]]
* [[Three-body problem]]
* [[Algoritma Barnes-Hut]]
* [[Simulasi numerik]]
* [[Komputasi paralel]]

== Referensi ==

* Aarseth, S. J. (2003). ''Gravitational N-Body Simulations: Tools and Algorithms''. Cambridge University Press.
* Hairer, E., Lubich, C., & Wanner, G. (2006). ''Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations''. Springer.
* Nyland, L., Harris, M., & Prins, J. (2007). [https://developer.nvidia.com/gpugems/gpugems3/part-v-physics-simulation/chapter-31-fast-n-body-simulation-cuda "Fast N-Body Simulation with CUDA"]. ''GPU Gems 3'', NVIDIA Corporation.