Zon + planeet simulator; hoe?

[Archer]

Ik bereken inderdaad alles met de 'echte' eenheden, en daarna vertaal ik ze pas naar scherm coördinaten. Ik denk dat ik de formule nu wel begrijp. Maar dan heb ik toch nog een vraagje. Als ik met die formule die vector heb uitgerekend, kan ik dan de planeet verplaatsen met die cijfers van de vector, of dient er nog iets anders te gebeuren?

Want zoals ik het nu zie, zit er geen component van tijd of iets dergelijks in. Dus als ik de formule voor elke planeet 60 keer per seconde uit laat voeren, worden - volgens mij - altijd dezelfde waarden teruggegeven. En als hij elke keer - zoals in Brinx' voorbeeld - 3 x en 4 y gaat, dan verdwijnt die planeet snel van het scherm. Volgens mijn - vaak onnavolgbare logica - zou de formule toch ook ene richting of iets dergelijks moeten hebben? Als die planeet binnen geslingerd wordt van links, zou je andere waarden verwachten dan als hij binnen geslingerd wordt van rechts.

Laat het me alsjeblieft weten als ik onzin uitkraam, dan weet ik dat ook weer.

Archer

[Brinx]

Om die versnelling te gebruiken in je simulatie, zul je voor iedere tijdstap twee dingen moeten doen: de nieuwe snelheid van de planeet uitrekenen en de nieuwe positie. De meest eenvoudige manier is de volgende:

\(\vec{V}_{nieuw} = \vec{V}_{oud} + dt \cdot \vec{a}_{oud}\)

en

\(\vec{X}_{nieuw} = \vec{X}_{oud} + dt \cdot \vec{V}_{oud}\)

In deze vergelijkingen is 'dt' de tijdstap, 'a' de versnellingsvector, 'V' de snelheidsvector en 'X' de positievector. De versnelling die volgde uit de formule wordt dus alleen indirect gebruikt om de nieuwe positie te berekenen! Bovenstaande vergelijkingen vormen de meest simpele manier om numeriek te integreren (methode van Euler, zie http://en.wikipedia.org/wiki/Euler's_method ). Betere methodes zijn voorhanden (waarbij de numerieke fout veel minder snel toeneemt), zoals RK4 ( http://en.wikipedia.org/wiki/Runge-***_methods ) of nog veel hogere-orde methoden.

[aaargh]

Ik heb ooit eens een mooie simulator gemaakt. http://sciencetalk.nl/forum/index.php?showtopic=59454

Jammer genoeg wekt hij nu niet meer, alle planeten spiralen nu naar buiten. Ik kan de fout niet vinden, en ik ben nu ook bezig met andere projecten. (EM i.p.v gravitatie bijvoorbeeld)

[Bart]

Ik heb ooit eens een mooie simulator gemaakt. http://sciencetalk.nl/forum/index.php?showtopic=59454

Jammer genoeg wekt hij nu niet meer, alle planeten spiralen nu naar buiten. Ik kan de fout niet vinden, en ik ben nu ook bezig met andere projecten. (EM i.p.v gravitatie bijvoorbeeld)

klinkt alsof er ergens een minteken mist

[Brinx]

Ik ben op het moment ook zel aan het werken aan een eenvoudig baansimulatortje, waarbij de gebruiker via een grafische interface complete controle heeft over de positie en snelheid (ook tijdens het lopen van de simulatie).

De integratoren die ik tot nu toe heb ingebouwd zijn de 2e-orde leapfrog (symplectisch), de 4e-orde Runge-*** en de 4-orde Yoshida (symplectisch). Ik moet alleen een fout gemaakt hebben in de implementatie ervan, want alledrie de integratie-algoritmen gedragen zich als eerste-orde algoritmen! Met andere woorden, in alle drie de gevallen is de fout lineair afhankelijk van de gekozen tijdstap. Dat hoort uiteraard niet het geval te zijn!

Ik ben nu aan het uitzoeken waar deze fout vandaan komt. De integratiestappen lijken te kloppen, ook zit ik in het lineaire bereik van de integratoren qua gekozen tijdstap. De derde mogelijkheid, die ik nu moet uitzoeken, is die van de gekozen magnitude van de te integreren grootheden. Het is misschien namelijk mogelijk dat de afstanden in de simulatie (orde duizenden kilometers) leiden tot onvoldoende precisie bij bijvoorbeeld het uitrekenen van de optredende versnellingen.

Wordt vervolgd!

[aaargh]

Neenee, de planeten trekken nog normale banen, alleen gaan ze langzaamaan naar buiten, zelfs al draaien ze cirkelbaantjes, ze gaan steeds meer naar buiten. Mijn simulator heeft altijd wat verkeerd geweest, maar zo ongelooflijk sterk naar buiten gaand nooit.

Ik weet niet waardoor het komt, en het kan met toch niet zoveel schelen, ik ben nu bezig met een elektro(magnetisch)e simulatie.

Er is nog niet zoveel aan, maar zo ziet het er al uit.

[Brinx]

n.b.: ik krijg nu trouwens hetzelfde probleem als aaargh: planeetbanen spiralen bij mijn simulatirtje langzaam naar binnen! Dit komt doordat het lineaire gedrag van de simulator (euler-achtig gedrag) zorgt voor een geleidelijk verlies van de totale energie. Maar goed, dit is dus gedrag dat niet voor zou moeten komen bij symplectische integratoren!

[Bart]

Welke tijdstappen neem je voor je simulaties? Het kan zijn dat je tijdstap te groot is, waardoor je numerieke onnauwkeurigheden krijgt.

[aaargh]

Het onbegrijpbare is voor mij dat het vroeger perfect werkte, en pas na een aanpassing spiraliseerden ze naar buiten. Je moet maar eens kijken, ik zet het programma op het net. Wie erin geinterreseerd is, kijkt ernaar, ik kijk er niet meer naar om.

Download hier

Ik hoop dat het een beetje duidelijk is, anders zal ik nog wat uitleg erbij geven.

[Brinx]

Bart, de tijdstap is het probleem niet: ik zit in het 'lineaire bereik' van de integrator. Met andere woorden, er is een keurig nette lineaire relatie tussen de energiefout en de gebruikte tijdstap. Die energiefout kan ik op die manier in principe zo klein maken als ik wil (tot aan machineprecisie), maar het probleem is nu juist dat dit lineaire gedrag niet voor zou moeten komen. Voor een tweede-orde integrator geldt dat de energiefout 100 keer zo klein wordt wanneer je de tijdstap met een factor 10 verkleint, en voor een 4e-orde integrator zou die energiefout zelfs kleiner moeten worden met een factor 10000. Maar dat gebeurt dus niet...een raadsel!

Ik heb even naar je code gekeken aaargh, en ik zie dat je Euler-integratie gebruikt. In de functie 'Natuur()', in frmMain.frm.

Code: Selecteer alles

Private Sub Natuur()

Dim i As Integer

For i = 1 To obj

If b(i).fixed Then b(i).vx = 0: b(i).vy = 0

b(i).x = b(i).x + b(i).vx

b(i).y = b(i).y + b(i).vy

b(i).fax = b(i).ax

b(i).fay = b(i).ay

b(i).vx = b(i).vx + b(i).ax

b(i).vy = b(i).vy + b(i).ay

b(i).ax = 0

b(i).ay = 0

Next i

For i = 1 To obj

Call Grav(i)

Next i

End Sub

Dus: hier wordt de nieuwe positie gelijkgesteld aan de oude positie plus de oude snelheid, en de nieuwe snelheid is de oude snelheid plus de oude versnelling. Blijkbaar gebruik je impliciet een tijdstap van 1 seconde trouwens (omdat je v_niew = v_oud + a_oud gebruikt, en niet v_nieuw = v_oud + dt * a_oud)?

Euler integratie is eerste-orde, en het 'wegdrijven' van planeetbanen is te verwachten met zo'n integrator. De snelheid waarmee dat gebeurt is wel afhankelijk van de keuze van je parameters (massa's, afstanden, snelheden)!

