Probleem bij baansimulatie

Moderators: Michel Uphoff, jkien

Reageer
Gebruikersavatar
Lorentziaan
Berichten: 1.433

Probleem bij baansimulatie

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-*** (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 \(x_i(t)\), de snelheidsvector als \(v_i(t)\), en de versnelling van dat object als \(a_i(t)\). De leapfrog-integrator doet om de volgende positie en snelheid te berekenen de volgende dingen:
\(x_i(t+1) = x_i(t) + v_i(t) \cdot \Delta t + a_i(t) \cdot \frac{\Delta t^2}{2}\)
dan
\(a_i(t+1) = f(x_i(t+1))\)
(de functie 'f' stelt de berekening van de versnelling voor als functie van de positie) en vervolgens
\(v_i(t+1) = v_i(t) + \frac{a_i(t) + a_i(t+1)}{2}\cdot \Delta t\)
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:
\(x_i(t+1) = x_i(t) + v_i(t) \cdot \Delta t + a_i(t) \cdot \frac{\Delta t^2}{2}\)
\(x_j(t+1) = x_j(t) + v_j(t) \cdot \Delta t + a_j(t) \cdot \frac{\Delta t^2}{2}\)
dan
\(a_i(t+1) = f(x_i(t+1),x_j(t+1))\)
\(a_j(t+1) = f(x_i(t+1),x_j(t+1))\)
(beide versnellingen worden hier berekend gebruik makend van de nieuwe posities) en vervolgens
\(v_i(t+1) = v_i(t) + \frac{a_i(t) + a_i(t+1)}{2}\cdot \Delta t\)
\(v_j(t+1) = v_j(t) + \frac{a_j(t) + a_j(t+1)}{2}\cdot \Delta t\)
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?

Gebruikersavatar
Berichten: 305

Re: Probleem bij baansimulatie

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

Gebruikersavatar
Lorentziaan
Berichten: 1.433

Re: Probleem bij baansimulatie

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.

Gebruikersavatar
Berichten: 305

Re: Probleem bij baansimulatie

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?

Gebruikersavatar
Lorentziaan
Berichten: 1.433

Re: Probleem bij baansimulatie

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

Gebruikersavatar
Berichten: 305

Re: Probleem bij baansimulatie

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;


Gebruikersavatar
Lorentziaan
Berichten: 1.433

Re: Probleem bij baansimulatie

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.artcompsci.org/kali/vol/two_bod...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.edu/~dwhysong/prog/symp.cpp

(c++-implementatie van Yoshida's algoritmen tot en met orde 8)

Gebruikersavatar
Lorentziaan
Berichten: 1.433

Re: Probleem bij baansimulatie

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

Energiefouten voor Euler-integratie:

Code: Selecteer alles

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:

Code: Selecteer alles

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.

Gebruikersavatar
Berichten: 305

Re: Probleem bij baansimulatie

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.

Code: Selecteer alles

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

Code: Selecteer alles

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?

Gebruikersavatar
Lorentziaan
Berichten: 1.433

Re: Probleem bij baansimulatie

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.

Gebruikersavatar
Berichten: 305

Re: Probleem bij baansimulatie

Afbeelding

Gebruikersavatar
Moderator
Berichten: 51.259

Re: Probleem bij baansimulatie

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

Gebruikersavatar
Lorentziaan
Berichten: 1.433

Re: Probleem bij baansimulatie

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.

Reageer