Bolygómozgás
Számítógépes szimulációk 1n4i11/1 Csabai István ELTE Komplex Rendszerek Fizikája Tanszék 5.102
Email: csabaiΘcomplex.elte.hu
2009 tavasz
A bolygómozgás Kepler törvényei
1
2
3
A bolygók pályája ellipszis, és annak egyik gyújtópontjában van a Nap. A bolygók vezérsugara (a bolygót a Nappal összeköt® szakasz) azonos id®k alatt azonos területet súrol. A bolygók Naptól való átlagos távolságainak, azaz a pálya fél nagytengelyeinek (a) köbei úgy aránylanak egymáshoz, mint a 3 keringési idejük (T ) négyzetei. Tehát a Ta 2 hányados minden bolygó esetén ugyanakkora (ha azok ugyanabban a naprendszerben keringenek).
Kepler az el®z® törvényeket nem elméleti úton vezette le, hanem Tycho Brahe csillagászati meg gyelései alapján találta meg. A törvényeket kés®bb Isaac Newton vezette le a gravitációs elmélete alapján. A gravitációs elmélet alapján több általánosítás is tehet®: • a törvények nem csak egy bolygó-csillag párosra, hanem bolygó
körül kering® holdakra és m¶holdakra, illetve bármely nagy tömeg¶ égitest körül kering® más égitestekre is igazak
• a természetben nem csak kötött ellipszis alakú pályák lehetésgesek,
hanem parabola és hiperbola is lehetséges
Newton gravitációs törvénye
Kepler törvényei Newton gravitációs törvényéb®l vezethet®ek le: F12 = −
Gm1 m2 r12 , 3 r12
F21 = −
Gm1 m2 r21 3 r21
A középpontba mutató er®k következtében a perdület megmarad, így a pályák síkban helyezkednek el.
A Kepler probléma egzakt megoldása • • •
a Nap tömege legalább 1000-szerese bármely bolygó tömegénél feltehetjük, hogy a Nap mozdulatlan, illetve elhanyagolhatjuk a bolygók tömegét bolygó-hold, csillag-csillag rendszerek esetén a tömegek nem hanyagólhatóak el, ekkor a tömegközéppont lesz mozdulatlan (redukált tömeg!) µ=
•
m1 m2 , r = r1 − r2 , m1 + m2
m1 r1 + m2 r2 = 0
Kepler pálya egyenlete ( a pálya excentricitása) r(θ) =
a(1 − 2 ) , 1 − cos θ
b=a
p 1 − 2
A Kepler probléma egzakt megoldása •
a pálya sebesség s v=
•
G(m1 + m2 )
2 1 − r a
a teljes energia = kinetikus + potenciális energia 1 Gm1 m2 Gm1 m2 E = µv 2 − =− 2 r 2a
•
keringési id® - Kepler 3-ik törvényéb®l T2 =
•
4π 2 a3 G(m1 + m2 )
a perdület
L = r × (µv)
A Kepler probléma numerikus számítása Két bolygó kering a közös tömegközéppontjuk körül. Határozzuk meg a pálya paramétereit!
Csillagászati egységek •
a Föld fél nagytengelye:
1AU = 1.496 · 1011 m •
Kepler 3-ik törvénye G(MSun + mEarth ) =
•
4π 2 a3 t2
a távolságot mérjük csillagászati egységben (AU ) és az id®t években (yr) G(MSun + mEarth ) = 4π 2 AU 3 /yr2
A Kepler probléma numerikus számítása Kezdeti értékek • • • •
válasszuk a nagytengelyt x tengelynek rögzítsük x(t = 0)-t a közelpontban: x(0) = a(1 + ) a közelpontban a sebesség y irányú, adjuk meg vy (t = 0)-t a sebesség a pálya egyenletb®l s v=
•
G(m1 + m2 )
2 1 − r a
az el®z® összefüggések alapján t = 0-ban meghatározhatjuk a megadott paraméterek mellett az excentricitást és a fél nagytengely méretét = x(0)vy (0)2 /(G(M + m)) − 1,
a = x(0)/(1 + )
C++ program #include
#include #include using namespace std; #include "vector.hpp" #include "odeint.hpp" using namespace cpl; const double pi = 4 * atan(1.0); const double GmPlusM = 4 * pi * pi; bool switch_t_with_y = false; // to interpolate to y = 0 // Derivative vector for Newton's law of gravitation Vector f(const Vector& x) { double t = x[0], r_x = x[1], r_y = x[2], v_x = x[3], v_y = x[4]; double rSquared = r_x*r_x + r_y*r_y; double rCubed = rSquared * sqrt(rSquared); Vector f(5); f[0] = 1; f[1] = v_x; f[2] = v_y; f[3] = - GmPlusM * r_x / rCubed; f[4] = - GmPlusM * r_y / rCubed; if (switch_t_with_y) { // use y as independent variable for (int i = 0; i < 5; i++) f[i] /= v_y; } return f; }
A pálya vizsgálatához ki akarjuk íratni az x-tengely metszeteket. Mivel a véges lépések miatt nem lenne mindig kiértékelés itt, amikor keresztezzük az x-tengelyt, visszaléptatünk oda. A következ® függvény vissza-interpolálja a pályát egy x-tengely metszése után az y = 0-hoz. • ehhez áttérünk t-r®l y változóra: dx dx/dt = dy dy/dt •
a Runge-Kutta lépésköz legyen -y
// Change independent variable from t to y and step back to y = 0 void interpolate_crossing(Vector x, int& crossing) { ++crossing; switch_t_with_y = true; RK4Step(x, -x[2], f); cout << " crossing " << crossing << "\t t = " << x[0] << "\t x = " << x[1] << endl; switch_t_with_y = false; }
//The code contains both simple and adaptive Runge-Kutta method int main() { cout << " Kepler orbit comparing fixed and adaptive Runge-Kutta\n" << " -----------------------------------------------------\n" << " Enter aphelion distance in AU, and eccentricity: "; double r_ap, eccentricity, a, T, v0; cin >> r_ap >> eccentricity; a = r_ap / (1 + eccentricity); T = pow(a, 1.5); v0 = sqrt(GmPlusM * (2 / r_ap - 1 / a)); cout << " Enter number of periods, step size, and adaptive accuracy: "; double periods, dt, accuracy; cin >> periods >> dt >> accuracy; Vector x0(5); x0[0] = 0; x0[1] = r_ap; x0[2] = 0; x0[3] = 0; x0[4] = v0; ofstream dataFile("fixed.data"); Vector x = x0; int steps = 0, crossing = 0; cout << "\n Integrating with fixed step size" << endl; do { for (int i = 0; i < 5; i++) dataFile << x[i] << '\t'; dataFile << '\n'; double y = x[2]; RK4Step(x, dt, f); ++steps; if (y * x[2] < 0) interpolate_crossing(x, crossing); } while (x[0] < periods * T); cout << " number of fixed size steps = " << steps << endl; cout << " data in file fixed.data" << endl; dataFile.close();
}
dataFile.open("adaptive.data"); x = x0; steps = crossing = 0; double dt_max = 0, dt_min = 100; cout << "\n Integrating with adaptive step size" << endl; do { for (int i = 0; i < 5; i++) dataFile << x[i] << '\t'; dataFile << '\n'; double t_save = x[0]; double y = x[2]; adaptiveRK4Step(x, dt, accuracy, f); double step_size = x[0] - t_save; ++steps; if (step_size < dt_min) dt_min = step_size; if (step_size > dt_max) dt_max = step_size; if (y * x[2] < 0) interpolate_crossing(x, crossing); } while (x[0] < periods * T); cout << " number of adaptive steps = " << steps << endl; cout << " step size: min = " << dt_min << " max = " << dt_max << endl; cout << " data in file adaptive.data" << endl; dataFile.close();
Három test probléma Egzakt megoldások • a probléma egyszer¶nek t¶nik: három tömeg (m1 , m2 , m3 ), három pozició (r~1 , r~2 , r~3 ) és három egyenlet : ~r1 − ~r2 ~r1 − ~r3 d2~r1 = −Gm2 − Gm3 dt2 |~r1 − ~r2 | |~r1 − ~r3 |
• • •
(ez az 1. testre vonatkozik, hasonló egyenletek írhatóak fel a 2. és 3. testre) a megoldások azonban er®sen függenek a kezdeti feltételekt®l! analitikusan csak nagyon kevés eset kezelhet® kezelhet® egyszer¶bb esetek: • síkbeli 3 test probléma - a 3 test egy síkban mozog • korlátozott 3 test probléma - m3 = 0, ekkor m1 , m2 Kepler
pályán mozog, míg a 3-ik körülöttük kering
Síkbeli 3 test probléma
általában a 3 test problémák nem síkbeliek, de pl. a Naprendszer kezelhet® síkbeliként • gyorsulások ~a1 = − •
Gm3 Gm2 ~r12 − 2 ~r13 2 r12 r13
relatív poziciók ~s1 = ~r3 − ~r2 , ~s2 = ~r1 − ~r3 , ~s3 = ~r2 − ~r1 ~s1 + ~s2 + ~s3 = 0
•
mozgás egyenletek d2~si mG ~ = − 2 ~si + mi G 2 dt s ~ =G m = m1 + m2 + m3 , G
3 X ~si 2 s i=1 i
Linkek
Nap-Föld-Jupiter applet: t/T : keringési id®, r: távolság, v : sebesség állíthatóak: mi : a Jupiter és a Föld tömege, ei excentricitásuk, a: a Föld fél nagytengelye, theta: a Föld kezdeti pálya szöge, a keringési iránya, valamint a egy csúszkával az animáció sebessége.
Feladatok Használjuk a kepler.cpp programot! A nagytengely iránya és nagysága elvileg állandó. Teszteljük, hogy a numerikus hibák miatt ez változik-e, és ha igen, akkor mennyire. Vizsgáljuk meg, hogy a lépéshossz és az alkalmazott integrál módszer mennyire befolyásolja ezt. (Keressük meg a trajektória perihéliumát és mérjük le a nagytengely szögét és hosszát a kiindulási helyzethez képest.) Vizsgáljuk meg, hogy az adaptív lépéshossz hogyan változik a pálya mentén! Ábrázoljuk és magyarázzuk meg a kapott eredményt! Szorgalmi feladat: Hasonlítsuk össze, hogy ugyanakkora precizitás megkövetelésekor (pl. 10 teljes pálya után ugyanakkora hiba) mennyi futási id® szükséges a normál és az adaptív Runge-Kutta szimulációhoz. A Merkúr perihélium precessziója (lapozz!) 1
2
3
Feladatok 3
A Merkúr perihélium precessziója. A Merkúr pályája er®sen elnyúlt (excentricitása: = 0.2056), így a pályán hosszú id® alatt meg gyelhet® annak elfordulása. Erre az elfordulásra az általános relativitás elmélet adott magyarázatot. Ennek alapján a következ® korrekció szükséges az er®k megadásakor: F =
4
GmSum mM ercury α 1 + r2 r2
Írjuk át a kódot ennek a korrekciónak a gyelembevételére. Vizsgáljuk meg, hogy α ≈ 1.1 · 10−8 AU 2 érték mellett mekkora perihélium mozgást tapasztalunk. A kapott értéket vessük össze a a Merkúr pályájának évszázadonkénti 43 ívmásodperces mért elfordulásával. (tipp: kis α esetén a precesszó rátája ∝ α, azaz meghatározhatjuk a meredekséget és extrapolálhatunk α ≈ 1.1 · 10−8 AU 2 -ra) Alakítsd át a kepler.cpp forráskódot, hogy 3 test problémát kezeljen. Próbálj ki különféle paramétereket, ábrázold a pályákat, vizsgáld a stabilitást!
Feladatok
5
Szorgalmi feladat: Largrange pontok (lásd Wikipédia) stabilitása. Szimuláljuk a Nap-Jupiter rendszert, és vizsgáljuk meg az L1, L2, L3, L4, L5 Lagrange pontokba helyezett kis (a bolygó tömegénél sokkal kisebb) tömegek mozgását, igazoljuk, hogy (csak) a L4-L5 pontok stabilak. A Jupiter esetében az utóbbi pontokban valóban összegy¶lnek aszteroibák, az ún Trójai kisbolygók (linkek: trójaiak az SDSS-ben, wikipedia).
