Célom egy olyan tudományos kísérlet szimulálása, ami 3 bolygó gravitációs problémáját modellezi, az úgynevezett N-test problémát szimulálja numerikus integrálással. Ez a probléma a csillagászat, asztrofizika és számítógépes fizika egyik alapvető kutatási területe. A program a Runge–Kutta 4. rendű módszerrel (RK4) számolja ki az égitestek mozgását.A program alkalmazza a Newton-féle gravitációs törvényt, numerikus differenciálegyenlet-megoldással, egy Runge–Kutta 4. rendű integrátorral, az N-test dinamikát alkalmazva, egy asztrofizikai szimulációban. A program továbbfejleszthető relativisztikus korrekciókkal, szimplektikus integrátorokkal. A Nap-Föld-Jupiter rendszer egy jó példája a fizikában és csillagászatban vizsgált N-test problémának (általában háromtest-problémának). Mivel az égitestek kölcsönösen vonzzák egymást, pályáik matematikai leírása rendkívül összetett, és gyakran kaotikus rendszert eredményez. A szimulációt remound nélkül valósítottam meg, ami igazi kihívás a programozók számára.
--------------
import numpy as np
import matplotlib.pyplot as plt
G = 6.67430e-11 # gravitációs állandó
class Body:
def __init__(self, mass, position, velocity, name):
self.mass = mass
self.position = np.array(position, dtype=float)
self.velocity = np.array(velocity, dtype=float)
self.name = name
def accelerations(bodies):
n = len(bodies)
acc = [np.zeros(2) for _ in range(n)]
for i in range(n):
for j in range(n):
if i == j:
continue
r = bodies[j].position - bodies[i].position
dist = np.linalg.norm(r)
if dist > 0:
acc[i] += G * bodies[j].mass * r / dist**3
return acc
def rk4_step(bodies, dt):
original_positions = [b.position.copy() for b in bodies]
original_velocities = [b.velocity.copy() for b in bodies]
# k1
a1 = accelerations(bodies)
k1v = [a * dt for a in a1]
k1x = [b.velocity * dt for b in bodies]
# k2
for i, b in enumerate(bodies):
b.position = original_positions[i] + k1x[i] / 2
b.velocity = original_velocities[i] + k1v[i] / 2
a2 = accelerations(bodies)
k2v = [a * dt for a in a2]
k2x = [b.velocity * dt for b in bodies]
# k3
for i, b in enumerate(bodies):
b.position = original_positions[i] + k2x[i] / 2
b.velocity = original_velocities[i] + k2v[i] / 2
a3 = accelerations(bodies)
k3v = [a * dt for a in a3]
k3x = [b.velocity * dt for b in bodies]
# k4
for i, b in enumerate(bodies):
b.position = original_positions[i] + k3x[i]
b.velocity = original_velocities[i] + k3v[i]
a4 = accelerations(bodies)
k4v = [a * dt for a in a4]
k4x = [b.velocity * dt for b in bodies]
# végső frissítés
for i, b in enumerate(bodies):
b.position = (
original_positions[i]
+ (k1x[i] + 2*k2x[i] + 2*k3x[i] + k4x[i]) / 6
)
b.velocity = (
original_velocities[i]
+ (k1v[i] + 2*k2v[i] + 2*k3v[i] + k4v[i]) / 6
)
# Nap–Föld–Jupiter rendszer
sun = Body(
mass=1.989e30,
position=[0, 0],
velocity=[0, 0],
name="Nap"
)
earth = Body(
mass=5.972e24,
position=[1.496e11, 0],
velocity=[0, 29780],
name="Föld"
)
jupiter = Body(
mass=1.898e27,
position=[7.785e11, 0],
velocity=[0, 13070],
name="Jupiter"
)
bodies = [sun, earth, jupiter]
trajectories = {b.name: [] for b in bodies}
dt = 24 * 3600 # 1 nap
steps = 3650 # 10 év
for _ in range(steps):
for b in bodies:
trajectories[b.name].append(b.position.copy())
rk4_step(bodies, dt)
# Ábrázolás
plt.figure(figsize=(10, 10))
for name, traj in trajectories.items():
traj = np.array(traj)
plt.plot(traj[:, 0], traj[:, 1], label=name)
plt.scatter([0], [0], s=200, label="Nap")
plt.xlabel("x [m]")
plt.ylabel("y [m]")
plt.title("Háromtest-probléma szimulációja (RK4)")
plt.legend()
plt.axis("equal")
plt.grid()
plt.show()
print("N-test gravitációs szimuláció indítása...")
print(f"Égitestek száma: {len(bodies)}")
print(f"Szimuláció hossza: {steps} nap")
print()
for step in range(steps):
# 100 lépésenként állapotjelzés
if step % 100 == 0:
print(f"Folyamat: {100 * step / steps:.1f}%")
for b in bodies:
trajectories[b.name].append(b.position.copy())
rk4_step(bodies, dt)
print("\nSzimuláció befejeződött.")
print("Pályaadatok kiszámítva.")
print("Grafikon megjelenítése...")
print("N-test gravitációs szimuláció indítása...")
print(f"Égitestek száma: {len(bodies)}")
print(f"Szimuláció hossza: {steps} nap")
print()
for step in range(steps):
# 100 lépésenként állapotjelzés
if step % 100 == 0:
print(f"Folyamat: {100 * step / steps:.1f}%")
for b in bodies:
trajectories[b.name].append(b.position.copy())
rk4_step(bodies, dt)
print("\nSzimuláció befejeződött.")
print("Pályaadatok kiszámítva.")
print("Grafikon megjelenítése...")
print("=" * 60)
print(" NUMERIKUS GRAVITÁCIÓS SZIMULÁTOR")
print("=" * 60)
for b in bodies:
print(f"{b.name:10} | tömeg = {b.mass:.3e} kg")
print("-" * 60)
print("Runge–Kutta 4. rendű integrátor aktiválva")
print(f"Időlépés: {dt/3600:.1f} óra")
print("-" * 60)
print("\nVégső pozíciók:")
for b in bodies:
print(
f"{b.name:10} "
f"x = {b.position[0]:.3e} m, "
f"y = {b.position[1]:.3e} m"
)
print("\nMinden számítás sikeresen lezajlott.")
----------------
N-test gravitációs szimuláció indítása...
Égitestek száma: 3
Szimuláció hossza: 3650 nap
Folyamat: 0.0%
Folyamat: 2.7%
Folyamat: 5.5%
Folyamat: 8.2%
Folyamat: 11.0%
Folyamat: 13.7%
Folyamat: 16.4%
Folyamat: 19.2%
Folyamat: 21.9%
Folyamat: 24.7%
Folyamat: 27.4%
Folyamat: 30.1%
Folyamat: 32.9%
Folyamat: 35.6%
Folyamat: 38.4%
Folyamat: 41.1%
Folyamat: 43.8%
Folyamat: 46.6%
Folyamat: 49.3%
Folyamat: 52.1%
Folyamat: 54.8%
Folyamat: 57.5%
Folyamat: 60.3%
Folyamat: 63.0%
Folyamat: 65.8%
Folyamat: 68.5%
Folyamat: 71.2%
Folyamat: 74.0%
Folyamat: 76.7%
Folyamat: 79.5%
Folyamat: 82.2%
Folyamat: 84.9%
Folyamat: 87.7%
Folyamat: 90.4%
Folyamat: 93.2%
Folyamat: 95.9%
Folyamat: 98.6%
Szimuláció befejeződött.
Pályaadatok kiszámítva.
Grafikon megjelenítése...
N-test gravitációs szimuláció indítása...
Égitestek száma: 3
Szimuláció hossza: 3650 nap
Folyamat: 0.0%
Folyamat: 2.7%
Folyamat: 5.5%
Folyamat: 8.2%
Folyamat: 11.0%
Folyamat: 13.7%
Folyamat: 16.4%
Folyamat: 19.2%
Folyamat: 21.9%
Folyamat: 24.7%
Folyamat: 27.4%
Folyamat: 30.1%
Folyamat: 32.9%
Folyamat: 35.6%
Folyamat: 38.4%
Folyamat: 41.1%
Folyamat: 43.8%
Folyamat: 46.6%
Folyamat: 49.3%
Folyamat: 52.1%
Folyamat: 54.8%
Folyamat: 57.5%
Folyamat: 60.3%
Folyamat: 63.0%
Folyamat: 65.8%
Folyamat: 68.5%
Folyamat: 71.2%
Folyamat: 74.0%
Folyamat: 76.7%
Folyamat: 79.5%
Folyamat: 82.2%
Folyamat: 84.9%
Folyamat: 87.7%
Folyamat: 90.4%
Folyamat: 93.2%
Folyamat: 95.9%
Folyamat: 98.6%
Szimuláció befejeződött.
Pályaadatok kiszámítva.
Grafikon megjelenítése...
============================================================
NUMERIKUS GRAVITÁCIÓS SZIMULÁTOR
============================================================
Nap | tömeg = 1.989e+30 kg
Föld | tömeg = 5.972e+24 kg
Jupiter | tömeg = 1.898e+27 kg
------------------------------------------------------------
Runge–Kutta 4. rendű integrátor aktiválva
Időlépés: 24.0 óra
------------------------------------------------------------
Végső pozíciók:
Nap x = 1.477e+09 m, y = 1.198e+10 m
Föld x = 1.508e+11 m, y = 2.156e+10 m
Jupiter x = -7.695e+11 m, y = -1.043e+11 m
Minden számítás sikeresen lezajlott
Pályaadatok kiszámítva.
Ready---------------
Futató környezet; https://onecompiler.com/python#draft-6krb
Az N-test gravitációs szimuláció olyan számítógépes modell, amely sokaságok (égitestek vagy részecskék) egymásra gyakorolt gravitációs vonzását számolja ki lépésről lépésre. A szimuláció minden egyes időpillanatban meghatározza a testekre ható összes erőt, majd ez alapján frissíti a pozíciójukat és sebességüket.Alapelv és működésA szimuláció alapja a fizika törvényei által meghatározott mozgásegyenlet. Két test közötti gravitációs erő (F) a Newton-féle gravitációs törvény szerint számítható:
a gravitációs állandóm₁, m₂: a két test tömeger: a testek közötti távolságA szimuláció lépései egy adott dt (időlépcső) időtartamra a következők:Erőhatások kiszámítása: A rendszerben lévő összes testre ható erőt összegezzük a többi testtől.Gyorsulás meghatározása: A F = m ⋅ a képlet alapján minden test gyorsulását kiszámítjuk.Sebesség frissítése: A gyorsulásból és a meglévő sebességből meghatározzuk az új sebességet.Pozíció frissítése: Az új sebesség alapján a testek elmozdulnak új pozíciójukba. Nehézségek (Komplexitás)Lépésszám növekedése: Ha N darab test van, minden egyes testre N-1 másik test hat. Egyetlen lépésben az elvégzendő műveletek száma nagyságrendileg N². Ez rengeteg test (pl. galaxisok több millió csillaga) esetén rendkívül lassú lehet.Optimalizáció: Nagy számú objektumnál a szimuláció gyorsítására olyan algoritmusokat használnak, mint a Barnes–Hut, amely a távoli testeket összevonva, egyetlen pontként kezeli.Gyakorlati alkalmazásokAsztrofizika: Galaxisok ütközése, bolygópályák stabilitása, csillaghalmazok fejlődése és a sötét anyag tágulási folyamatainak vizsgálata.Kozmológia: A Világegyetem nagyléptékű struktúrájának evolúciója. Kipróbálható szimulációkHa szeretné saját maga is tesztelni vagy vizuálisan megfigyelni, az interneten számos interaktív és nyílt forráskódú opció érhető el: Alapszintű oktató modell: Interaktív és vizuális számításokhoz használható a Desmos N-test szimuláció webes kalkulátor Asztrofizikai részecskeszimulátor: A Gadget egy széles körben használt kozmológiai kód, amellyel sötét anyag és gázok gravitációs dinamikáját modellezik szuperszámítógépeken. (https://pypi.org/project/rebound/)
Nincsenek megjegyzések:
Megjegyzés küldése