//Goes through all the planets + sun. Gets total kinetic energy.
for(unsigned short i=0;i<iObject;i++){
dEnergy+=0.5*point.dMass*(pow(point.dVeloX,2)+pow(point.dVeloY,2)+pow(point.dVeloZ,2));
}
//Goes through all the objects and calculates the potential energy of the system.
for(unsigned short i=0;i<iObject;i++){
for(unsigned short j=0;j<iObject;j++){
if(i!=j){
dEnergy+=-(6.67e-11*point.dMass*point[j].dMass)/(Pytheg(point,point[j]));
}
}
}
double CDraw::Pytheg(gPOINT point,gPOINT point2){
return sqrt(pow(point2.dPosX-point.dPosX,2)+pow(point2.dPosY-point.dPosY,2)+pow(point2.dPosZ-point.dPosZ,2));
}
Total energy of a system cycles?
Author
Hello.
I was trying to check correctness of my solar system simulation and i decided to calcualte total energy and momentum. For momentum there is no problem, i get a nuber that never changes and is fine, but for my total energy i get a number that cycles forward and backwards......it is rather interesting since it is for the whole system.
A/w here is a small snippet of the code that i use to find total energy (potential and kinetic)
I reset the total energy each new calculations. This is for Euler;s method , i haven;t finish with the above simulation in rk4 yet. I was wondering if i was forgeting something when i was calculating the total energy that makes it cycles depending on the planets position.
well, I see that when calculating the potential energy you are taking all the objects into account, do you do the same when simuluating?
In other words, is it possible that when simulating the planets' movements only the sun's gravity was taken into account?
In other words, is it possible that when simulating the planets' movements only the sun's gravity was taken into account?
Author
nop, all forces are taken into the account and simulation is pretty precies, close to 1-2 days of being up to Saturn a little more later, but no more then 3 % error in time.
here si teh whole function which i proly posted like 10 times by know around here but it is easy to read :
Basicly momentum is good, never changes but the energy does change....so i think the error is in energy calcs and not the system.I was wondering if there is a way to find out if this is an erro due to Euleler's algo that i use or no.
here si teh whole function which i proly posted like 10 times by know around here but it is easy to read :
//Null the forces. for(unsigned short i=0;i<iObject;i++){ point.dForceX=0; point.dForceY=0; point.dForceZ=0; } //Calculate the new forces. for(unsigned short i=0;i<iObject;i++){ for(unsigned short j=0;j<iObject;j++){ if(i!=j){ point.dForceX+=ForceX(point,point[j]); point.dForceY+=ForceY(point,point[j]); point.dForceZ+=ForceZ(point,point[j]); } } } //Update the veoloctiy. for(unsigned short i=0;i<iObject;i++){ point.dVeloX=Velocity(point.dVeloX,Acceleration(point.dForceX,point.dMass),Time(1)); point.dVeloY=Velocity(point.dVeloY,Acceleration(point.dForceY,point.dMass),Time(1)); point.dVeloZ=Velocity(point.dVeloZ,Acceleration(point.dForceZ,point.dMass),Time(1)); } //Update the positions. for(unsigned short i=0;i<iObject;i++){ point.dPosX+=point.dVeloX*Time(1); point.dPosY+=point.dVeloY*Time(1); point.dPosZ+=point.dVeloZ*Time(1); } //Reset the checks. dMomentum=0; dMomentumX=0; dMomentumY=0; dMomentumZ=0; dEnergy=0; for(unsigned short i=0;i<iObject;i++){ dMomentumX+=point.dVeloX*point.dMass; dMomentumY+=point.dVeloY*point.dMass; dMomentumZ+=point.dVeloZ*point.dMass; dEnergy+=0.5*point.dMass*(pow(point.dVeloX,2)+pow(point.dVeloY,2)+pow(point.dVeloZ,2)); } //Calculate total momentum dMomentum=dMomentumX+dMomentumY+dMomentumZ; //Calculate total energy for(unsigned short i=0;i<iObject;i++){ for(unsigned short j=0;j<iObject;j++){ if(i!=j){ dEnergy+=-(6.67e-11*point.dMass*point[j].dMass)/Pytheg(point,point[j]); } } } //Do the time scaling using Time factordouble CDraw::Time(long double dTimeMs){ return (dTimeMs*TimeFactor);}//calculate the new velocity. newVelocity=oldVelocioty +AccelerationTimedouble CDraw::Velocity(double dLastVelocity,double dAcceleration,double dTime){ if(bFirst){ return dLastVelocity+dAcceleration*0.5*dTime; } return dLastVelocity+dAcceleration*dTime;}//Do acceleration from f=madouble CDraw::Acceleration(double dForce,double dMass){ return dForce/dMass;}//Calculate force in the X direction//Fx=forceGravity*changeinX/radiusdouble CDraw::ForceX(gPOINT point,gPOINT point2){ return (Gravity(point,point2)*(point2.dPosX-point.dPosX)/Pytheg(point,point2));} //Calculate force in the X direction//Fy=forceGravity*changeiny/radiusdouble CDraw::ForceY(gPOINT point,gPOINT point2){ return (Gravity(point,point2)*(point2.dPosY-point.dPosY)/Pytheg(point,point2));}double CDraw::ForceZ(gPOINT point,gPOINT point2){ return (Gravity(point,point2)*(point2.dPosZ-point.dPosZ)/Pytheg(point,point2));}//Gravity, Fg=Gmm/r^2double CDraw::Gravity(gPOINT point,gPOINT point2){ return ((6.67e-11*point.dMass*point2.dMass)/pow(Pytheg(point,point2),2));}//r^2=a^2+b^2//r=sqrt(a^2+b^2)double CDraw::Pytheg(gPOINT point,gPOINT point2){ return sqrt(pow(point2.dPosX-point.dPosX,2)+pow(point2.dPosY-point.dPosY,2)+pow(point2.dPosZ-point.dPosZ,2));} Basicly momentum is good, never changes but the energy does change....so i think the error is in energy calcs and not the system.I was wondering if there is a way to find out if this is an erro due to Euleler's algo that i use or no.
Quote:Original post by rever
I was wondering if there is a way to find out if this is an erro due to Euleler's algo that i use or no.
i dont know what kind of fluctuations you are talking about, but a bit of fluctuation is to be expected with euler integration, yes.
i believe the term for this is conservatism, and i believe verlet is conservative. its not that big a deal though, usually.
Author
well i am taalkign about fluctuation withing 3 percent up or down from the average...is somethign like that to be expected fro Euler's method?
cuz it is a bit of..by as much as 4% on the outmost planets.
cuz it is a bit of..by as much as 4% on the outmost planets.
I'm not sure, but I think when calculating the gravitational potential energy between two bodies, the pair should only be considered once. Your loop is considering the pairs (x,y) and (y,x) separately, so the potential energy is being counted twice. To consider each pair once, the loop should be something like for (i=1; i<
Don't you have to talk relativety into account for Mercurey since it's so close to the sun?
Quote:Original post by Anonymous Poster
I'm not sure, but I think when calculating the gravitational potential energy between two bodies, the pair should only be considered once. Your loop is considering the pairs (x,y) and (y,x) separately, so the potential energy is being counted twice. To consider each pair once, the loop should be something like for (i=1; i<
with regard to that:
in order to test the correctness of your energy calculations, to see if you miss a factor two somewhere or something, just place two stationary bodys, and let them attract. if that causes much bigger energy errors, then your energy calculation is wrong.
Author
I don;t think you should do the :
for (i=1; i<n; ++i) for (j=0; j<i; ++j)
but something like this should be used instead:
thx for pointing that out:D
Now that i use the above formula i still get some variation but it is much smaller.
When i did the calculations for 2 points i got the energy still changin but not so much as when i had 10 points.
[Edited by - rever on January 8, 2006 6:54:57 PM]
for (i=1; i<n; ++i) for (j=0; j<i; ++j)
but something like this should be used instead:
for(unsigned short i=0;i<iObject;i++){ for(unsigned short j=i+1;j<iObject;j++){ dEnergy+=-(6.67e-11*point.dMass*point[j].dMass)/Pytheg(point,point[j]); } } thx for pointing that out:D
Now that i use the above formula i still get some variation but it is much smaller.
When i did the calculations for 2 points i got the energy still changin but not so much as when i had 10 points.
[Edited by - rever on January 8, 2006 6:54:57 PM]
This topic is closed to new replies.
Advertisement
Popular Topics
Advertisement
