Apa cara yang benar untuk mengintegrasikan dalam simulasi astronomi?

15

Saya membuat simulator astronomi sederhana yang harus menggunakan fisika Newton untuk mensimulasikan pergerakan planet dalam suatu sistem (atau benda apa pun, dalam hal ini). Semua benda adalah lingkaran dalam bidang Euclidean, yang memiliki sifat seperti posisi, kecepatan, massa, jari-jari dan gaya yang dihasilkan.

Saya ingin memperbarui alam semesta dalam langkah-langkah waktu yang kecil, biasanya beberapa milidetik, tetapi saya tidak yakin bagaimana cara menghitung perubahan posisi dengan benar.

Gaya adalah sederhana: fr = sum(G * body.m * bodyi.m / dist(body, bodyi)^2).

Tetapi bagaimana saya melanjutkan dari sana?

Saya bisa melakukan ini:

a = Fr/body.m
v += a*dt
position += v*dt

Tapi tentu saja itu salah. Mungkin jika saya menambahkan 0,5 sebagai faktor dalam perhitungan posisi?

jcora
sumber
Terlalu lucu untuk tidak berkomentar: Ini memang masalah astronomi yang umum untuk mensimulasikan pergerakan "tanaman" ;-)
Wolfgang Bangerth

Jawaban:

17

Anda pada dasarnya mendapatkan jawabannya - tidak perlu faktor 0,5.

Pada dasarnya Anda memiliki sistem dua dimensi ODE orde pertama: mana semuanya adalah fungsi waktu kecuali mungkinm, dan titik menunjukkan turunan waktu. Jika Anda melakukan diferensiasi orde pertama sederhana, maju-Euler-esque ini, Anda menemukan x n + 1 -xn

x˙=vv˙=Fm,
m atau xn+1
xn+1-xnΔt=vnvn+1-vnΔt=Fnm,
Di sini saya mengindeks timestep dengann.
xn+1=xn+Δtvnvn+1=vn+ΔtFnm.
n

Namun, forward-Euler secara inheren tidak stabil. Untungnya, ada metode simplektik tepat di sudut. (Itu artikel terkait lebih dari rintisan, tapi mungkin mengandung beberapa link yang berguna.) Kuncinya adalah untuk posisi muka dari ke t n + 1 menggunakan kecepatan di t n + 1 / 2 . Artinya, misalkan Anda diberi x 0 dan v 1 / 2 untuk setiap partikel. Kemudian Anda bisa menggunakan x n + 1tntn+1tn+1/2x0v1/2

xn+1=xn+Δtvn+1/2vn+1/2=vn-1/2+ΔtFnm

v1/2v0

ω=GM.r3,
M.r

sumber
Hei, bisakah Anda menjelaskan mengapa saya tidak membutuhkan 0.5faktornya? Tampaknya melakukan hal yang sama seperti mengambil kecepatan n-1/2dtdetik yang lalu, itulah yang Anda sarankan.
jcora
(n-1)vnvn+1vn0