[C#]Omzetten wiskundige formule naar (pseudo) code

[CRiMiNaL]

Ik wil voor een hobby project het effect van zwaartekracht van de zon op planeten berekenen, op wikipedia vond ik de volgende formule (op basis van vectoren):

En omdat ik vroeger niet goed heb opgelet op school, lukt het mij niet om deze om te zetten naar werkende (pseudo)code in C#.

Wat ik tot nu toe heb staan:

C#:

class OrbitalBody
    {
        public const double GRAV = 6.6720e-08;

        public double mass;
        public Vector2 position;

        public OrbitalBody(double m, Vector2 p)
        {
            this.mass = m;
            this.position = p;
        }

        public Vector2 CalculateGravForce(OrbitalBody ob1)
        {
            double a = ob1.mass * this.mass;
            Vector2 distAB = this.position - ob1.position;

            Vector2 unitAB = distAB;
            unitAB.Normalize();

            double temp = a / (double)distAB.LengthSquared();
            Vector2 retVec = -(Vector2.Multiply(unitAB, (float)(temp * GRAV)) );
            return retVec;

        }

    }

Ik weet dat de fout in de laatste 2 regels zit. Ik heb de twee massa's met elkaar vermenigvuldigd en deel ze hier door de afstand tussen de 2 objecten in het kwadraat, maar ik vermoed dat ik de formule daar verkeerd begrijp. Iemand die mij hier een duwtje in de rug kan geven?

offtopic:
Ik weet dat ik een hoop precisie verlies door het gebruik van een op (float) gebaseerde Vector class en ben ook van plan dit te gaan wijzigen als ik het bovenstaande eenmaal werkend heb

[Woy]

Je kunt het toch gewoon debuggen?

Kijk eens per stap wat je verwacht dat er uit komt ( Met de rekenmachine ), en vergelijk dat in de debugger met wat er daadwerkelijk gebeurt.

[Armageddon_2k]

En probeer de namen van de formule terug te laten komen in je code.
Want nu (met variablen als: a, temp, abdist) is het hard onoverzichtelijk.

code:

 public Vector2 CalculateGravForce(OrbitalBody ob1) 
        { 
            double m1m2 = ob1.mass * this.mass; 
            Vector2 r12 = this.position - ob1.position; 
             r12 .Normalize(); 

            Vector2 r12_2 = distAB.LengthSquared(); 
        }

En probeer het samenvoegen van regels pas te doen als je code werkt.
Dit houdt het overzichtelijk en zorgt ervoor dat je je fouten veel sneller ziet.

[ Voor 14% gewijzigd door Armageddon_2k op 18-05-2010 15:43 ]

[CRiMiNaL]

Maar wat nu als ik de formule verkeerd interpreteer? Moet ik wel delen door lengte² van de vector? Ik ben hier dus niet zeker van en ditzelfde probleem heb ik ook op een rekenmachine. Dat maakt het misschien een meer wiskundig gerelateerd probleem, maargoed.

Ik was net blij dat de computer alle Vector math voor mij deed, maargoed, je hebt een punt, rekenmachine en papier en een paar uur aankloten kunnen misschien het e.e.a. verduidelijken.

---

Nogmaals de functie, met hopenlijk duidelijkere variabelen namen:

C#:

public Vector2 CalculateGravForce(OrbitalBody ob1)
        {
            double m1m2 = ob1.mass * this.mass;
            Vector2 r12 = this.position - ob1.position;

            Vector2 unit_r12 = r12;
            unit_r12.Normalize();

            double temp = m1m2 / (double)r12.LengthSquared();
            Vector2 F12 = -(Vector2.Multiply(unit_r12, (float)(temp * G)) );
            return F12;

        }

@ hieronder, dit was een reactie op Woy

[ Voor 38% gewijzigd door CRiMiNaL op 18-05-2010 15:54 ]

[Armageddon_2k]

Zoals ik al aangaf:
Probeer je code op te delen in meer regels.
Ik zie sowiso al een paar fouten

Je moet namelijk uiteindelijk:
-G * (Uitkomst deling mass/afstand enz.) * (Unit Vector Obj1 tot Obj2).

Loop je code maar eens stap voor stap door, wat je precies moet doen kan ik zo niet zeggen maar ik weet wel dat er minder code staat dan de formule hoort te zijn
Nogmaals, deel je code netjes op.

[Armageddon_2k]

Oke ik zie iets wat sowiso misgaat:

Vector2 F12 = -(Vector2.Multiply(unit_r12, (float)(temp * G)) );

je gaat hier vanuit dat unit_r12 gelijk is aan ^unit_r12.
Kijk maar even goed in de link die je stuurde.

[leuk_he]

die |r₁₂| wil zeggen dat je de absolute waarde neemt

(20 wordt 20)
-20 wordt 20
0 wordt 0 (maar delen door 0 levert een zwart gat op)


/Edit ik dacht dat normalize de vector op 1 zetten, maar goed, beetje debuggen zal meer helopen.
:
http://msdn.microsoft.com...ows.vector.normalize.aspx

[ Voor 37% gewijzigd door leuk_he op 18-05-2010 16:06 . Reden: normizle. ]

[Woy]

CRiMiNaL schreef op dinsdag 18 mei 2010 @ 15:45:
Maar wat nu als ik de formule verkeerd interpreteer? Moet ik wel delen door lengte² van de vector? Ik ben hier dus niet zeker van en ditzelfde probleem heb ik ook op een rekenmachine. Dat maakt het misschien een meer wiskundig gerelateerd probleem, maargoed.

Dat staat nou juist aardig duidelijk uitgelegd op de pagina die je aanhaalt

F12 is the force applied on object 2 due to object 1,
G is the gravitational constant,
m1 and m2 are respectively the masses of objects 1 and 2,
|r12| = |r2 − r1| is the distance between objects 1 and 2, and
is the unit vector from object 1 to 2.

Dus inderdaad F12 = -G * ( ( M1 * M2 ) / ( (r2-r1).LengthSquared ) ) * (r2-r1).UnitVector

leuk_he schreef op dinsdag 18 mei 2010 @ 15:54:
die |r₁₂| wil zeggen dat je de absolute waarde neemt

Nee dat betekent dat het de lengte van de vector is, en een ABS zou niks uithalen aangezien je hem daarna toch kwadrateert, en dus sowieso positief is.

[ Voor 17% gewijzigd door Woy op 18-05-2010 15:56 ]

[Armageddon_2k]

leuk_he schreef op dinsdag 18 mei 2010 @ 15:54:
die |r₁₂| wil zeggen dat je de absolute waarde neemt

(20 wordt 20)
-20 wordt 20
0 wordt 0 (maar delen door 0 levert een zwart gat op)

Nope das een normalize,
En dat heeftie met Normalize() al gefixed.

Woy!, nu lijk ik weer traag.

[ Voor 8% gewijzigd door Armageddon_2k op 18-05-2010 15:58 ]

[Woy]

Armageddon_2k schreef op dinsdag 18 mei 2010 @ 15:56:
[...]

Nope das een normalize,
En dat heeftie met Normalize() al gefixed.

Woy!, nu lijk ik weer traag.

Maar het klopt ook niet wat je zegt, want dat is geen Normalize. Normalize is schalen met de lengte van de vector, zodat je een vector van lengte 1 krijgt.

[Armageddon_2k]

Woy schreef op dinsdag 18 mei 2010 @ 16:12:
[...]
Maar het klopt ook niet wat je zegt, want dat is geen Normalize. Normalize is schalen met de lengte van de vector, zodat je een vector van lengte 1 krijgt.

My bad, je hebt gelijk, was even geconfuseerd.
Die hele normalize kan er uit, gaat sowiso rare resultaten geven.
Desalniettemin wordt in de code de hele unit vector vergeten.

[Woy]

Armageddon_2k schreef op dinsdag 18 mei 2010 @ 16:18:
[...]

My bad, je hebt gelijk, was even geconfuseerd.
Die hele normalize kan er uit, gaat sowiso rare resultaten geven.
Desalniettemin wordt in de code de hele unit vector vergeten.

Nee die Normalize is wel nodig, want dat is de richtingsvector van de kracht, die vermenigvuldig je dan met de schaal van de kracht ( G*M1*M2/|r12|² ).

De Normalize levert juist die unit vector op

[ Voor 6% gewijzigd door Woy op 18-05-2010 16:21 ]

[CRiMiNaL]

Woy schreef op dinsdag 18 mei 2010 @ 16:12:
[...]

Maar het klopt ook niet wat je zegt, want dat is geen Normalize. Normalize is schalen met de lengte van de vector, zodat je een vector van lengte 1 krijgt.

Als je een Vector normalized krijg je een Unit Vector en een UnitVector heeft een lengte van 1, dus ik begrijp het verschil dat jij tracht aan te duiden niet meer.

Daarnaast word er op wikipedia ook aangegeven dat || gebruikt word om de lengte aan te geven en | voor een absolute waarde. Maar | word ook af en toe gebruikt om de lengte van de Vector aan te geven. Het is dat de absolute waarde nemen hier geen nut heeft, omdat het in het kwadraat word gezet.

Na de suggesties hier is mijn code inmiddels veranderd tot onderstaande:

C#:

 public Vector2 CalculateGravForce(OrbitalBody ob1)
        {
            double m1m2 = ob1.mass * this.mass;
            Vector2 r12 = this.position - ob1.position;

            Vector2 unit_r12 = r12;
            unit_r12.Normalize();


            Vector2 v_unit_r12 = r12 / unit_r12;

            double temp = m1m2 / (double)r12.LengthSquared();

            Vector2 F12 = -(float)G * (float)temp * v_unit_r12;
            
            //Vector2 F12 = -(Vector2.Multiply(unit_r12, (float)(temp * G)) );
            return F12;

        }

[.oisyn]

Ik zie geen fout in de code. Waar Armageddon en Woy het allemaal over hebben ontgaat me compleet.

[ Voor 4% gewijzigd door .oisyn op 18-05-2010 16:26 ]

[Woy]

CRiMiNaL schreef op dinsdag 18 mei 2010 @ 16:22:
[...]
Als je een Vector normalized krijg je een Unit Vector en een UnitVector heeft een lengte van 1, dus ik begrijp het verschil dat jij tracht aan te duiden niet meer.

Ik probeer geen verschil uit te leggen. Armageddon_2k reageerde op de post van leuk_he en zei dat dat een normalize was, wat niet klopte.

Het delen door de lengte, en dus een Unit vector maken is inderdaad Normaliseren.

Daarnaast word er op wikipedia ook aangegeven dat || gebruikt word om de lengte aan te geven en | voor een absolute waarde. Maar | word ook af en toe gebruikt om de lengte van de Vector aan te geven. Het is dat de absolute waarde nemen hier geen nut heeft, omdat het in het kwadraat word gezet.

Wiskundig gezien kan het inderdaad ook gebruikt worden als absolute waarde, maar in dit geval betekent het gewoon lengte.

Na de suggesties hier is mijn code inmiddels veranderd tot onderstaande:

Maar wat is de reden dat je denkt dat het niet goed is? Komt er wat anders uit dan je denkt?

.oisyn schreef op dinsdag 18 mei 2010 @ 16:25:
Ik zie geen fout in de code. Waar Armageddon en Woy het allemaal over hebben ontgaat me compleet.

Dat was als reactie op leuk_he die stelde dat |r12| betekent dat je de absolute waarde van r12 neemt, want hier dus niet het geval is.

Overigens zie ik nog steeds niet wat er nou fout zou moeten zijn aan je initiële code, en waarom je nu in je code de v_unit_r12 geïntroduceerd hebt.

[ Voor 16% gewijzigd door Woy op 18-05-2010 16:33 ]

[Armageddon_2k]

Woy schreef op dinsdag 18 mei 2010 @ 16:20:
[...]

Nee die Normalize is wel nodig, want dat is de richtingsvector van de kracht, die vermenigvuldig je dan met de schaal van de kracht ( G*M1*M2/|r12|² ).

De Normalize levert juist die unit vector op

Ik verdwaal helemaal in dat 1e stuk code, ik had geinterperteerd dat de normailze al voor de deling van M1M2 werd gedaan. Even de code opnieuw gelezen

Comment op de nieuwe code:

code:

           unit_r12.Normalize(); 
            Vector2 v_unit_r12 = r12 / unit_r12;

Is die deling nog nodig? wat doet de normalize functie, als ik em zo even snel doorlees is:
Ik weet het niet zeker, hangt er even vanaf wat c# ermee doet.
"The term normalized vector is sometimes used as a synonym for unit vector."

dat zou dan betekenen dat je die deling weg kan laten.

Maar!!!!
.Normalize();
Past deze je vector aan? of levert deze een vector?
bv: String.Replace() laat de orginele string ook intact, en geeft een nieuwe terug.

code:

           Vector2 v_unit_r12 = unit_r12.Normalize();

[ Voor 41% gewijzigd door Armageddon_2k op 18-05-2010 16:35 ]

[Woy]

Armageddon_2k schreef op dinsdag 18 mei 2010 @ 16:27:
[...]

Maar!!!!
.Normalize();
Past deze je vector aan? of levert deze een vector?
bv: String.Replace() laat de orginele string ook intact, en geeft een nieuwe terug.

code:

           Vector2 v_unit_r12 = unit_r12.Normalize();

Ik ga ervan uit dat hij de Vector2 uit het XNA framework gebruikt en die past inderdaad de huidige vector aan, dus dat zou geen probleem moeten zijn
http://msdn.microsoft.com...(v=XNAGameStudio.31).aspx

[ Voor 9% gewijzigd door Woy op 18-05-2010 16:38 ]

[.oisyn]

CRiMiNaL schreef op dinsdag 18 mei 2010 @ 16:22:
Na de suggesties hier is mijn code inmiddels veranderd tot onderstaande:

De suggesties zijn suf, je originele code uit de startpost was gewoon goed. Wat je nu hebt gemaakt is alleen maar onduidelijker. Daarnaast deel je een vector door een andere vector, dat kan wiskundig niet eens (compileert dat überhaupt?)

[Woy]

.oisyn schreef op dinsdag 18 mei 2010 @ 16:39:
[...]

De suggesties zijn suf, je originele code uit de startpost was gewoon goed. Wat je nu hebt gemaakt is alleen maar onduidelijker. Daarnaast deel je een vector door een andere vector, dat kan wiskundig niet eens (compileert dat überhaupt?)

Ik zie in die wijziging anders niks van de suggesties terug, want de suggestie was dat hij eens moest debuggen en kijken waar zijn code een ander resultaat had dan hij verwacht.

Overigens zal dat denk ik wel compileren: http://msdn.microsoft.com...(v=XNAGameStudio.31).aspx
Het deelt blijkbaar gewoon de componenten door elkaar. Maar het nut ontgaat me ook een beetje.

[ Voor 19% gewijzigd door Woy op 18-05-2010 16:46 ]

[.oisyn]

Ik zie in die wijziging anders niks van de suggesties terug

Ik wel, de wijziging aan de code is namelijk net zo onsamenhangend als bepaalde comments in deze topic .

Anyway, bottom line is, die code is goed, het probleem zal wel ergens anders liggen. Zo heb ik de TS nog altijd niet horen noemen waarom hij denkt dat de code verkeerd is, oftewel wat de symptomen zijn van wat hij ziet, en wat hij verwacht te zien (hey, waren dit geen punten uit de quickstart? ).

Overigens ben ik wel benieuwd wat de TS verder met deze berekening wil. Het is namelijk een momentopname. Wil je de boel animeren dan zul je een numerieke integrator moeten gebruiken. Als je gewoon elke frame de dingen uitrekent en die simpelweg gebruikt als kracht-vectoren dan valt je planetenstelsel uit elkaar (of hij stort in). Een stabiele orbit zul je iig niet krijgen

Woy schreef op dinsdag 18 mei 2010 @ 16:43:
Overigens zal dat denk ik wel compileren: http://msdn.microsoft.com...(v=XNAGameStudio.31).aspx
Het deelt blijkbaar gewoon de componenten door elkaar. Maar het nut ontgaat me ook een beetje.

Ah. Brrr. Doet HLSL ook (waar het overigens wat practischer is). Ik zie a*b liever als dotproduct dan als een component-multiply. De Vector klassen van managed Direct3D doen dit overigens niet (ook geen dotproduct trouwens, je kunt een vector alleen met een scalar vermenigvuldigen en andersom), vreemd

[ Voor 55% gewijzigd door .oisyn op 18-05-2010 16:53 ]

[CRiMiNaL]

Die deling was inderdaad niet nodig, het was het resultaat van mijn verwarring over het unit vector verhaal, maar gewoon vermenigvuldigen met de unit vector is hier voldoende.

Mijn reden tot posten was dat de uitkomst van mijn initiële code een vector gaf met een langere lengte dan de afstand tussen de twee.

C#:

OrbitalBody sun = new OrbitalBody(1.9891e30, new Vector2(0, 0));
OrbitalBody earth = new OrbitalBody(5.9736e24, new Vector2(147098290000, 0));

Eerste argument is massa in kg, tweede is de positie vector (Met afstand tussen zon - aarde in meters voor x).
Dus vector zon - aarde is = {X:1,470983E+11 Y:0}
Vector van de zwaartekracht op de aarde volgens de huidige berekening: {X:-3,663815E+25 Y:0}

Ik verwachte een vector met een kortere lengte, die als hij vermenigvuldigd werd met de positie & "velocity" resulteerde in een effect waarbij je een orbit om de zon krijgt. Maar het is mij inmiddels duidelijk dat ik de materie nogal onderschat heb en me nog iets meer in de achterliggende wiskunde moet verdiepen.

---

offtopic:
Ik besef me nu ook dat ik bovenstaande al had moeten vermelden in de TS en dat door het niet vermelden van de informatie het topic nogal een zooitje geworden is, mijn excuses

[ Voor 8% gewijzigd door CRiMiNaL op 18-05-2010 16:57 ]

[.oisyn]

CRiMiNaL schreef op dinsdag 18 mei 2010 @ 16:50:
Mijn reden tot posten was dat de uitkomst van mijn initiële code een vector gaf met een langere lengte dan de afstand tussen de twee

Nou en? De resulterende vector is een kracht. De vector wordt dus langer als de kracht toeneemt. De kracht neemt tevens toe als de afstand afneemt. Er is dus, per definitie, naarmate de afstand kleiner wordt, een moment waarbij de lengte van de krachtvector langer wordt dan de afstand.

Ik verwachte een vector met een kortere lengte, die als hij vermenigvuldigd werd met de positie & "velocity" resulteerde in een effect waarbij je een orbit om de zon krijgt.

Hmm, hoe je dit dan doet is me niet helemaal duidelijk, maar wat je wilt is een kracht omzetten in een versnelling. Via de beroemde F=m∙a kom je op a=F/m, ergo, je moet de krachtvector eerst delen door de massa van het object, en dan krijg je z'n versnelling. Die versnelling moet je vervolgens weer verwerken in de snelheid van het object (met simpele euler integration vermenigvuldig je de versnelling gewoon met de timestep (constante bij een vaste framerate, of tijd tussen huidige en vorig frame bij variabele framerate) en tel je dat bij de snelheid op - maar Euler's method is net zo instabiel als dat hij simpel is ).