Welke aanpassing had je gedaan vlak voordat de simulator dit gedrag ging vertonen?

[aaargh]

Dus dit noemt de Euler-methode? Ik heb het nochtans zelf gevonden.

Mijn laatste aanpassing was in de sub "Grav(i)". Ik maakte de zwaartekracht in 2 richtingen i.p.v. in 1 richting.

[Schwartz]

Men moet een type nemen die past bij het item:

10 bytes is echt nodig als men een goede berekening wilt.

De computer gaat bij 4 bytes integer als aardig de fout in bij het resultaat.

Na 1 berekening is het weining maar bij 200 per seconde maakt met na 1 minuut al 1200 keer die fout.

Omdat negatieve fouten en positieve fouten op elkaar weg worden gestreept moet men wel een gebalanceerde berekening maken.

Bij integer bewerking gaat de computer uit van een ongebalanceerde berekening.

Of men moet 0.499999999999 toevoegen maar dat kan alleen bij (real of extended in pascal).

Vroeger had ik een simpele simulatie gemaakt op de BBC acorn.

Ik had bij 20 sterretjes, al een aardig resultaat na 2 uur draaien.

Bij 2 sterren kon ik live zien hoe ze rond elkaar zaten te draaien.

Soms bleef het 1 minuut stabiel om toch uiteindelijk na wat gezwier instabiel te worden.

Men kan technieken toepassen om de getallen in de computer op een juist nivo te houden voor de berekening.

Bijvoorbeeld *100 voor de afstand in km en dan /100 voor de uiteindelijke optelling.

Ook kan men afwijken van een formule met een complexere berekening maar die technisch geheel gelijk is om de snelheid van de berekening hoog te houden (ook inzake de afwijking).

Niet elke berekening hoeft getoond te worden, men moet het aantal berekeningen instellen die men heeft bij 1 toning aan de gebruiker.