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 $N\ge 3$, 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 $m_1$ dan $m_2$ yang berjarak $r$, terdapat gaya tarik gravitasi yang dirumuskan oleh Isaac Newton sebagai:

$$F=G\frac{m_1{m}_2}{r^2}$$

di mana:

• $G=6,674\times 10^{-11}\mathrm{N}\mathrm{\cdot }{\mathrm{m}}^{\mathrm{2}}/\mathrm{k}{\mathrm{g}}^{\mathrm{2}}$ adalah konstanta gravitasi universal
• $r$ adalah jarak antara kedua body (meter)
• $F$ adalah besar gaya (Newton)

Gaya gravitasi memiliki dua sifat penting:

1. Selalu menarik (attractive), tidak pernah menolak.
2. 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 $F=G{m}_1{m}_2/r^2$ hanya memberikan besar gaya. Untuk simulasi numerik, kita memerlukan informasi arah gaya.

Misalkan body $i$ berada pada posisi ${\overrightarrow{r}}_i=(x_i,y_i,z_i)$ dan body $j$ pada posisi ${\overrightarrow{r}}_j$. Maka:

Vektor jarak dari $i$ ke $j$:

$${\overrightarrow{r}}_{ij}={\overrightarrow{r}}_j-{\overrightarrow{r}}_i=(x_j-x_i,y_j-y_i,z_j-z_i)$$

Magnitudo jarak:

$$r_{ij}=\Vert {\overrightarrow{r}}_{ij}\Vert =\sqrt{(x_j-x_i{)}^2+(y_j-y_i{)}^2+(z_j-z_i{)}^2}$$

Gaya gravitasi pada body $i$ yang disebabkan oleh body $j$:

$${\overrightarrow{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}{\overrightarrow{r}}_{ij}$$

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

Generalisasi ke N Body

Untuk sistem dengan $N$ body, total gaya pada body $i$ adalah jumlah vektor dari semua gaya yang dikerjakan oleh body-body lain:

Failed to parse (unknown function "\substack"): {\displaystyle \vec{F}_i = \sum_{\substack{j=1 \\ j \neq i}}^{N} G \frac{m_i m_j}{r_{ij}^3} \, \vec{r}_{ij}}

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

Hukum II Newton dan Integrasi Numerik

Setelah total gaya pada body $i$ diketahui, percepatan diperoleh dari hukum II Newton:

$${\overrightarrow{F}}_i=m_i{\overrightarrow{a}}_i⟹{\overrightarrow{a}}_i=\frac{{\overrightarrow{F}}_i}{m_i}$$

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

$${\overrightarrow{v}}_i(t+\mathrm{\Delta }t)={\overrightarrow{v}}_i(t)+{\overrightarrow{a}}_i(t)\cdot \mathrm{\Delta }t$$

$${\overrightarrow{r}}_i(t+\mathrm{\Delta }t)={\overrightarrow{r}}_i(t)+{\overrightarrow{v}}_i(t+\mathrm{\Delta }t)\cdot \mathrm{\Delta }t$$

Metode yang lebih akurat seperti Leapfrog atau Runge-Kutta umum digunakan dalam simulasi production karena memiliki sifat konservasi energi yang lebih baik (lihat bagian Template:Section link).

Contoh Perhitungan: 3-Body

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

Kondisi Awal

Body	Massa (kg)	Posisi (m)	Kecepatan (m/s)
A	$1\times 10^{10}$	$(0,0)$	$(0,0)$
B	$1\times 10^{10}$	$(1,0)$	$(0,0)$
C	$1\times 10^{10}$	$(0,1)$	$(0,0)$

Ketiga body berada dalam keadaan diam, membentuk segitiga siku-siku. Gunakan $\mathrm{\Delta }t=1$ detik dan abaikan softening ($ϵ=0$) untuk kebersihan perhitungan.

Langkah 1: Vektor Jarak Antar Pasangan

Pasangan A–B:

$${\overrightarrow{r}}_{AB}=(1-0,0-0)=(1,0),r_{AB}=1,r_{AB}^3=1$$

Pasangan A–C:

$${\overrightarrow{r}}_{AC}=(0-0,1-0)=(0,1),r_{AC}=1,r_{AC}^3=1$$

Pasangan B–C:

$${\overrightarrow{r}}_{BC}=(0-1,1-0)=(-1,1),r_{BC}=\sqrt{2}\approx 1,414,r_{BC}^3=2\sqrt{2}\approx 2,828$$

Langkah 2: Perhitungan Gaya

Karena semua body bermassa sama, faktor $G\cdot m_i\cdot m_j$ bernilai konstan untuk semua pasangan:

$$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}\mathrm{\cdot }{\mathrm{m}}^{\mathrm{2}}$$

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

$${\overrightarrow{F}}_A=\frac{K}{r_{AB}^3}{\overrightarrow{r}}_{AB}+\frac{K}{r_{AC}^3}{\overrightarrow{r}}_{AC}=\frac{K}{1}(1,0)+\frac{K}{1}(0,1)$$

$${\overrightarrow{F}}_A=(6,674\times 10^9,6,674\times 10^9)\mathrm{N}$$

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 $(-1,0)$, dan vektor dari B ke C adalah $(-1,1)$:

$${\overrightarrow{F}}_B=\frac{K}{1}(-1,0)+\frac{K}{2\sqrt{2}}(-1,1)$$

$${\overrightarrow{F}}_B=(-6,674\times 10^9,0)+(-2,360\times 10^9,2,360\times 10^9)$$

$${\overrightarrow{F}}_B=(-9,034\times 10^9,2,360\times 10^9)\mathrm{N}$$

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 $(0,-1)$, dan vektor dari C ke B adalah $(1,-1)$:

$${\overrightarrow{F}}_C=\frac{K}{1}(0,-1)+\frac{K}{2\sqrt{2}}(1,-1)$$

$${\overrightarrow{F}}_C=(2,360\times 10^9,-9,034\times 10^9)\mathrm{N}$$

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:

$${\overrightarrow{F}}_A+{\overrightarrow{F}}_B+{\overrightarrow{F}}_C\stackrel{?}{=}\overrightarrow{0}$$

Komponen $x$:

$$6,674+(-9,034)+2,360=0✓$$

Komponen $y$:

$$6,674+2,360+(-9,034)=0✓$$

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 $10^{10}$ kg):

Body	$a_x(\mathrm{m}/{\mathrm{s}}^{\mathrm{2}})$	$a_y(\mathrm{m}/{\mathrm{s}}^{\mathrm{2}})$
A	$0,6674$	$0,6674$
B	$-0,9034$	$0,2360$
C	$0,2360$	$-0,9034$

Langkah 5: Update Kecepatan dan Posisi

Dengan $\mathrm{\Delta }t=1$ detik dan kecepatan awal nol, metode Euler semi-implisit menghasilkan:

Body	Kecepatan baru (m/s)	Posisi baru (m)
A	$(0,6674,0,6674)$	$(0,6674,0,6674)$
B	$(-0,9034,0,2360)$	$(0,0966,0,2360)$
C	$(0,2360,-0,9034)$	$(0,2360,0,0966)$

Satu time step telah selesai. Untuk time step berikutnya, ulangi prosedur dari Langkah 1 dengan posisi dan kecepatan baru ini.

Implementasi dalam C++

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

#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;
    }
}

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:

1. Reset gaya: ${\overrightarrow{F}}_i=\overrightarrow{0}$ untuk semua $i$.
2. Loop atas pasangan body: untuk setiap $i,j$ dengan $j\ne i$, akumulasi:

   ${\overrightarrow{F}}_i+=G\frac{m_i{m}_j}{r_{ij}^3}{\overrightarrow{r}}_{ij}$

3. Hitung percepatan: ${\overrightarrow{a}}_i={\overrightarrow{F}}_i/m_i$.
4. Update kecepatan: ${\overrightarrow{v}}_i+={\overrightarrow{a}}_i\cdot \mathrm{\Delta }t$.
5. Update posisi: ${\overrightarrow{r}}_i+={\overrightarrow{v}}_i\cdot \mathrm{\Delta }t$.

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

Catatan Implementasi

Kompleksitas dan Paralelisasi

Kompleksitas perhitungan langsung (direct sum) adalah $O(N^2)$ per time step. Untuk $N=10^3$, terdapat $10^6$ perhitungan pasangan per step; untuk $N=10^6$, terdapat $10^{12}$ perhitungan per step. Kebutuhan komputasi yang masif ini memotivasi penggunaan komputasi paralel menggunakan teknologi seperti CUDA, OpenCL, dan MPI.

Algoritma alternatif dengan kompleksitas lebih rendah meliputi:

• Algoritma Barnes-Hut — kompleksitas $O(N\log N)$, menggunakan struktur octree untuk meng-aproksimasi gaya dari kelompok body yang jauh.
• Fast Multipole Method — kompleksitas $O(N)$, menggunakan ekspansi multipolar.

Softening Parameter

Jika dua body mendekati jarak nol, faktor $1/r^3$ menyebabkan gaya menuju tak hingga (singularitas), yang menyebabkan simulasi tidak stabil. Untuk mengatasinya, digunakan softening parameter $ϵ$:

$$r_{ij}^{\text{soft}}=\sqrt{r_{ij}^2+{ϵ}^2}$$

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 $E=E_{\text{kinetik}}+E_{\text{potensial}}$ harus kekal:

$$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}$$

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 yang bersifat simplektik memiliki drift yang bounded (berosilasi tetapi tidak tumbuh), sehingga lebih disukai untuk simulasi jangka panjang.

Optimisasi Hukum III Newton

Karena ${\overrightarrow{F}}_{i\leftarrow j}=-{\overrightarrow{F}}_{j\leftarrow i}$, secara teoritis kita dapat menghemat 50% perhitungan dengan iterasi $j>i$ saja dan mengakumulasi gaya ke kedua body sekaligus. Namun optimisasi ini menyulitkan paralelisasi pada GPU karena menyebabkan write conflict antar thread. Pada implementasi paralel, perhitungan $N^2$ 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). "Fast N-Body Simulation with CUDA". GPU Gems 3, NVIDIA Corporation.