[ Voor 50% gewijzigd door .oisyn op 18-05-2010 17:08 ]

[G33rt]

.oisyn schreef op dinsdag 18 mei 2010 @ 16:56:
Nou en? De resulterende vector is een kracht. De vector wordt dus langer als de kracht toeneemt. De kracht neemt tevens toe als de afstand afneemt. Er is dus, per definitie, naarmate de afstand kleiner wordt, een moment waarbij de lengte van de krachtvector langer wordt dan de afstand.

Het vervelende van dit soort code als je het niet gewend bent (ik weet dat jij dat wel bent .oisyn) is dat het gewoon keiharde fysica blijft en al die grootheden eenheden hebben. Je bent in feite (weliswaar statische) computationele fysica aan het doen.

Om de verwarring voor de TS te verduidelijken: het is net zoiets als zeggen (artificieel voorbeeld; denk aan een kubus of een cilinder met grondoppervlak)

C++:

double area = 1457.8;
double volume = 50.3;

en dan de gedachte hebben dat het niet kan, omdat je een groter oppervlak dan volume hebt. Je vergeet echter dat je m³ helemaal niet kan vergelijken met m²!

In dit geval ben je verbaasd dat de kracht - in N - groter is dan de afstand - in m! Die twee getallen zijn echter totaal niet aan elkaar te relateren op liniaal (dwz: groter of kleiner dan elkaar)!

