Springen naar inhoud

Probleem bij baansimulatie


  • Log in om te kunnen reageren

#1

Brinx

    Brinx


  • >1k berichten
  • 1433 berichten
  • Lorentziaan

Geplaatst op 21 mei 2007 - 15:15

Ik ben op het moment bezig met het schrijven van een relatief eenvoudig programmaatje dat de gebruiker in staat stelt via een grafische interface een simulatie te runnen van een n-lichamenprobleem (de beweging van n massa's als gevolg van hun wederzijdse gravitatiekrachten). De gebruiker kan met die interface objecten plaatsen met posities, snelheden, massa's en radii naar keuze, en hun bewegingen vervolgens simuleren met een druk op de knop.

Voor de daadwerkelijke simulatie zelf ben ik bezig geweest om een aantal verschillende integratoren te implementeren, om te kijken welke van de mogelijke opties de meest geschikte is. Tot zover heb ik de volgende integratoren gebruikt:

- Euler
- Leapfrog (2e orde)*
- Runge-Kutta (4e orde)
- Yoshida (4e orde)*

De twee integratoren die ik aangaf met een sterretje zijn symplectisch: reversibel en energiebehoudend, voor situaties als deze. Om de 'symplecticiteit' te testen houd ik voor iedere integratiestap de totale energie van het systeem bij (de som van de kinetische energie van alle massa's plus de som van alle relatieve potentialen), en kijk ik of die inderdaad praktisch constant blijft (of eigenlijk slechts kleine periodieke schommelingen vertoont, en geen drift).

Het rare is nu dat ik, na een serie testjes, merk dat bijvoorbeeld de leapfrog-integrator slechts eerste-orde gedrag vertoont (de fout in de totale energie is recht evenredig met de tijdstapgrootte), waar hij uiteraard tweede-orde gedrag (dus energiefout is evenredig met het kwadraat van de stapgrootte) zou moeten vertonen. Dit gedrag komt alleen maar voor wanneer ik een simulatie laat lopen waarin twee of meer lichamen vrij in een vlak bewegen. Wanneer ik de simulatie run zodanig dat er een van de twee massa's vastgehouden wordt (dus wanneer er slechts een enkele massa door de integrator 'onder handen genomen' wordt) vertoont de leapfrog-integrator gewoon het gewenste (en verwachte) gedrag.

Ik heb het vermoeden gekregen dat het te maken heeft met de manier waarop meerdere posities en versnellingen in een tijdstap berekend worden. Als toelichting zal ik in het kort beschrijven wat de leapfrog-integrator precies doet.

Stel de positievector van een massa 'i' in de simulatie op tijdstip 't' voor als LaTeX , de snelheidsvector als LaTeX , en de versnelling van dat object als LaTeX . De leapfrog-integrator doet om de volgende positie en snelheid te berekenen de volgende dingen:

LaTeX

dan

LaTeX

(de functie 'f' stelt de berekening van de versnelling voor als functie van de positie) en vervolgens

LaTeX

Deze procedure werkt dus prima wanneer er slechts sprake is van een bewegend object en een vaste centrale aantrekker. Wanneer er echter meerdere bewegende objecten zijn ontstaat er een probleem. De volgende code voor twee objecten levert bijvoorbeeld niet het gewenste resultaat op:

LaTeX
LaTeX

dan

LaTeX
LaTeX

(beide versnellingen worden hier berekend gebruik makend van de nieuwe posities) en vervolgens

LaTeX
LaTeX

Ook wanneer ik eerst alleen de positie en snelheid van massa 'i' integreer en daarna pas de positie en snelheid van massa 'j' integreer (gebruik makend van de situatie op tijdstip 't'), komt er niet het gewenste (tweede-orde) gedrag uit.

Wie heeft er een idee van de correcte implementatie van dit algoritme wanneer het om meerdere objecten gaat?

Dit forum kan gratis blijven vanwege banners als deze. Door te registeren zal de onderstaande banner overigens verdwijnen.

#2

Dr. Who?

    Dr. Who?


  • >250 berichten
  • 305 berichten
  • Ervaren gebruiker

Geplaatst op 21 mei 2007 - 15:28

Hum, even een eerste ideetje: zou het kunnen dat de richting van de versnelling niet helemaal correct berekend wordt?

#3

Brinx

    Brinx


  • >1k berichten
  • 1433 berichten
  • Lorentziaan

Geplaatst op 21 mei 2007 - 15:44

Snelle reactie! :(

De richting van de versnellingen klopt wel: objecten worden precies naar elkaar toe getrokken, volgens de welbekende relaties. De afwijking die ik krijg in totale energie manifesteert zichzelf via langzaam naar binnen spiraliserende banen voor de configuratie die ik nu gebruik bij wijze van test (klein object in lage, vrijwel circulaire aardbaan). Ook de vereiste nauwkeurigheid zit volgens mij wel goed: ik heb nergens in de berekeningen te maken met extreem kleine of grote getallen die tegen de limiet van de variabelen aanzitten.

#4

Dr. Who?

    Dr. Who?


  • >250 berichten
  • 305 berichten
  • Ervaren gebruiker

Geplaatst op 21 mei 2007 - 16:02

Hehe - hum, ja, ik, eh, zit even wat te relaxen na een telefoongesprek van een uur met Thales Alenia... :( & :(

Nee, maar goed - als ze naar binnen bewegen, dan zou je haast verwachten dat een component van de versnelling in de negatieve along-track richting wijst. Het eerste wat dus in me opkwam, was dat op de ťťn of andere manier een toekomstige positievector x wordt gebruikt bij het berekenen van de versnellingsvector. In dat geval zou de berekende versnellingsvector dus iets naar achteren wijzen en opzichte van de correcte versnellingsvector - aan de andere kant zou dat probleem ook moeten optreden als er maar ťťn lichaam beweegt...

Enniewee, misschien is het een idee om een testcase te maken met twee even grote massa's die in cirkelbanen om een gemeenschappelijk zwaartepunt draaien en voor twee opeenvolgende tijdstappen de snelheids-, en versnellingsvectoren te plotten, plus de verwachte baan?

Veranderd door Dr. Who?, 21 mei 2007 - 16:02


#5

Brinx

    Brinx


  • >1k berichten
  • 1433 berichten
  • Lorentziaan

Geplaatst op 21 mei 2007 - 19:58

Goeie suggestie, daar was ik inderdaad al mee begonnen. :( Updates over voortgang volgen later!

#6

Dr. Who?

    Dr. Who?


  • >250 berichten
  • 305 berichten
  • Ervaren gebruiker

Geplaatst op 21 mei 2007 - 22:48

Ik heb je algoritme geimplementeerd als MATLAB-script, en het lijkt redelijk te werken. :( Het zijn wel fascinerende dingen, die symplectische integratoren - heb je misschien ook het algoritme voor de Yoshida integrator? Heeft die variabele stapgrootte? En, last but not least, kan je er dense output mee genereren?

function test_symplectic
% Functie voor het testen van een symplectische integrator


% Initialisatie
clc
clear
close all

% Simulatie
T = 1; % Baanperiode
steps = 400; % aantal stappen per orbit
N_orb = 10; % aantal orbits
step_tot = N_orb*steps; % totaal aantal stappen
dt = T/steps; % stapgrootte
dt2 = .5*dt^2;
time = [0:dt:N_orb-dt]';

% Aanmaken van de variabelen
r1 = zeros(step_tot,3); r2 = r1; v1 = r1; v2 = r1; a1 = r1; a2 = r1;
H_bar = r1; H = zeros(step_tot,1); E = H;

% Omgevingsconstanten

mu = 4*pi^2; % Zwaartekrachtsparameter
m1 = .7; % Massa 1
m2 = 1-m1; % Massa 2

% Simulatievariabelen
d = 1; % afstand van de massa's tot het gemeenschappelijke zwaartepunt
v = sqrt(mu/d); % snelheid van de massa's t.o.v. het gemeenschappelijke zwaartepunt (2*pi)

% Begincondities
r1(1,: ) = [0 -d 0]*m2; r2(1,: ) = [0 d 0]*m1; % Positie op t = 0;
v1(1,: ) = [v 0 0]*m2; v2(1,: ) = [-v 0 0]*m1; % Snelheid op t = 0;
a1(1,: ) = (get_acc(r1(1,: ),r2(1,: ),mu,m1,m2))/m1; a2(1,: ) = -a1(1,: )*m2/m1; % Versnelling op t = 0;

H_bar(1,: ) = m1*cross(r1(1,: ),v1(1,: ))+m2*cross(r2(1,: ),v2(1,: )); % Impulsmomentvector op t = 0;
H(1) = norm(H_bar(1,: )); % Impulsmoment op t = 0;

r = norm(r1(1,: )-r2(1,: )); % Relatieve positie
E(1) = .5*m1*v1(1,: )*v1(1,: )' + .5*m2*v2(1,: )*v2(1,: )' - m1*m2*mu/r; % Energie op t = 0;
E_norm = ones(step_tot,1)*E(1); % Verondersteld constante energie

% ----------
% Integratie
% ----------


for n = 1:step_tot-1

% Positie update
r1(n+1,: ) = r1(n,: ) + v1(n,: )*dt + a1(n,: )*dt2;
r2(n+1,: ) = r2(n,: ) + v2(n,: )*dt + a2(n,: )*dt2;

% Versnelling update
a1(n+1,: ) = (get_acc(r1(n+1,: ),r2(n+1,: ),mu,m1,m2))/m1;
a2(n+1,: ) = -a1(n+1,: )*m1/m2;

% Snelheid update
v1(n+1,: ) = v1(n,: ) + .5*(a1(n,: )+a1(n+1,: ))*dt;
v2(n+1,: ) = v2(n,: ) + .5*(a2(n,: )+a2(n+1,: ))*dt;

% Impulsmomentvector en impulsmoment update
H_bar(n+1,: ) = m1*cross(r1(n+1,: ),v1(n+1,: ))+m2*cross(r2(n+1,: ),v2(n+1,: ));
H(n+1) = norm(H_bar(n,: ));

r = norm(r1(n+1,: )-r2(n+1,: )); % Relatieve positievector

% Update energie
E(n+1) = .5*m1*v1(n+1,: )*v1(n+1,: )' + .5*m2*v2(n+1,: )*v2(n+1,: )' - m1*m2*mu/r;

end

% -------
% Figuren
% -------


figure
axis equal
hold on
plot(r1(:,1),r1(:,2),'b')
plot(r2(:,1),r2(:,2),'r')
title('Posities van massa 1 & 2 in inertiele ruimte')
xlabel('X')
ylabel('Y')
legend('Positie massa 1','Positie massa 2')

figure
plot(time,H)
grid on
title('Impulsmoment als functie van de tijd')
xlabel('tijd [omwentelingen]')
ylabel('Impulsmoment')

figure
hold on
grid on
title('Energie als functie van de tijd')
plot(time,E_norm,'r')
plot(time,E)
legend('Correct','Na integratie')
xlabel('tijd [omwentelingen]')
ylabel('Energie')

% --------------------------------
% Subfunctie voor de zwaartekracht
% --------------------------------


function F = get_acc(r1,r2,mu,m1,m2);
% Berekening van de zwaartekracht
% -------------------------------


% Relatieve positie
r12_bar = r2-r1;

% Afstand van massa 1 tot massa 2
r_12 = norm(r12_bar);

% Zwaartekrachtsvector
F = mu*m1*m2*r12_bar/r_12^3;

Veranderd door Dr. Who?, 21 mei 2007 - 22:48


#7

Brinx

    Brinx


  • >1k berichten
  • 1433 berichten
  • Lorentziaan

Geplaatst op 22 mei 2007 - 10:08

Bedankt voor de test! Ik heb nu een redelijk vermoeden van wat er mankeert aan mijn implementatie, maar ik probeer het eerst even uit voordat ik van de toren blaas. :(

Info over symplectische integratoren, specifiek Yoshida:

http://www.artcompsc...lem_2/ch07.html
(handige tekst, wordt verwezen naar het oorspronkelijke paper uit 1990, en de tekst bevat nog veel meer info over andere typen integratoren)

http://www.aoc.nrao....g/prog/symp.cpp
(c++-implementatie van Yoshida's algoritmen tot en met orde 8)

#8

Brinx

    Brinx


  • >1k berichten
  • 1433 berichten
  • Lorentziaan

Geplaatst op 22 mei 2007 - 14:46

Juist, nu heb ik de boel fatsoenlijk aan het lopen. Een stukje output:

Energiefouten voor Euler-integratie:

Initial energy = -3.04
time step = 0.1
dE = 0.25931799078552675, t = 1.0 (10 steps)

Initial energy = -3.04
time step = 0.01
dE = 0.027934765632875447, t = 1.0 (100 steps)

Initial energy = -3.04
time step = 0.001
dE = 0.0028174131077802755, t = 1.0 (1000 steps)

Initial energy = -3.04
time step = 1.0E-4
dE = 2.8198534160850386E-4, t = 1.0 (10000 steps)

Initial energy = -3.04
time step = 1.0E-5
dE = 2.820143415505072E-5, t = 1.0 (100000 steps)

Het is hier goed te zien hoe de energiefout dE telkens een factor 10 kleiner wordt (over eenzelfde integratie-interval) wanneer de stapgrootte een factor 10 verkleind wordt.

Energiefouten voor Leapfrog-integratie:

Initial energy = -3.04
time step = 0.1
dE = 0.0015886737403736362, t = 1.0 (10 steps)

Initial energy = -3.04
time step = 0.01
dE = 1.6105538297672695E-5, t = 1.0 (100 steps)

Initial energy = -3.04
time step = 0.0010
dE = 1.6107747846660914E-7, t = 1.0 (1000 steps)

Initial energy = -3.04
time step = 1.0E-4
dE = 1.610767963455828E-9, t = 1.0 (10000 steps)

Initial energy = -3.04
time step = 1.0E-5
dE = 1.6142642778049776E-11, t = 1.0 (100000 steps)

Hier is goed te zien hoe de energiefout bij het 10x kleiner maken van de tijdstap een factor 100 kleiner wordt: tweede-orde gedrag, zoals je zou verwachten.

Voor de volledigheid de beginsituatie van de simulatie:

positievectoren:
r1(t = 0) = ( 1.0, 0.0)
r2(t = 0) = (-1.0, 0.0)

snelheidsvectoren:
v1(t = 0) = ( 0.0, -0.8)
v2(t = 0) = ( 0.0, 0.4)

massa's:
m1 = 2.0
m2 = 4.0

De gravitatieconstante G is voor het gemak op 1 gezet voor deze simulatie.

#9

Dr. Who?

    Dr. Who?


  • >250 berichten
  • 305 berichten
  • Ervaren gebruiker

Geplaatst op 22 mei 2007 - 15:38

En nu de plaatjes d'rbij, hahaha! :(

Ik heb m'n simulatie ook maar even wat uitgebreid met Yoshida-8 code en een fout in de beginconditie voor de versnelling verbeterd.

a2(1,:) = -a1(1,:)*m2/m1;

moet worden:

a2(1,:) = -a1(1,:)*m1/m2;

Het lijkt me op zich niet verkeerd om zoiets te gebruiken als voorbeeld in de minicursus baanmechanica in spe...

Kan je m-files eigenlijk ook als attachment bij je bericht invoegen?

#10

Brinx

    Brinx


  • >1k berichten
  • 1433 berichten
  • Lorentziaan

Geplaatst op 22 mei 2007 - 15:51

Over die attachment-mogelijkheden weet ik niet zoveel eigenlijk...

Maar wat misschien wel leuk zou zijn is om het java-programmaatje dat ik aan het maken ben eventueel als applet in te voegen op een uitlegpagina over baanmechanica. Ik denk niet dat zoiets in een minicursus kan (java is denk ik niet toegestaan ivm veiligheid), maar op een externe pagina zou dat geen probleem moeten zijn.

#11

Dr. Who?

    Dr. Who?


  • >250 berichten
  • 305 berichten
  • Ervaren gebruiker

Geplaatst op 22 mei 2007 - 19:24

Geplaatste afbeelding

#12

Jan van de Velde

    Jan van de Velde


  • >5k berichten
  • 44863 berichten
  • Moderator

Geplaatst op 22 mei 2007 - 22:20

joi, interactieve minicursusjes. :grin:
off-topic: waarom zou het trouwens onveilig zijn om applets te hosten op wsf?
(dus gooi voor deze vraag dit topic niet overhoop)
ALS WIJ JE GEHOLPEN HEBBEN....
help ons dan eiwitten vouwen, en help mee ziekten als kanker en zo te bestrijden in de vrije tijd van je chip...
http://www.wetenscha...showtopic=59270

#13

Brinx

    Brinx


  • >1k berichten
  • 1433 berichten
  • Lorentziaan

Geplaatst op 23 mei 2007 - 11:18

Pfff, EINDELIJK gevonden waar de fout lag. Voor degenen die het willen weten: wanneer je een Vector gelijkstelt aan een andere Vector in Java (dus wanneer je "Vector1 = Vector2;" uit laat voeren) wordt er slechts een referentie gemaakt naar de originele Vector, en wordt er geen daadwerkelijke kopie gemaakt. Duurde wel even voordat ik dat doorhad...

Maar nu werkt het dus perfect! Ik poets nog wat meer aan de interface, en dan zal ik kijken hoe ik het hier beschikbaar kan maken.

Ik zal overigens eens vragen in het relevante subforum hoe het zit met het invoegen van java applets in posts. Ik geef 't niet veel kans - aangezien de browser vooraf niet kan 'ruiken' wat de applet precies gaat doen kan dat nare gevolgen hebben.





0 gebruiker(s) lezen dit onderwerp

0 leden, 0 bezoekers, 0 anonieme gebruikers

Ook adverteren op onze website? Lees hier meer!

Gesponsorde vacatures

Vacatures