Gravitation attraction proportioned

[bkelly]

Begin with two stars in space, A is 1 solar mass in size, B is two solar masses. The distance between them is such that the equation for gravitational attraction works out to 3 meters per second squared. It seems to me that with A being smaller, it will accelerate towards B at the rate of 2 m/sec^2 while B will accelerate towards A at 1 m/sec^squared. Is this correct?

 

[Hootenanny]

Begin with two stars in space, A is 1 solar mass in size, B is two solar masses. The distance between them is such that the equation for gravitational attraction works out to 3 meters per second squared. It seems to me that with A being smaller, it will accelerate towards B at the rate of 2 m/sec^2 while B will accelerate towards A at 1 m/sec^squared. Is this correct?

No, that is not correct. Consider Newton's law of Gravitation,

$$F = G\frac{m_1m_2}{r^2}$$

And Newton's second law,

$$F = ma \Rightarrow a = \frac{F}{m}$$

Hence, if one calculate the acceleration a₁ of the first body with mass m₁ one finds,

$$a_1 = \frac{1}{m_1}\cdot G\frac{m_1m_2}{r^2} = \frac{Gm_2}{r^2}$$

And similarly for the second mass,

$$a_2 = \frac{1}{m_2}\cdot G\frac{m_1m_2}{r^2} = \frac{Gm_1}{r^2}$$

As you can see, the acceleration dye to gravity of a body does not depend on it's own mass, rather it depends on the mass of the body it is accelerating towards.

Do you follow?

 

[D H]

No, that is not correct. ... As you can see, the acceleration dye to gravity of a body does not depend on it's own mass, rather it depends on the mass of the body it is accelerating towards.

Hoot, that is exactly what the OP has. He is correct.

Begin with two stars in space, A is 1 solar mass in size, B is two solar masses. The distance between them is such that the equation for gravitational attraction works out to 3 meters per second squared. It seems to me that with A being smaller, it will accelerate towards B at the rate of 2 m/sec^2 while B will accelerate towards A at 1 m/sec^squared. Is this correct?

As Hoot mentioned, the gravitational accelerations of A toward B and B toward A as seen by an inertial observer are

$$\begin{aligned}<br /> a_A &= G\frac {m_B}{r^2} = G\frac {2m_A}{r^2}\\<br /> a_B &= G\frac {m_A}{r^2}<br /> \end{aligned}$$

Right off the bat one can see that that inertial acceleration of A is twice that of the inertial acceleration of B because B is twice as massive as is A. Since A accelerates toward B and B accelerates toward A, the acceleration of A relative to B (or B relative to A), given as 3 m/sec², is the sum of these accelerations:

$$\begin{aligned}<br /> a_{\text{rel}} &= G\frac {2m_A}{r^2} + G\frac {m_A}{r^2} = G\frac {3m_A}{r^2} \equiv 3\,\text{m}/{\text{sec}^2} \\<br /> \text{Thus} \\<br /> a_B &= G\frac {m_A}{r^2} = 1\,\text{m}/{\text{sec}^2} \\<br /> a_A &= 2 a_B = 2\,\text{m}/{\text{sec}^2}<br /> \end{aligned}$$

 

[Hootenanny]

Hoot, that is exactly what the OP has. He is correct.

 

[bkelly]

DH & Hoot,
Thanks for the confirmation. I do have a follow up. Start with:

$$F = G\frac{m_1m_2}{r^2}$$

I started out with that equation for force. Finding a value for the mass of the earth, I plugged in the values and simplified down to

$$F = \frac{9.803 meters}{Kg sec^2}$$

I now put this into words to see if my understanding is correct.

That given calculated force will cause one Kg of mass to accelerate at 9.803 meters per second squared. But part of saying that is general knowledge rather than a pure reading of the equation.

After rewriting this a few times and not getting what I want, I will try to simplify. The 9.8 meters I understand, as do I the per second squared. But the per Kg part seems awkward and is not making sense. If the mass is doubled to 2 Kg, then the accelration will half to 4.9 meters / second squared.

But just the same, the meters per kilogram part does not set well in my brain. I don't know what else to say, but am hoping you can say something that will help me alter my perspective.

Thanks for your replies.

 

[D H]

The force due to gravity near the Earth's surface is $F = mg$ where $g=9.8\,\text{m}/\text{s}^2$. It is not $F = 9.8\,\text{m}/\text{kg}/\text{s}^2$. The acceleration is (nearly) constant.

 

[bkelly]

I start out with
F = (G m1 m2 ) / r^2

If I expand G and drop the constants (makes it easier to track the types) I get:

G = meters^3 / Kg seconds^2