.oisyn schreef op dinsdag 18 mei 2010 @ 16:56:
Hmm, hoe je dit dan doet is me niet helemaal duidelijk, maar wat je wilt is een kracht omzetten in een versnelling. Via de beroemde F=m∙a kom je op a=F/m, ergo, je moet de krachtvector eerst delen door de massa van het object, en dan krijg je z'n versnelling. Die versnelling moet je vervolgens weer verwerken in de snelheid van het object (met simpele euler integration vermenigvuldig je de versnelling gewoon met de timestep (constante bij een vaste framerate, of tijd tussen huidige en vorig frame bij variabele framerate) en tel je dat bij de snelheid op - maar Euler's method is net zo instabiel als dat hij simpel is ).

Met een simpele stap als Eulerintegratie vermoed ik dat dat haast wel goed moet gaan (in ieder geval voor 1 enkel rondje): je hebt geen leipe dempingseffecten of zoiets (de differentiaalverlijking bevat geen eerste afgeleide). Een kubus op een aan de muur vastzittende veer schieten is bijvoorbeeld notoir omdat die veer we voor zo'n term (en dus demping) zorgt, waardoor je als je niet oppast dwars door de muur heen knalt (wanneer je blokje niet goed stil staat bij de veer zal zijn snelheid niet omkeren). In dit geval draait een object in de ruimte rondjes onder invloed van een zwaartekrachtsveld dus levert een klein discretisatiefoutje niet in een keer (na een omloop) zo'n overduidelijk fout effect als het blokje door de muur op - tenzij je hele groffe benaderingen doet.