Adding in the m1m2/r^2 I get

F = (meters^3) ( Kg ) ( Kg ) / ( Kg ) ( sec^2 ) ( meters^2)

F = meters Kg / sec^2

I see one error in that Kg should have been in the numerator instead of the denominator. But my equation is still wrong. Where am I getting or not dropping the extra Kg?

DH said: The acceleration is (nearly) constant.
If an object falls towards Earth from 100 feet above the surface, the acceleration at the surface will be a tiny bit faster than at 100 feet. Was that your intent?

 

[D H]

F = meters Kg / sec^2

You can verify this via the SI version of Newton's second law, F=ma. The units of force in SI units is kg-m/sec².

In other words, your dimensional analysis is just fine.

(Note: The generic form of Newton's second law is F=kma. SI is constructed by design such that the constant of proportionality in Newton's second law is exactly one. The proportionality constant is not necessarily one in other systems. For example, in the English system with mass expressed in pounds-mass and force expressed in pounds-force the constant of proportionality is 1/32.1740486 lbf-sec²/lbm/ft.)

I see one error in that Kg should have been in the numerator instead of the denominator. But my equation is still wrong. Where am I getting or not dropping the extra Kg?

Your error was using little g, Earth standard gravitational acceleration, as a force. g is 9.803 m/sec². It is not just 9.803 (it has units), and it is not a force (it is an acceleration).

If an object falls towards Earth from 100 feet above the surface, the acceleration at the surface will be a tiny bit faster than at 100 feet. Was that your intent?

That is one of two reasons I said "(nearly) constant". Another reason is that larger masses do accelerate toward the surface of the Earth slightly faster than do smaller ones. Think of Newton's third law. Gravity makes the object accelerate toward the Earth and the Earth accelerate toward the object. The relative acceleration of an object with mass m toward the surface of the Earth is

$$G\left(\frac{M_e+m}{r_e^{\;2}}\right)$$

However, [itex]M_e\ggg m[/tex] for even the largest man-made object. The difference in acceleration between a one gram object and a one metric ton object is immeasurably small.[/itex]

 

[bkelly]

I am still missing something, but for the moment I would like to cut to the chase. I am writing a program that will provide the ability to enter the position, mass, and velocity of a few stars, then watch them zip around on the screen and see what happens. Going back to Hoot’s post it looks like I can use:


$$a_1 = = \frac{Gm_2}{r^2}$$

And

$$a_2 = \frac{Gm_1}{r^2}$$

This will calculate the acceleration imposed by any two stars on each other. (direction calculations omitted for now) I will do a round robin comparing every star in the simulator against every other star.

As I know the units, the calculation will be:

$$a_1 = \frac{6.67428 * 10^-11 * m2 }{r^2}$$

And if I have this right, I can express mass in tens (or hundreds) of million kilograms and radius in tens (or hundreds) of million meters and it should work out ok.

That's just the acceleration in a straight line between the two. There is more work to do. My plan is to do two passes through. First calculate and sum the acceleration vectors on each star, then on the second pass, calculate the new velocity and positions. If I can get that to work (looking more difficult all the time), then I may consider the effects of general relativity.

Thanks again for your time.

 

[D H]

That's just the acceleration in a straight line between the two. There is more work to do. My plan is to do two passes through. First calculate and sum the acceleration vectors on each star, then on the second pass, calculate the new velocity and positions. If I can get that to work (looking more difficult all the time), then I may consider the effects of general relativity.

Here is an even better way:
You will need a pair of functions (plus other functions to do your graphics) and a driver. One of the functions computes the acceleration vectors of all of the stars. Call this function the derivative function. The other function is your state propagator. This function takes as inputs a time step and state descriptions of each star (state meaning position, velocity, mass). The state propagator will call the derivative function and use this information to propagate the states of all of the stars to the end of the time step.

Here are two simple state propagators; both call the derivative function once per time step. The accelerations depend only upon position.In the math that follows, I have denoted the accelerations produced by your derivative function as f(x). (The bold indicates vectors).

Basic Euler method
$$\begin{aligned}<br /> \mathbf x(t+\Delta t) &= \mathbf x(t) + \Delta t\,\mathbf v(t) \\<br /> \mathbf v(t+\Delta t) &= \mathbf v(t) + \Delta t\,\mathbf f(\mathbf x(t))<br /> \end{aligned}$$

Euler-Cromer method
$$\begin{aligned}<br /> \mathbf v(t+\Delta t) &= \mathbf v(t) + \Delta t\,\mathbf f(\mathbf x(t)) \\<br /> \mathbf x(t+\Delta t) &= \mathbf x(t) + \Delta t\,\mathbf v(t+\Delta t)<br /> \end{aligned}$$

These two methods differ only in the order in which the positions and velocities are updated. In the basic Euler method, position and velocity are updated simultaneously. In the Euler-Cromer method, the position is updated using the updated velocity. If you want to use one of these simplistic techniques, use Euler-Cromer. Energy is more or less in Euler-Cromer conserved but not with the basic Euler method.

The advantage of writing your simulator in the fashion I outlined is that doing so enables you to use more advanced propagation algorithms without changing any other parts of your code. How you propagate state should have nothing to do with how you initialize the system, compute acceleration, or display things graphically.

A slightly more advanced integrator is the velocity verlet:
$$\begin{aligned}<br /> \Delta \mathbf v_1(t) &= \Delta t\,\mathbf f(\mathbf x(t)) \\<br /> \mathbf x(t+\Delta t) &= \mathbf x(t) + \Delta t(\mathbf v(t) +<br /> \frac 1\ 2 \mathbf \Delta v_1(t)) \\<br /> \Delta \mathbf v_2(t) &= \Delta t\,\mathbf f(\mathbf x(t+\Delta t)) \\<br /> \mathbf v(t+\Delta t) &= \mathbf v(t) +<br /> \frac 1 2(\Delta \mathbf v_1(t) + \Delta \mathbf v_2(t))<br /> \end{aligned}$$

Here you call the derivative function twice per time step. This is *a lot* (not just twice) more accurate than either of the Euler methods.

Even more accurate techniques include Runge-Kutta methods (RK4 being the most famous), various Adams-Bashforth-Moulton methods (ABM3 being the most widely used of these methods), and many, many other techniques.

 

[bkelly]

Hello DH,
I did not want to go into much detail on program concepts, but we are clearly in agreement. I am creating a class for a star object. The class will contain the mass, position and velocity. It will be able to query other star objects for their mass and position, then calculate the imposed acceleration. On the first pass, every star will query every other star to determine a composite acceleration rate. On the second pass, the acceleration will be applied to velocity then to position over the time interval. Then I start over. As you noted, all the processes of initialization, display, calculating accelerations, velocities, and positions will be as independent as possible.

To your last set of equations. I would like to verify that I understand what you have taken the time to write. I understand that these equations are applied after all the composite accelerations have been calculated. The first of the last four:

$$\Delta \mathbf v_1(t) &= \Delta t\,\mathbf f(\mathbf x(t)) \\$$

This says that the change in velocity and a given time is equal to the change in time multiplied by the acceleration rate. The f(x(t)) referring to the earlier discussions in this thread for acceleration.

Pass 1, part 1: Calculate the acceleration rate (vector) generated by the relationship between each pair of stars.
Part 2: Determine the angle of the vector
Part 3: Break the acceleration down into x, y, and z components
Part 4: Add the components to the current velocity. (delay until pass 2?)

This does imply that I intend to keep the acceleration and velocity vectors in x, y, and z components rather than a polar type component.

I think the next three equations are applied in part 2. I will state how I interpret the first one and see where I stand.

$$\mathbf x(t+\Delta t) &= \mathbf x(t) + \Delta t(\mathbf v(t) +<br /> \frac 1\ 2 \mathbf \Delta v_1(t)) \\$$

The new position of the star is equal to:
The current position + the change in time multiplied by the current velocity + negative 2 times the change in velocity at time t.

I don’t understand the superscripted 1 in front of the last delta symbol, but I suspect that you really wanted the fraction ½ delta v1 at time t. That would make that phrase: plus one half delta v1 at t.

As I read this I get the idea that the star object might need to keep the current velocity and acceleration separate from the delta velocities and acceleration until the first step is complete, then I will have both for the second step when applying the acceleration to the current velocity and position.

I don’t know what the next two equations do yet or how they work together, but I would like to go one stage at a time. When I understand them individually, then, hopefully, I can put them together in my head.

 

[D H]

Hi, kbelly.

On the first pass, every star will query every other star to determine a composite acceleration rate. On the second pass, the acceleration will be applied to velocity then to position over the time interval. Then I start over.

Don't do it this way. This architecture will restrict you to the simplest of integrators. It is much better if the integrator calls the derivative function. This will allow you to advance to integration techniques that are much more accurate the the primitive Euler techniques.

To your last set of equations. I would like to verify that I understand what you have taken the time to write.

Oops. I wrote too quickly writing and spent even less time reviewing what I wrote.