Het enige probleem is dat als je niet helemaal uitkomt op je beginpunt je een debiele spiraal kunt programmeren als je niet oplet (in het geval dat je bij elke omloop consequent een stapje verder naar buiten of binnen zou komen), maar vermoedelijk kun je dat met hogere order termen (zie bijvoorbeeld Wikipedia: Runge-Kuttamethode) al grotendeels vermijden.

Wat je in feite wil doen is een setje differentiaalvergelijkingen numeriek oplossen, maar ik ben bang dat de meeste literatuur redelijk hoog gegrepen zal zijn als je vectoren e.d. al lastig vind - er komt een redelijke hoeveelheid klassieke mechanica bij om de hoek. Via Wikipedia: Computational mechanics kun je in ieder geval wel een eindje komen en wat links vinden.

[ Voor 50% gewijzigd door G33rt op 18-05-2010 18:36 ]

[NickThissen]

Ik heb in mijn opleiding wel een keer een orbit moeten simuleren. We lostten het stelsel op met de runge-kutta 45 methode, waarmee hij, met de juiste begincondities (tot op zo'n 8 a 9 cijfers achter de komma!) een stabiele ellips kon vormen. Een iets afwijkende beginconditie zorgt al voor een zichtbare afwijking na een paar rondjes.

Ik weet niet zeker of de Euler methode ook gaat werken; aangezien het juist de bedoeling was om bekend te raken met de RK-45 methode heb ik de Euler methode niet gebruikt. Maar het valt te proberen.

Ik kan de opdracht wat ik bedoel wel posten, met mijn oplossingen, maar er zit een hoop fysica in en vrij weinig programmeer werk (de RK-45 methode kregen we gewoon kado), en de code is in matlab gemaakt dus niet echt te vergelijken met C#.

Ook ging het bij onze opdracht niet om verschillende massa's ofzo, alles werd massaloos uitgedrukt, dus hiermee zul je niet verschillende planeten kunnen simuleren.


Meer on-topic: ik heb niet het hele topic gelezen, maar ik zie wel wat rare dingen:
- de 'absoluut strepen' om een vector betekenen nooit absolute waarde; dat bestaat niet eens voor een vector. Voor een vector betekent het gewoon de lengte van de vector. En dubbele strepen worden ook wel eens gebruikt, is gewoon een andere notatie maar betekent hetzelfde.
- je deelt een vector door een andere vector. Dat kan wiskundig niet, dus ik vind het vreemd dat je code dit toelaat. Is Vector2 een klasse van jezelf, of van het XNA framework oid? Kun je eens uitzoeken wat Vector2 / Vector2 precies doet? Waarschijnlijk components-gewijs delen, wat zeker niet de bedoeling is.


Ten slotte raad ik je aan de formule op een andere manier te bekijken. Het feit dat hij in vectoren is uitgedrukt maakt het alleen maar verwarrend als je daar niet bekend mee bent.

Er staat feitelijk dat de grootte van de kracht van massa 1 op massa 2 gelijk is aan (G m1 m2) / r^2, waar m1 en m2 de massa's zijn, r de afstand tussen de twee massa's (GEEN VECTOR!), en G de gravitatieconstante.
De rest (het minteken en de vector met het hoedje) zeggen feitelijk alleen dat deze kracht gericht is van massa 2 naar massa 1. Want, de vector met het hoedje is een eenheidsvector (dat geeft het hoedje aan) en deze heeft dus lengte 1. Het is een vector van lengte 1 die van massa 1 naar massa 2 wijst. Met het minteken is het dus een vector van lengte 1 die van massa 2 naar massa 1 wijst.
Omdat deze vector lengte 1 heeft hoef je hem eigenlijk helemaal niet te gebruiken in je berekening; hij verandert alleen de richting van de kracht en niet de grootte. Maar de richting weet je al: die is gewoon tegengesteld aan de vector van massa 1 naar massa 2.


Dus in feite kun je het in twee delen splitsen:

- De grootte is F12 = G m1 m2 / r^2
- De richting is van massa 2 naar massa 1. Dit kun je ook logisch beredeneren: de kracht die massa 1 op massa 2 uitoefent is aantrekkend (zwaartekracht is altijd aantrekkend) dus van massa 2 naar massa 1 gericht. Let wel dat massa 2 ook dezelfde kracht op massa 1 uitoefent (maar dan precies de andere kant op). Dus massa 1 (de zon) zal ook richting massa 2 gaan bewegen. Omdat de massa van de zon echter zoveel groter is dan de massa van een planeet, zal de zon veel langzamer gaan bewegen, en vaak wordt (bij benadering) aangenomen dat de zon niet beweegt, wat de berekening een stuk makkelijker maakt.

[ Voor 7% gewijzigd door NickThissen op 18-05-2010 19:20 ]

[.oisyn]

G33rt schreef op dinsdag 18 mei 2010 @ 18:13:
Met een simpele stap als Eulerintegratie vermoed ik dat dat haast wel goed moet gaan (in ieder geval voor 1 enkel rondje): je hebt geen leipe dempingseffecten of zoiets (de differentiaalverlijking bevat geen eerste afgeleide).

Het probleem met orbits op basis van de Euler methode is dat de zwaartekracht te weinig wordt toegepast. Als we even de aarde en de maan nemen, dan beweegt de maan zich steeds weg van de aarde. De zwaartekracht trekt de maan continu naar de aarde toe. Echter, met een vaste timestep beweeg je de maan een stukje opzij, maar in al die tijd is er geen kracht die loodrecht op de beweging zit - die komt er de volgende timestep pas bij, maar dan is het al te laat en is de maan een stukje verder van de aarde dan ie eerst was. Dit stapelt zich op en zo verliest de maan z'n baan om de aarde.

NickThissen schreef op dinsdag 18 mei 2010 @ 19:10:
Er staat feitelijk dat de grootte van de kracht van massa 1 op massa 2 gelijk is aan (G m1 m2) / r^2, waar m1 en m2 de massa's zijn, r de afstand tussen de twee massa's (GEEN VECTOR!)

Het grappige aan die vergelijking is dat het ook werkt als r wél een vector is . r² is namelijk een dotproduct met zichzelf, en dus lengte-kwadraat.

[ Voor 19% gewijzigd door .oisyn op 18-05-2010 19:36 ]

[NickThissen]

.oisyn schreef op dinsdag 18 mei 2010 @ 19:34:
[...]

Het grappige aan die vergelijking is dat het ook werkt als r wél een vector is . r² is namelijk een dotproduct met zichzelf, en dus lengte-kwadraat.

Dat snap ik ook wel, maar als de OP niet veel verstand van vectoren heeft kan hij ze beter vergeten en gewoon grootte en richting apart behandelen. In het geval van zwaartekracht is dat triviaal, omdat je gewoon 'common sense' kunt gebruiken om te weten welke richting de kracht moet hebben. Aan zijn code te zien weet hij hoe hij de vector van massa 1 naar massa 2 krijgt (r₁₂). De richting van F₁₂ is dan gewoon de andere kant op. De grootte krijgt hij door alle vector zooi te vergeten in de vergelijking (behalve dan LengthSquared van r₁₂ maar die had hij al).

[ Voor 21% gewijzigd door NickThissen op 18-05-2010 19:50 ]

[G33rt]

.oisyn schreef op dinsdag 18 mei 2010 @ 19:34:
[...]

Het probleem met orbits op basis van de Euler methode is dat de zwaartekracht te weinig wordt toegepast. Als we even de aarde en de maan nemen, dan beweegt de maan zich steeds weg van de aarde. De zwaartekracht trekt de maan continu naar de aarde toe. Echter, met een vaste timestep beweeg je de maan een stukje opzij, maar in al die tijd is er geen kracht die loodrecht op de beweging zit - die komt er de volgende timestep pas bij, maar dan is het al te laat en is de maan een stukje verder van de aarde dan ie eerst was. Dit stapelt zich op en zo verliest de maan z'n baan om de aarde.

Ik kan me goed vergissen, maar in dat geval heb je op t₀ de zwaartekracht niet meegenomen toch? Daardoor zou je je maan inderdaad horizontaal zien verschuiven. Er is ook een a₀ en als je die wel meerekent zie ik niet direct de reden waarom het fout zou moeten gaan. Alhoewel ik dit probleem niet in concreto uitgeprogrammeerd (wel andere differentiaalvergelijkingen, en niet met de Eulermethode) heb, dus ben wel benieuwd wat ik dan hier over het hoofd zou zien.

[.oisyn]

Het probleem is dat het een tweede orde afgeleide is, en Euler een 'drag' heeft per orde.
x_i+1 = x_i + v_i∙dt
v_i+1 = v_i + a_i∙dt
Wat je a₀ ook is, hij wordt niet in x₁ verwerkt, maar pas in x₂. De orbit wordt op die manier steeds groter.

Overigens is hij wel redelijk stabiel als je de boel omdraait (eerst de snelheid updaten, daarna de positie). Wiskundig kom je dan op
v_i+1 = v_i + a_i∙dt
x_i+1 = x_i + v_i+1∙dt
Oftewel
x_i+1 = x_i + v_i∙dt + a_i∙dt²
v_i+1 = v_i + a_i∙dt

De orbit is dan wel stabiel, maar strict gezien is dit niet meer Euler.

[ Voor 60% gewijzigd door .oisyn op 19-05-2010 12:17 ]

[CRiMiNaL]

.oisyn schreef op woensdag 19 mei 2010 @ 11:50:
[......]
Oftewel
x_i+1 = x_i + v_i∙dt + a_i∙dt²
v_i+1 = v_i + a_i∙dt

De orbit is dan wel stabiel, maar strict gezien is dit niet meer Euler.

Dit begrijp ik niet helemaal, de lengte van vector v_i+1 word op de door jou beschreven manier toch alleen maar langer? (vertaald naar velocity: het object gaat steeds sneller).
Besef me net dat a_i natuurlijk ook negatief kan zijn.

En ik neem aan dat je met a_i∙dt bedoeld de acceleratie te vermenigvuldigen met het tijdsverschil (delta Time) ?

[.oisyn]

idd
Zie trouwens http://oisyn.nl/orbit.html voor een voorbeeld van die laatste methode. Het hele systeem neigt naar verloop van tijd wat naar rechts door de door mij gegeven v₀ snelheden van beide objecten, maar op zich is wel te zien dat dit gewoon stabiel blijft. Implementeer je het op de euler manier, dan zul je zien dat de orbit steeds langer wordt. Ik wil nog even een versie maken waarin je euler, mijn aanpassing daarop, en RK4 met elkaar kunt vergelijken.

Deze pagina is voor jou wellicht nog interessant om te lezen, daar staat RK4 ook meteen uitgelegd: http://gafferongames.com/game-physics/integration-basics/

[ Voor 162% gewijzigd door .oisyn op 19-05-2010 14:33 ]

[G33rt]

.oisyn schreef op woensdag 19 mei 2010 @ 11:50:
x_i+1 = x_i + v_i∙dt + a_i∙dt²
v_i+1 = v_i + a_i∙dt

De orbit is dan wel stabiel, maar strict gezien is dit niet meer Euler.

Ah, daar komt mijn verwarring vandaan - ik zou onmiddelijk met deze twee vergelijkingen gaan werken (het is immers een 2^e-orde differentiaalvergelijking), zonder er bij stil te staan dat je dan strikt gesproken niet meer aan het Eulerintegreren bent. Dank voor het verhelderen

[.oisyn]

.oisyn schreef op woensdag 19 mei 2010 @ 14:31:
idd
Zie trouwens http://oisyn.nl/orbit.html voor een voorbeeld van die laatste methode. Het hele systeem neigt naar verloop van tijd wat naar rechts door de door mij gegeven v₀ snelheden van beide objecten, maar op zich is wel te zien dat dit gewoon stabiel blijft. Implementeer je het op de euler manier, dan zul je zien dat de orbit steeds langer wordt. Ik wil nog even een versie maken waarin je euler, mijn aanpassing daarop, en RK4 met elkaar kunt vergelijken.

Done. Op die pagina zie je nu 3 orbits, de rode is de euler methode en loopt al vrij snel de soep in. De groene is wat ik maar even "reverse" euler heb genoemd, waarbij de snelheid eerst wordt geupate en die weer wordt gebruikt om de positie te updaten. De blauwe is een RK4 implementatie. Groen en blauw lopen redelijk met elkaar op, maar toch is te zien dat de groene iets langzamer gaat.

[CRiMiNaL]

Er is mij in ieder geval een hoop duidelijk geworden, graag wil ik mijn conclusie nog even met jullie verifiëren. Relevante code:

C#:

            bool test = false;
            if (test)
            {
                OrbitalBody.G = 10000;
                OrbitalBody.scaler = 1;
                speed = 5;
                //                    mass,      position,          velocity   
                sun = new OrbitalBody(100, new Vector(0, 0), new Vector(0,0));
                earth = new OrbitalBody(1, new Vector(0, 200), new Vector(75, 0));
            }
            else
            {
                OrbitalBody.G = 6.6720e-08;
                OrbitalBody.scaler = 1000000000;
                speed = 5000;
                sun = new OrbitalBody(1.9891e30, new Vector(0, 0), new Vector(0, 0));
                earth = new OrbitalBody(5.9736e24, new Vector(0, 149598261000), new Vector(2.978589e4, 0));
            }

Draai ik de simulatie met test = true, dan krijg ik een werkende orbit (waardes gestolen van .oisyn's javascript)
Pak ik echter de test = false, dan stort de aarde al snel naar de zon, terwijl de waardes voor massa / velocity / positie allemaal lijken te kloppen. Mijn conclusie is dat dit komt door de hoge waarde van de variabele speed, maargoed, met speed = 1 betekend het dus gewoon mooi dat ik een jaar mag wachten op een complete orbit. Is deze conclusie juist, of klopt een van mijn parameters niet? (Berekeningen zijn verder identiek aan die van .oisyn)

Bedankt voor alle uitleg en links trouwens.

edit:
de "scalar" variabel gebruik ik alleen om de uiteindelijke positie binnen de schermcoordinaten te krijgen en word verder niet gebruikt bij de berekeningen

[ Voor 8% gewijzigd door CRiMiNaL op 19-05-2010 17:05 ]

[.oisyn]

speed=5000 geeft denk ik idd een te grote timestep, dan performt RK4 waarschijnljik wat beter. Wat je denk ik beter kan doen (ik neem aan dat je die reverse euler implementatie gebruikt?) is de speed omlaag zetten, maar gewoon meerdere keren per frame de simulatiestap te doen (bijv. speed=50 en dan 100 simulatiestappen per render). Dat verhoogt de precisie van het algoritme.

Overigens zijn 2 bodies ook nog analytisch uit te rekenen, wellicht wil je daar ook wat mee doen zodat je referentiemateriaal hebt.

.edit: als je trouwens in mijn test het grote object een beginsnelheid van (-0.75, 0) meegeeft dan blijft het geheel op z'n plek.

[ Voor 26% gewijzigd door .oisyn op 19-05-2010 17:11 ]

[NickThissen]

.oisyn schreef op woensdag 19 mei 2010 @ 16:44:
[...]


Done. Op die pagina zie je nu 3 orbits, de rode is de euler methode en loopt al vrij snel de soep in. De groene is wat ik maar even "reverse" euler heb genoemd, waarbij de snelheid eerst wordt geupate en die weer wordt gebruikt om de positie te updaten. De blauwe is een RK4 implementatie. Groen en blauw lopen redelijk met elkaar op, maar toch is te zien dat de groene iets langzamer gaat.

Leuk weergegeven. Misschien is het goed om eens te kijken wat er met de energie en het impulsmoment gebeurt tijdens de orbit. Die zouden constant moeten zijn maar door errors (verwaarlozing hogere orde termen en in mindere mate afrond fouten) zul je zien dat ze fluctueren en ook nog langzaam toe of afnemen. Hoeveel ze toe of afnemen is een indicatie van de fout die je maakt in de berekening. Ik neem aan dat je dan goed kan zien dat de RK4 methode het beste is.

De energie is
E = (1/2) ( v_x² + v_y² ) - 1 / sqrt( x² + y² )

Impulsmoment is
L = x v_y - y v_x

Zou je makkelijk elke stap moeten kunnen berekenen denk ik. Ik zou het zelf doen maar ik weet niet veel van javascript en niet eens hoe ik het zou moeten runnen ook al weet ik je code aan te passen

[.oisyn]

Goed idee, dat zal ik vanavond eens inbouwen.

Ik zou het zelf doen maar ik weet niet veel van javascript en niet eens hoe ik het zou moeten runnen ook al weet ik je code aan te passen

gewoon een kwestie van de .html aanklikken zodat ie in je browser opent

[ Voor 14% gewijzigd door .oisyn op 19-05-2010 21:09 ]

[NickThissen]

.oisyn schreef op woensdag 19 mei 2010 @ 19:21:
[...]

gewoon een kwestie van de .html aanklikken zodat ie in je browser opent

Ok dat was makkelijker dan gedacht Hij laat alleen de beginposities altijd zien, bovenop de rond draaiende vierkantjes.

Maar dan nog weet ik niet hoe ik de uitkomst zou moeten tonen... Weet niet eens hoe ik gewoon elke stap de waarde van E zou kunnen printen (heb gewoon "print" geprobeerd maar dat leek niet te werken), laat staan een mooi tabelletje (E als functie van aantal stappen / tijd ofzo).

[.oisyn]

Overigens

NickThissen schreef op woensdag 19 mei 2010 @ 19:12:
[...]

Leuk weergegeven. Misschien is het goed om eens te kijken wat er met de energie en het impulsmoment gebeurt tijdens de orbit. Die zouden constant moeten zijn

Volgens mij gaat dat alleen op als de orbit perfect circulair is. Dat is in mijn testgeval iig niet het geval, ik heb het kleinere object gewoon simpelweg ergens onder het grotere object laten beginnen en 'm een willekeurige zijwaartse snelheid meegegeven. De orbit lijkt redelijk circulair, maar dat is omdat ik de waarden zo heb getweakt, maar in werkelijkheid is het waarschijnlijk wel meer een ellips. In het extreme geval is de beginsnelheid van beide objecten 0 en krijg je een slingereffect. Het is logisch dat er dan geen impulsmoment is, en dat de kinetische energie bij de uiteinden 0 is en maximaal als ze door elkaar heen gaan.

[NickThissen]

Nee, de energie (niet alleen de kinetische energie, ook de potentiele energie erbij!) en het impulsmoment zijn altijd constant in de tijd, ook bij ellipsbanen en ik geloof zelfs ook bij een hyperbolische "baan". Schrijf de bewegingsvergelijkingen maar eens op en vul E en L dan in, je zult zien dat die niet van de tijd afhangen.

EDIT
Ik had wel de massa vergeten mee te nemen in mijn post.

Zal het even uittypen (misschien niet zo makkelijk te lezen). De massaloze bewegingsvergelijkingen zijn:

x' ' = - x / r³,
y' ' = -y / r³

Zoals ik zei is E de kinetische energie (1/2 |v|²) plus de potentiele energie (-1/r):
E = [ ( x' )² + ( y' )² ] / 2 - 1 / sqrt( x² + y² )

De tijdsafgeleide (een accentje is afleiden naar de tijd uiteraard) is dus
E' = x' x' ' + y' y' ' + r' / r²

met r' = dr/dx x' + dr/dy y' = x x' / r + y y' / r

In totaal is E' dus
E' = x' x' ' + y' y' ' + (x x' + y y') / r³
= - (x x' + y y') / r³ + (x x' + y y') / r³ = 0

dus E' = 0, maakt niet uit welke beweging je hebt (daar heb ik geen uitspraak over gedaan).


Hetzelfde geldt voor L' :

L = x y' - y x'
L' = x' y' + x y' ' - y' x' - y x' ' = - (xy) / r³ + (xy) / r³ = 0


In jou geval moet je wel de massa meenemen omdat je die in je berekening ook meeneemt.. Maar dan zou het in principe constant moeten zijn.

[ Voor 63% gewijzigd door NickThissen op 19-05-2010 21:47 ]

[Mischa_NL]

Interessant, zelf ook binnenkort nodig...!
Via wiki wordt je niet veel wijzer en dit is toch wel even handig. (soort verkapte tvp)

[.oisyn]

NickThissen schreef op woensdag 19 mei 2010 @ 21:39:
Nee, de energie (niet alleen de kinetische energie, ook de potentiele energie erbij!)

Ja, daar zat ik net in de auto idd nog aan te denken. Kinetische energie kan afnemen, maar dat wordt "opgeslagen" als potentiele energie.

.edit: en je hebt idd gelijk: http://oisyn.nl/orbit.html laat idd zien dat de energie bij de reverse euler wat schommelt, maar bij RK4 constant blijft (althans, hij neemt heel langzaam af). Impulsmoment blijft bij de reverse euler overigens wel exact constant.

.edit2: ik heb 'm nu een tijdje lopen. RK4 blijkt toch nog wat instabiel - de energie neemt naar verloop van tijd af terwijl het impulsmoment toeneemt. De reverse euler implementatie blijft echter constant. Ja, de energie schommelt, maar wel per periode van rotatie, en hij blijft gewoon strak tussen de -2139 en -2170. Dus hoewel de afgelegde ellips wellicht niet exact klopt, blijft hij wel steeds datzelfde rondje draaien.

CRiMiNaL schreef op woensdag 19 mei 2010 @ 16:46:
Pak ik echter de test = false, dan stort de aarde al snel naar de zon, terwijl de waardes voor massa / velocity / positie allemaal lijken te kloppen.

Het viel me net trouwens op dat je G er een factor 1000 naast zit. Hij is niet in de orde van grootte van e-8, maar van e-11. Ik vermoed eerlijk gezegd dat dat dan ook de grote boosdoener zit (je waarde is te groot waardoor de aantrekkingskracht effectief groter is dan dat hij moet zijn, en de zaak dus al snel instort)

[ Voor 121% gewijzigd door .oisyn op 20-05-2010 02:19 ]

[NickThissen]

Bij mijn RK4 namen E en L ook langzaam af, met wat kleine fluctuaties (paar cijfers achter de komma, volgens mij meer dan jij laat zien), dus dat zal wel kloppen

[ Voor 12% gewijzigd door NickThissen op 20-05-2010 10:33 ]

[.oisyn]

Javascript gebruikt overigens wel doubles, geen floats.

[CRiMiNaL]

.oisyn schreef op donderdag 20 mei 2010 @ 00:14:
[...]

Het viel me net trouwens op dat je G er een factor 1000 naast zit. Hij is niet in de orde van grootte van e-8, maar van e-11. Ik vermoed eerlijk gezegd dat dat dan ook de grote boosdoener zit (je waarde is te groot waardoor de aantrekkingskracht effectief groter is dan dat hij moet zijn, en de zaak dus al snel instort)

Je hebt gelijk, overal staat G vermeld als 6.67428 X 10^-11. Begrijp niet dat ik dat verkeerd heb overgenomen. Nu blijkt een speed van 604800 (1 seconde simulatie = 1 week) geen probleem te zijn.

Nogmaals bedankt voor de uitleg en hulp.

[.oisyn]

Cool . Welke integratiemethode gebruik je nu?

[ Voor 13% gewijzigd door .oisyn op 20-05-2010 13:01 ]

[MSalters]

Ik zou moeten gaan graven, maar uit m'n hoofd is de (intuitieve) keuze voor een carthetisch coordinatenstelsel (x-y-z) niet zo heel erg handig in numerieke baanmechanica. Een polair coordinatenstelsel werkt een stuk handiger. Verder blijf je de methodiek van NickThissen gebruiken.

Het voordeel van een polair systeem is dat je relevante variabelen veel minder varieren. Zo zie je dat bij een perfect cirkelvormige baan dat dR/dt=0 en d²φ/dt²=0

[NickThissen]

Ik had inderdaad poolcoordinaten gebruikt, mede omdat je je stelsel differentiaalvergelijkingen dan om kan zetten in een enkele DV (ofzoiets, weet het niet meer uit m'n hoofd), en dat was de enige waar de (voorgeprogrameerde) RK-45 integrator iets mee kon. Met behulp van Langrange mechanica kun je heel makkelijk de DV afleiden in poolcoordinaten.

[G33rt]

NickThissen schreef op donderdag 20 mei 2010 @ 14:29:
Ik had inderdaad poolcoordinaten gebruikt, mede omdat je je stelsel differentiaalvergelijkingen dan om kan zetten in een enkele DV (ofzoiets, weet het niet meer uit m'n hoofd), en dat was de enige waar de (voorgeprogrameerde) RK-45 integrator iets mee kon. Met behulp van Langrange mechanica kun je heel makkelijk de DV afleiden in poolcoordinaten.

Klopt natuurlijk als een bus voor de cirkelbeweging; als je poolcoordinaten zou substitueren (vul maar in: x=r*cos theta(t), y=r*sin theta(t) ) is de enige interessante vergelijking die je overhoudt eentje voor tijdsafgeleiden van theta - en inderdaad is het net zo makkelijk de Euler-Lagrange vergelijkingen te gebruiken. De formulering op de Nederlandstalige Wikipedia is trouwens wel beroerd, voor de klassieke mechanica blijft het gewoon dit. In dit geval zijn de generaliseerde coordinaten r en theta - in het algemeen moet je die altijd interpreteren als "variabele die handig is om je probleem in uit te drukken", de rest is gewoon een fancy naampje.

Wat je eigenlijk aan het doen bent is een tweetal gekoppelde differentiaalvergelijkingen ontkoppelen door een handige substitie. (De definitie van "handig" hangt natuurlijk van je probleem af: in dit geval is een keuze voor poolcoordinaten overduidelijk de beste.)

In het geval dat je geen cirkel zou gebruiken maar - realistischer! - de zon in een van de twee brandpunten van een ellips zet krijg je natuurlijk wel tijdsafgeleides van r (omdat je nu x(t) = r(t) cos theta(t) gebruikt). Dat is overigens ook die energieverwarring van .oisyn - en bijvoorbeeld de reden dat een komeet (met zijn enorme eccentriciteit) de zon letterlijk als slingshot gebruikt: de straal wordt zo klein dat 'ie enorm snel de bocht om moet vliegen vanwege de perkenwet

Voor de geinteresseerden: de perkenwet eist dat het oppervlak per tijdseenheid een constante is. Met andere woorden:

de tijd tussen A-B en C-D moet hetzelfde zijn als de oppervlaktes even groot zijn.

[CRiMiNaL]

.oisyn schreef op donderdag 20 mei 2010 @ 13:00:
Cool . Welke integratiemethode gebruik je nu?

Ik gebruik nu de reverse euler implementatie, als ik alsnog tegen problemen aanloop bij het verhogen van dT zal ik een poging wagen de RK4 methodiek te gebruiken.
Zoals het voorbeeld op jou javascript al aangeeft, is de reverse euler goed genoeg voor dit projectje.