Corrected velocity verlet equations:
$$\begin{aligned}<br /> \Delta \mathbf v_1(t) &= \Delta t\,\mathbf f(\mathbf x(t)) \\<br /> \mathbf x(t+\Delta t) &= \mathbf x(t) + \Delta t(\mathbf v(t) +<br /> \frac 1\ 2 \mathbf \Delta v_1(t)) \\<br /> \Delta \mathbf v_2(t+\Delta t) &= \Delta t\,\mathbf f(\mathbf x(t+\Delta t)) \\<br /> \mathbf v(t+\Delta t) &= \mathbf v(t) +<br /> \frac 1 2 (\Delta \mathbf v_1(t) + \Delta \mathbf v_2(t+\Delta t))<br /> \end{aligned}$$

In pseudocode, assuming Star contains public members position, velocity, acceleration,

VelocityVerlet::integrate (
   double dt,
   Star * stars,
   int nstars)
{
   compute_accelerations (stars, nstars);
   for (int ii = 0; ii < nstars; ii++) {
      delta_v1[ii] = dt * stars[ii].acceleration;
      stars[ii].position += deltat *(stars[ii].velocity + 0.5*delta_v1[ii]);
   }
   compute_accelerations (stars, nstars);
   for (int ii = 0; ii < nstars; ii++) {
      delta_v2[ii] = dt * stars[ii].acceleration;
      stars[ii].velocity += 0.5*deltat *(delta_v1[ii] + delta_v2[ii]);
   }
}

The velocity verlet integration function evaluates the derivative at the start and end of the integration interval. Note that you can eliminate the first call to the derivative function by calling the derivative function at startup time. (This optimization will not be valid if you add in a general relativistic correction, as the acceleration becomes a function of position and velocity.)

BTW, here is RK4 (here written with the acceleration function taking position and velocity as inputs):
$$\begin{aligned}<br /> \mathbf x_0 &= \mathbf x(t) \\<br /> \mathbf v_0 &= \mathbf v(t) \\<br /> \mathbf a_0 &= \mathbf f(\mathbf x_0, \mathbf v_0) \\<br /> \mathbf x_1 &= \mathbf x_0 + \frac{\Delta t} 2 \mathbf v_0 \\<br /> \mathbf v_1 &= \mathbf v_0 + \frac{\Delta t} 2 \mathbf a_0 \\<br /> \mathbf a_1 &= \mathbf f(\mathbf x_1, \mathbf v_1) \\<br /> \mathbf x_2 &= \mathbf x_0 + \frac{\Delta t} 2 \mathbf v_1 \\<br /> \mathbf v_2 &= \mathbf v_0 + \frac{\Delta t} 2 \mathbf a_1 \\<br /> \mathbf a_2 &= \mathbf f(\mathbf x_2, \mathbf v_2) \\<br /> \mathbf x_3 &= \mathbf x_0 + \frac{\Delta t} 2 \mathbf v_2 \\<br /> \mathbf v_3 &= \mathbf v_0 + \frac{\Delta t} 2 \mathbf a_2 \\<br /> \mathbf a_3 &= \mathbf f(\mathbf x_3, \mathbf v_3)<br /> \end{aligned}$$

The updated position and velocity are a weighted average of these intermediate states:

$$\begin{aligned}<br /> \mathbf x(t) &= \mathbf x_0 +<br /> \frac{\Delta t} 6 (\mathbf v_0 + 2\mathbf v_1 + 2\mathbf v_2 + \mathbf v_3) \\<br /> \mathbf x(t) &= \mathbf v_0 +<br /> \frac{\Delta t} 6 (\mathbf a_0 + 2\mathbf a_1 + 2\mathbf a_2 + \mathbf a_3)<br /> \end{aligned}$$

 

[bkelly]

Originally Posted by bkelly
On the first pass, every star will query every other star to determine a composite acceleration rate. On the second pass, the acceleration will be applied to velocity then to position over the time interval. Then I start over.

DH wrote: Don't do it this way. This architecture will restrict you to the simplest of integrators. It is much better if the integrator calls the derivative function. This will allow you to advance to integration techniques that are much more accurate the the primitive Euler techniques.

As we carry on this discussion I am building the foundation of this program. I have created the basic galactic body class and am working on dialogs to view and enter data. I have not begun the work on calculating the values so I am completely flexible in that regard.

As I envision this work, I see a list of galactic bodies (i.e. stars). In order to calculate the next step (new velocities and positions), each star must be compared with all the remaining starts to determine the cumulative gravitational effect. It seems to me that I must calculate these gravitational effects for all the stars before I update the position of any star. Otherwise, a star would move in the middle of the calculation and all calculations following that move would be a half breed of two time intervals.

That is why I see two passes. The first pass determines the force exerted on each star at the beginning of a time interval. The second updates the velocities and positions of each star to what is expected at the end of the interval, which is the starting point of the next interval.

I can handle the programming part of this, (and do very much appreciate any coding suggestions you provide) but it is very clear that you (and most the people in this forum) know much more about the gravitational and motion equations than I. I am a software/systems engineer, but I had to struggle mightily to pass Calculus Two.

When you say “It is much better if the integrator calls the derivative function.” I lose track of which calculation you are referring to. Can you back up and state that in more of a layman’s terms.

 

[D H]

It seems to me that I must calculate these gravitational effects for all the stars before I update the position of any star. Otherwise, a star would move in the middle of the calculation and all calculations following that move would be a half breed of two time intervals.

That is correct.

That is why I see two passes. The first pass determines the force exerted on each star at the beginning of a time interval. The second updates the velocities and positions of each star to what is expected at the end of the interval, which is the starting point of the next interval.

This is where I am having an issue. Developing your simulation in this way will inherently constrain you to using a propagator with low accuracy.

When you say “It is much better if the integrator calls the derivative function.” I lose track of which calculation you are referring to. Can you back up and state that in more of a layman’s terms.

The "derivative function" to which I referred is that function which calculates the accelerations of each star due to gravitational attraction to all of the other stars.

Here is a suggested architecture:

class Star {
   public:
   double mass;
   double position[3];
   double velocity[3];
   double acceleration[3];
};

class Stars {
   Star * stars; // Or some other mechanism (e.g. vector) for storing star data
   int nstars;

   // Initializer function. This function allocates data for the stars and sets each star's
   // mass, initial position, and initial velocity.
   void compute_accelerations();

   // The derivative function. This function computes accelerations for all of the stars.
   void compute_accelerations();
};

class Propagator {
   void initialize (Stars & stars);
   void update (Stars & stars);
};

void simulation_driver (void)
{
   Stars stars;
   Propagator propagator;
   Graphics graphics;
   double deltat;
   double time;

   // Initialize the system.
   time = 0.0;
   deltat = ...;
   stars.initialize ();
   propagator.initialize (stars);
   graphics.initialize (stars);

   // Main loop. This loop propagates the stars to the next time step and portrays the
   // updated star positions graphically.
   while (! finished) {
      propagator.update (stars, deltat);
      time += deltat;
      graphics.update (stars);
   }

   // Cleanup code. You should play nice and free memory and other resources.
}

Some things to note:

• All of the accelerations must be computed before any stellar position is updated.
• The Stars::compute_accelerations method does just this.
• The Propagator::update method is responsible for computing the accelerations (via stars.compute_accelerations()), not the driver.
• Whether the propagator calls stars.compute_accelerations() once per time step or more is up to the propagator, not the driver. Better propagators will compute the accelerations multiple times per time step (i.e., one pass through the main loop).

 

[bkelly]

This is where I am having an issue. Developing your simulation in this way will inherently constrain you to using a propagator with low accuracy.

I don't understand why you say that. My comment further down on classes may clarify.

The "derivative function" to which I referred is that function which calculates the accelerations of each star due to gravitational attraction to all of the other stars.

Then I take it that the integral part of this operation is to apply the accelerations to produce new velocities and positions. Correct?

Regarding classes:
The class encapsulates the parameters of an object, and the methods that operate on those parameters. There is almost never a valid reason for making variables public.

The glactic object class will contains all the mass, position, and velocity variables. One object will be created for each star. All share the executable code, but do not share any of their variables. The base object contains the code to perform updates. To determine how star 1 is affected by star 2, I will call a procedure in star1 and pass it a pointer to star 2. The star 1 object will then query star 2 for its mass and position. Star 2 will provide a copy of the values it contains for mass and position. The star 1 object will use that data provided by star 2 to to determine the accelerations induced by star 2. The variables contained within any star object are never exposed to access from outside the object.

A dialog box will provide the means to update values such as initial mass, position, and velocity. That dialog box will call methods of the star object and give those methods the values to update. The class code will make the update with a call typified by:

Star_Pointer->Set_X_Pos( Value );

The dialog box never directly accesses the internal variables. In this manner, the procedure Set_X_Pos() can examine the values passed in, determine that they are within specified limits, and apply them to the internal variables.

I browsed around a bit and was not able to find a FAQ or other directions concerning formatting text in posts. Some controls are fairly common such as puting "quote" within brackets, but others are not. Where do I find that information about how to create a code box with and without scroll bars.

 

[D H]

Regarding formatting of posts: It always helps to read the FAQ, https://www.physicsforums.com/faq.php" . For example [noparse]

code ...

[/noparse] will display the "code ..." as code (monospaced font, etc). The display generator determines whether scrollbars are needed.

Each thread has a new "New Reply" button displayed near the bottom of the page. You have been you clicking that button. This brings up a form in which you type your response. A set of vBcode shortcut buttons lies above the entry panel. One of these (the "#" button) wraps CODE tags around highlighted text (or creates an empty set of CODE tags if no text was highlighted.)

One other thing of note: vBcode as-is doesn't display math very well. This is not good for a physics and math forum. This forum displays math written in LaTeX in the form [noparse]$$LaTeX code ...$$[/noparse].

 

[D H]

Regarding classes:
The class encapsulates the parameters of an object, and the methods that operate on those parameters. There is almost never a valid reason for making variables public.

I know about object oriented programming. Unfortunately, it like many other programming paradigms, has a religious aspect to it. The debate over public/private is, to some extent, just religion. We don't condone religious debates at this forum.

A dialog box will provide the means to update values such as initial mass, position, and velocity. That dialog box will call methods of the star object and give those methods the values to update.

Consider reading the star specifications from a file. The dialog box mechanism is not very good if you want to have more than a handful of stars in the system or if you want to explore variations on a theme (e.g., Monte-Carlo analyses).

I browsed around a bit and was not able to find a FAQ or other directions concerning formatting text in posts. Some controls are fairly common such as puting "quote" within brackets, but others are not. Where do I find that information about how to create a code box with and without scroll bars.

I mentioned this in my previous post. See the hyperlinks I posted there. Regarding code boxes, the user doesn't have control over the presence of scroll bars. The user specifies the content of a code box. The display software determines whether scroll bars need to be used. It adds horizontal scroll bars if the basic code box is too wide, vertical scroll bars if the basic code box is too tall.

 

[bkelly]

DH wrote:

This is where I am having an issue. Developing your simulation in this way will inherently constrain you to using a propagator with low accuracy.

I remain concerned about this post as I don’t understand why you say that. I would really like to back up to that point in our conversation and hear what you have to say on this.

DH wrote:

I know about object oriented programming. Unfortunately, it like many other programming paradigms, has a religious aspect to it. The debate over public/private is, to some extent, just religion. We don't condone religious debates at this forum.

Writing from the perspective of a professional engineer with significant experience in writing software that must perform properly and must be maintainable by others, I impose stricter guidelines on the code that I and my peers write than most people do. To me, the guidelines do not take on a religious aspect as they are backed up by evidence. However, as you say, and noting that this is not a programming forum, its best to not go there. I mentioned those items only to provide the concept of the environment in which my class will operate.

DH wrote:

Consider reading the star specifications from a file. The dialog box mechanism is not very good if you want to have more than a handful of stars in the system or if you want to explore variations on a theme (e.g., Monte-Carlo analyses).

An excellent point. Before I go there I want get the simulator running with two, then three or four stars. I will provide some default values for the first few stars. After I get the fundamentals working, then it will be time to add the ability to save simulation parameters to a file and restore them. OTH, even when entering two or three stars, it may well prove tedious and error prone to make the entries by hand. I just may implent the save and restore functions earlier than planned.

 

[D H]

I remain concerned about this post as I don’t understand why you say that. I would really like to back up to that point in our conversation and hear what you have to say on this.

So, let's back it up.

My problem is your wanting to do things in two stages. I took you to mean doing something like this in your driver:

while (! finished) {
   galaxy.compute_acceleration();
   galaxy.update_positions_and_velocities(deltat);
   time += deltat;
   graphics.update(galaxy);
}

That is fine for the simplest of propagators, i.e., the Euler methods. The general name for problems of this sort is "initial value problems". An initial value problem comprises the initial state of some set of parameters (in this case, the initial positions and velocities of a bunch of stars) and a set of ordinary differential equations that governs the time evolution of the state of the system (in this case, $\frac{d}{dt}(\mathbf x(t)) = \mathbf v(t)$ and $\frac{d}{dt}(\mathbf v(t)) = \mathbf a_{\text{grav}}(\mathbf x(t))$. Mathematicians, scientists, and engineers have investigated solving these initial value problems for the past two or three hundred of years or so. Doing things "in two stages" describes the workings of the simplest of techniques. Better techniques do things in multiple stages. If you establish an architecture up-front that accommodates a wide range of techniques, tada, you will be able to use a wide range of techniques to solve the problem. If, on the other hand, you establish an architecture up-front that accommodates only the simplest of techniques, you might will be hard-pressed to move on to better techniques.

 

[bkelly]

Replying to DHs last post, et al:

Understood. I did not present a clear picture of my intentions. The basic concept of implementation will have an effect on the outcome so here is a high level pseudo code of the control.

Please note that my spelling is atrocious. This was edited in Microsoft Word which tends to enforce capitalization rules capriciously. Please ignore any errors of that nature.

There will be a class called Galactic_Body. It will contain all the variables and code necessary to track the position of a start in simulated space, and calculate and remember the effects of gravity induced by other star. N number of instances of this class will be connected via a linked list for easy traversal.

After considerable thought I have decided that further details would be inappropriate for this forum. This may already push the limits, but I think DH (and possibly any other reader) will be more comfortable with this conversation if I post this.

This is an outline of the code that conducts each time interval of the galactic simulator.

Do
{
      // stage 1
      // get the first body in the list
   Star_pointer = link_root->next;   
   Do
   {
   Star_pointer->Traverse_List_For_Gravity();
   Star_pointer = star_pointer->get_next_body();
   }(until end of list )

       // stage 2
      // get the first body in the list
   Star_pointer = link_root->next;   
   Do
   {
   Star_pointer->Integrate_Acceleration_into_Next_Position();
   Star_pointer = Star_pointer->get_next_body();
   }(until end of list )

// insert code for housekeeping, update time, iteration counts, etc.
// insert code to check for operator input, maybe several places

}(until done);

Within the class is a method to traverse the list of stars, check each star’s position and mass, and calculate the effect of gravity. It will go something like this.

Void Traverse_List_For_Gravity()
{
Other_pointer = link_root->next;
Do
   {
      If( Other_pointer == me) continue;
         // structure Other_body is a local structure for temporary storage of the values
         // retrieved from the other body.
      Other_body.mass = Other_pointer->query_mass();
      Other_body.position = Other_pointer->query_position();
      Calculate_gravitational_effect();
      Other_pointer = Other_pointer->next;
   }(until end of list);
} // end of function

Now to the function that seems to be of primary importance:

   Star_pointer->Integrate_Acceleration_into_Next_Position();

This is the essence of what I am referring to as stage two of a single time interval. This function, contained within the Galactic_Body class, can be as simple or as complex as desired. Within this function, the code can iterate through the data as many times as necessary. If the time interval is one hour, it could calculate intermediate values for each second. I can start with an overly simple function, then change it as desired. Other than consuming time, changes to this function will have no effect on the remainder of the program.

After the first stage, each body will be aware of the gravitational effect from each other body, but will not have implemented any positional changes. Then the second stage begins and does the positional updates. (I think velocity can be calculated in stage 1, but am not certain.) The concepts to be implemented in method Integrate_Acceleration_into_Next_Position(); is where I need help.

Now the proximate question becomes, does this resolve the questions that the reader may have with my method? Or does it show a flaw in how I plan on conducting the simulation?

In the meantime, I am new to Visual Studio and having difficulties in getting the framework operational. I will be working that and looking at the equations already presented.

I do wish to say to DH, Your comments have been extraordinarily helpful. I am aware of and appreciative of the time you have consumed in answering my questions.
Thank you.

 

[D H]

This is an outline of the code that conducts each time interval of the galactic simulator.

Do
{
      // stage 1 [ compute accelerations]
       // stage 2 [integrate]
      // get the first body in the list
}(until done);

This is exactly what I suggested you not to do. Dang. Frustrated.

OK. You are well-versed in the concepts of information hiding. Think of the computation of accelerations (your Traverse_List_For_Gravity() method) as something that should only be called by the Integrate_Acceleration_into_Next_Position() method. Whether this function needs to compute the derivatives once, twice, or more times per simulation time step is a function of the integration algorithm. Doing this "stage 1" / "stage 2" thing in the driver is bad; very, very bad. Just to give you an inking of how bad an idea I think this is, this architecture is just as bad as making all of your member data public. A far better approach is to make your integrator call Traverse_List_For_Gravity() rather than your driver.

There are a plethora of algorithms for solving this class of initial value problems. From an OO perspective, you can envision these different algorithms as different child classes of some basic InitialValueProblemIntegrator class. View this InitialValueProblemIntegrator as a pure virtual class with a pure virtual method Integrate_Acceleration_into_Next_Position() method. Start with one child class of this virtual class that implements the integrate method. You can then add different child classes that implement different integration algorithms. This will give you the ability to compare different integration algorithms. The stage 1 / stage 2 approach inherently limits you to what are widely viewed as very weak integrators.

 

[bkelly]

DH wrote:

This is exactly what I suggested you not to do. Dang. Frustrated.

I am hoping you don't get frustrated. Check the last statement of my previous post. Be aware that I am appreciative of your responses. You have me sitting on the edge of my seat. Know also that I have not implemented those phases. I am not ready for that implementation and if I was, I would hold off until I understand your position.

As I am not understanding your concept, here is mine as briefly as I can.

Assume N number of stars in a galaxy at time t1. Each has a gravitational on every other. When star 1 changes position, this change modifies the effect on all other stars. The reason I want two passes is to first, calculate the effect that each causes on the other, before changing the position of any. I can save those calculations somewhere and not change any positions until I have determined all the effects. This is understanding the affects at a given instant in time.

Only after I know those effects, can I begin moving the stars to t2. Once I move any star from its t1 position to a t2 position, then the overall environment at t1 has been disturbed. If I move any star before I have determined how it affects all other stars, then simulation accuracy will be degraded.

Possible thoughts:
Do you think I should evaluate each star, and during that single evaluation, calculate the effect of gravity by every other star, and create a new velocity and position the first time I visit that star. To protect the integrity of the t1 environment, I could save the new values in a ready location, then activate them at the end of t1 or the beginning of t2.

And I do have to ask: In the process of calculating the effect of one star on another, and the results of those effects, there are several stages of processing to complete. As computer simulations require finite amounts of time to do their calculations, every calculation is displaced from every other calculation in real time by some degree. [simulated time is frozen] If the gravity/acceleration/velocity/position update for star 1 is interrupted while something else goes on, that has no effect on star 1, how will the calculation for star 1 be spoiled?

This gets to the concept of a computer simulation. We effectively freeze simulated time for some amount of real time while we perform all the calculations (all that we know about, or want to implement) needed to determine the interactions that take place (effectively) instantly in the analog/real world. Sometimes the simulations are deliberately slower than real time, sometimes the same, and sometimes faster. It all depends on the project.

In summary, I just don’t see how separating the calculation in real time, as opposed to simulated time, generates any inherent limitations. OTH, I am certain there are things of which you are aware that I am not.

DH wrote:

OK. You are well-versed in the concepts of information hiding.

Thank you.

But there is a bit more to it than hiding. It is also pushing the details of a complex calculation down to the point where the only thing to do in a function is that one complicated process. When you know that the function does that one thing and that one thing only, and you know to throw out any code that does not contribute to that one single clear task, then you have a good function. That isolation reduces what is called coupling which in turn, reduces the unanticipated side effects that changes to it will generate.

That’s it for today, and maybe tomorrow. I will be looking at your equations from several posts back and trying to understand them one step at a time. Again, I do appreciate you taking the time to write to me.

 

[D H]

Assume N number of stars in a galaxy at time t1. Each has a gravitational [noparse][attraction][/noparse] on every other. When star 1 changes position, this change modifies the effect on all other stars. The reason I want two passes is to first, calculate the effect that each causes on the other, before changing the position of any. I can save those calculations somewhere and not change any positions until I have determined all the effects. This is understanding the affects at a given instant in time.

You are correct in that you need to calculate all of the accelerations before you update any of positions. The problem is that some propagators need to perform a bunch of mini-steps for each simulation step. For example, fourth-order Runge-Kutta method computes the derivatives four times per simulation step. The updated state is based on a weighted average of these four derivatives.

Possible thoughts:
Do you think I should evaluate each star, and during that single evaluation, calculate the effect of gravity by every other star, and create a new velocity and position the first time I visit that star. To protect the integrity of the t1 environment, I could save the new values in a ready location, then activate them at the end of t1 or the beginning of t2.

I think you should separate state propagation from everything else. I also think you should read up on techniques for solving initial value problems.

But there is a bit more to it than hiding. It is also pushing the details of a complex calculation down to the point where the only thing to do in a function is that one complicated process. When you know that the function does that one thing and that one thing only, and you know to throw out any code that does not contribute to that one single clear task, then you have a good function. That isolation reduces what is called coupling which in turn, reduces the unanticipated side effects that changes to it will generate.

The propagator does do exactly one thing: It takes the state of the system from one time step to the next. It obviously needs to compute derivatives, and these need to be computed from a consistent system state (and a good way to do this is by having a function that atomically computes all of the derivatives). However exactly how the propagator calls the derivative function and exactly what it does with the calculated derivatives is best hidden. By calling the derivative function outside of the propagator you are poking at the innards of the propagation function.

That said, here is an alternate approach to state propagation:

do {
   do {
      system.compute_derivatives ();
      propagation_finished = system.propagate_state (deltat);
   } until (propagation_finished);
   time += deltat;
   system.do_something_with_new_state ();
} until (simulation_finished);

The simplest of propagators always returns true (i.e., state has been propagated). More complex propagators will need multiple executions of the inner loop before indicating that the state has been propagated to the next simulation time step.