Three body problem: please help a willing-to-learn beginner.

[PerchloricAcid]

I've got this c++ homework. The idea is to code something, we got to pick our ideas by ourselves. After browsing the web, I got interested in the gravitational three body problem. The prof. says it's a good idea.

I need some help. I aren't even sure how to start!

After even more googling, I figured I would like to work with a somewhat simplified 3 body problem: a system that consists of one moving body and two fixed bodies on a 2D plane (something like this).

So, what am I supposed to do? Ok, I have to have some equations. So, I'll work assuming that energy and angular momentum are conserved, I'll have equations for velocity (dr/dt) and acceleration (d^2r/dt^2), and my objective is to solve for trajectories, r(t).
I'm not really sure what and how to do, though.

Just for the record, I'm not looking for an instant solution to my problem, or for someone to do my homework. I really want to get this done by myself, but I need some hints and directions. I'm almost a complete beginner and I believe that not learning to code earlier is one of my mistakes, but I'm totally willing to change that. I WANT to learn!!!

[Yakk]

Timeslice.

Determine where the bodies will go after a short period. Detect for collision. Check energy totals (and any other invariant) after the timeslice, and attempt to correct for it.

That last correction ends up being important. Repeated errors in each timeslice calculation is going to make even a stable system explode: by maintaining your known invariants, and re-injecting lost energy/angular momentum/momentum into the simulation, you might have some hope of things staying sane.

You might also want to pay attention to higher order terms, like jerk, to determine how big you can afford to make timeslices? Not sure. (Ie, like you'd use the higher order terms of a power series to bound the error of the first terms).

[PerchloricAcid]

Thanks for your reply!

By timeslice you mean this?
I haven't heard of this term yet.


Regarding jerk, it seems like a fine idea, I'll get back to that as soon as I code the basics of my simulation.

[Yakk]

No, I don't mean that.

I mean you have a time step size. You work out what the forces are at the start of the time period, then you pretend they are uniform over the period of time, and see where things are at the end. You then correct for the myriad of mistakes this results in -- maybe you recalculate the forces at the end, take the average, and repeat the timestep to see if that gives you a better value, maybe you calculate the energy and angular momentum of the system at the end and make sure it is the same as your invariant, maybe you look at higher order terms and decide to sub-divide the time step again.

There are probably people who do this professionally (or have trained to) on these forums. I'm just describing it from a general familiarity with the idea, rather than practical experience.

[PerchloricAcid]

Time step size is a familiar term. Thank you for the clarification.
I actually had that in mind.

[Sagekilla]

Time step size is a familiar term. Thank you for the clarification.
I actually had that in mind.

Since you're actually dealing with the simplified 3-body problem, you can actually formally expand the solutions
for r(t) and r'(t) in a power series that you can compute (theoretically) as many terms as you want.

Also, realize that r(t) is not a scalar value. It's vector valued. There's the x and y components.
You may be familiar with Newton's Law of Gravity as F = -G M1 M2 / r^2.

Try computing the vector form of this expression. Then, compute higher order derivatives. You'll realize that
because two bodies are fixed, something special happens with the derivatives.

When you want to actually solve the equations of motion for the lone particle, there's one of two routes to go.
You can choose to use well known integration schemes (Euler Method, Leapfrog Method, Runge-Kutta methods),
or you can choose to try the power series solution.

If you know how to calculate a Taylor Series, you should be able to compute r(t + dt) and v(t + dt) to solve your system.

If you want to go with the Taylor Series method, you should remember the definition of the derivative and you should
try to use that when you try to build your approximations. You'll notice that certain approximations are familiar if you've
read up on integration schemes.


Sorry for the wall of text, I did a lot of work on this in the past 4 years at my university. There's so much stuff to do in this

[PerchloricAcid]

Oh, thanks!

Which method do you think would be better for a beginner, the series or some of the integration methods?
I thought I'd work with RK4.

I tried to do some work tonight, and I must admit I got confused tons of times. C++ is what bothers me the most.
Btw, I'm a 2nd year physics student, and am therefore quite familiar with classical mechanics and introductory mathematical analysis, so no need to worry about whether I know about series etc.

[Zamfir]

To get it started, you might just use forward Euler at first. But try to write the program in a way that you could "plug in" another scheme relatively easily, by only replacing specific functions for example. That kind of separation is a good exercise anyway.

You can then setup the program first, create input, output and plotting options. Check how well invariants are preserved , play around with step sizes. If you get that working, replacing the scheme with Runge-Kutta can be the next step. If that still doesn't give you good results, look for example for a scheme that explicitly conserves energy.

That way you can first focus on the programming. And it's good to have personaly experience with how certain schemes fail. Otherwise you just implement something complicated that you don't understand, in order to prevent problems you didn't know existed.

[PerchloricAcid]

Also, realize that r(t) is not a scalar value. It's vector valued. There's the x and y components.
You may be familiar with Newton's Law of Gravity as F = -G M1 M2 / r^2.

Try computing the vector form of this expression. Then, compute higher order derivatives. You'll realize that
because two bodies are fixed, something special happens with the derivatives.

Higher order derivatives of what?
If you mean r, that is v=r'(t) and a=r''(t), they will be equal to zero as r doesn't change over time, but I believe that wasn't what you meant.

[PerchloricAcid]

Uh, I'm kinda stuck.

I have equations for r's, a's (x&z components) (from a=F/m), I have a function that says that f(t)=1, f(r)=v and f(v)=a (f is be d/dt), energy conservation equations. I don't know what to do next.

[Yakk]

Do you have a class or struct that represents "the state of the universe at a particular point in time"? If not, build that class or struct.

Produce a means to load that class or struct from, say, a file on disk. And maybe a way to save it. Or a way to set it from an interactive prompt. Ideally all 3, which lets you set up an initial state, save an intermediate or final state, and do a non-graphical run and look at the results later, etc.

Write "draw state to screen", by a method or function.

That'll give you the ability to see what is going on. Possibly you have already done these things.

Note that no evolution of state has occurred in any of the above stuff. Just load, save, data entry, and display.

---

Next, you'll want a crude way to introduce state evolution. KISS first. We'll ignore everything except forces, and have a fixed time step size.

Write "calculate forces on free body". This should be easy to do.
Write "apply force on free body for delta t time". This should be easy to do.

Now, write "crude time step", which takes a delta, calculates forces, applies them.

Next, "crude simulation with time step delta t for x iterations". It takes a state. It displays it. It runs crude time step(delta t) on it. It display it. It does this x times. Possibly it saves out the final result.

This will suck. But after this, you can improve it by keeping track of energy, angular momentum, etc. Or even build purely analytic solutions.

[Sagekilla]

Also, realize that r(t) is not a scalar value. It's vector valued. There's the x and y components.
You may be familiar with Newton's Law of Gravity as F = -G M1 M2 / r^2.

Try computing the vector form of this expression. Then, compute higher order derivatives. You'll realize that
because two bodies are fixed, something special happens with the derivatives.

Higher order derivatives of what?
If you mean r, that is v=r'(t) and a=r''(t), they will be equal to zero as r doesn't change over time, but I believe that wasn't what you meant.

Derivatives of position, beyond 2nd order. And also, the derivatives are not zero. The scalar form of Newton's Law of Gravity only
tells you the magnitude of gravity between two bodies. It tells you nothing of which direction it points.

The full equation should be written as F_ij(t) = G M_i M_j / (r_i - r_j)^2
Where F_ij is the magnitude of the force between bodies i and j. Try getting the vector form

Hints:

Spoiler:

F_ij is the magnitude of the force. What is the direction?

Spoiler:

The Force F points from body i to body j.

Spoiler:

What is the unit direction from body 1 to body 2?

Spoiler:

Given the magnitude and direction (separately) of a quantity, what is it's vector form?

[PerchloricAcid]

First of all, I would like to apologize for not writing for a few days (due to some personal issues). Thanks to everyone who contributed, your advice helped me quite a bit. I don't feel as lost as I felt when I posted my previous post anymore.

Yakk, could I keep the values that would represent "the state of the universe at a particular point in time" in a vector, without using classes or structs?
I did that before posting about being stuck. Of course, I could try to modify my code (and use classes or structs) if it would make further work much easier. I am highly inexperienced, so I'm unable to predict which would be better, although I am hoping that only vectors would be fine enough.

I wrote something that asks the user to enter the initial conditions, and then saves those values to a .txt file. Thanks for that direction. It seems like a pretty obvious thing, but didn't come to my mind back then when I had asked.
[I might add the option of taking the values from a file on my disk after I make something functional, but this user-interactive thing works for me for now.]
My intention is to make a program that would allow the user to choose initial conditions and see what happens when which initial conditions are chosen.

Yakk, are you basically suggesting that I should calculate (delta) momentum values and work with them? In other words, when you say "apply force on free body for delta t time", are you referring to dp=Fdt? I'm not quite sure whether I have understood you correctly. Sorry for having so much (elementary) questions.
[From this identity, I believe, I'd get the velocity as v=p/m.]


Sagekilla, are you talking about

Code: Select all

r12 = sqrt((x1 - x2)*(x1 - x2) + (z1 - z2)*(z1 - z2)),
r13 = sqrt((x1 - x3)*(x1 - x3) + (z1 - z3)*(z1 - z3)),
r23 = sqrt((x2 - x3)*(x2 - x3) + (z2 - z3)*(z2 - z3));

ax3 = - G*m1*(x3 - x1)/(r13*r13*r13) - G*m2*(x3 - x2)/(r23*r23*r23)
az3 = - G*m1*(z3 - z1)/(r13*r13*r13) - G*m2*(z3 - z2)/(r23*r23*r23);


something like this?
These would obviously be the second derivatives of position, that is, the x and z components of position; I'm working with a plane defined by the axes x and z.
The vector form would be a=ax*ex+az*ez, where ex and ez are unit vectors.


P. S.
Zamfir, although I hadn't replied to your post, I'd just like to say that I value that advice of yours and have decided to listen to it.
You're right, I would definitely learn coding better than if just applying RK4 or any other method as soon as I get to that part of work.

[Robert'); DROP TABLE *;]

Yakk, could I keep the values that would represent "the state of the universe at a particular point in time" in a vector, without using classes or structs?

Although I'm not Yakk, AFAIK, the answer to this is (like a lot of CS problems) "Yes, but it would be painful." This is a very object-orientated problem (since it involves actual objects in space) so it would make sense to use compound types like structs and classes to organise your data for you.

[PerchloricAcid]

[Yakk]

The reason you want to stuff it in a struct or class is so that you can add more state later on without having to rewrite your program.

If you want to keep track of invariant energy and angular momentum, for example, storing it in your state is the right spot.

[I might add the option of taking the values from a file on my disk after I make something functional, but this user-interactive thing works for me for now.]

Reading from disk is rather useful. It lets you halt a simulation and continue it later. Saving the state to disk is also useful, because it lets you compare the results of two simulations, among other things.

---

Yes, if your time slice is 0.1 seconds, and you have a force vector of v, then the change in momentum is going to be 0.1 * v.

I'm not sure what you mean by "vector" -- the physics, or the C++ std::vector?

The amount of code splatter there is rather bad.

Ie:

Code: Select all

struct vector_2 {
  double x, z;
};
double norm(vector_2 v) {
  return v.x * v.x + v.z * v*z;
}
double magnitude(vector_2 v) {
  return sqrt(norm(v));
}
vector_2 add( vector_2 left, vector_2 right )
{
  vector_2 result;
  result.x = left.x+right.x;
  result.z = left.z+right.z;
  return result;
}
vector_2 negate( vector_2 v )
{
  v.x = -v.x;
  v.z = -v.z;
  return v;
}
vector_2 subtract( vector_2 left, vector_2 right )
{
  return add( left, negate(right) );
}


Then:

Code: Select all

double r12 = norm( subtract(v1, v2) );
double r13 = norm( subtract(v1, v3) );
double r23 = norm( subtract(v2, v3) );


which makes what meaning of those lines much easier to read, and means that any bugs will be uniform over the 3 lines.

(Of course, operator overloading makes it even less verbose, but people get confused by that.)

We can then calculate gravity acting on the mobile body for one of them, and then on the other, then add them up. There is no need to do all of this on a single line of code -- that makes it harder to read, and hard to read code is more bug prone.

The rest of your math could be correct, but I am not all that interested in debugging it. It is definitely "too cute" -- optimize later, doing it before you have something working is just stupid.

And yes, even though acceleration due to gravity does not depend on the mass of the object in question, I would seriously calculate the force and then the acceleration. Make your code easy to read before you make it fast or efficient. Debugging a program will take you longer than writing it, so write a program to make it easy to be debugged as your primary goal. You can make it faster later, after it works.

[EvanED]

Debugging a program will take you longer than writing it, so write a program to make it easy to be debugged as your primary goal. You can make it faster later, after it works.

Kernighan had a great quote along the lines of "Debugging a program is twice as hard as writing it in the first place. If you write the code as cleverly as possible, you are, by definition, not smart enough to debug it."

[PerchloricAcid]

I'm not sure what you mean by "vector" -- the physics, or the C++ std::vector?

When I asked you about using vectors rather than classes or structs, I meant the C++ std::vector.

This alternative code you posted looks like something new to me.
I'll have to do a bit of googling to be able to fully comprehend it.

Thanks for the suggestions.

I'm not sure about when I'll be able to post back. As I've said, I've been going through some personal troubles lately, so now I must speed up my overall exam preparation, which means less time for coding. Just sayin' because I know that some people trying to help other people online really dislike the feeling that they've helped somebody, after which that particular somebody never replies, not even with a "thanks". I'll reply.

[PerchloricAcid]

Oh, golly. The thing is working. Working, I tell you

I feel so proud now that I've made my first simulation ever. It's still kinda crude, though, so there's still work to be done.

After I optimize my code, I might like to make an animation of the movement of the body.
However, I googled a bit, and it seems quite complicated to me.
So, uh, would that be hard/time consuming?

Thanks to everyone who helped!
Coding is fun

Yes, I'm euphoric!

[cemper93]

An animation? Do you ask about graphical output?
If yes, then the required work is heavily dependent on how complex you want it. If it is okay for you to simply display one dull square per object, without scale, rotation or user input - well, that should not be too difficult, even if you've never done it before. There are lots of easy-to-use, free graphics packages out there, where everybody's favorite seems to be the (not-so-much object oriented) SDL. Actually, the first example on the Wikibooks page might be all you need to know: Just render (draw) such a rectangle for every object you are simulating, measure the time you needed for that rendering (call SDL_GetTicks before and after rendering to measure the time) and use that as the time for your next simulation step, repeat. To limit framerate, you may want to make calls to SDL_Delay after rendering (before the second call to SDL_GetTicks), but that is not necessary.

[PerchloricAcid]

Yes, I'm asking about graphical output.
I would like to display my bodies as circles, preferably differently coloured, and I'd like to display the trajectory of the body that moves. I'd like to make it an animated trajectory, not just display the final state. Am I being clear enough?
The concept of my "program" is that the user is asked about the initial parameters, so I'd like to ask the user about the initial conditions, and then make an animated plot of the trajectory.

I'll take a look at the things you've posted about. Thanks

[EvanED]

One potential idea is the following. When I did a simulation like that (for a parallel programming class) what I did was have the actual simulation output a data file with the position of each object at each timestep, then had an entirely different program (I just used C# with some GDI calls like drawCircle or whatever) read that file and do the animation. Of course, I did this so I didn't have to do a GUI in C++ (which is what I wrote the simulation in), while you may be masochistic enough for that desire to not apply.

[PerchloricAcid]

Hm, if doing this in C++ would be masochistic, I would do as you did.
My simulation already outputs the data to a file, and from what you wrote, it seems to me that it wouldn't be so time-consuming to do it in C#. [Other exams have yet to be studied for, and I won't really have much time for other than studying.]
A colleague of mine [with much more experience than me] suggested I use C#, too.


Btw, I came up with an idea that seems pretty neat to me, but opinions from others are always welcome. While I was coding, I fucked around with different integration methods, and the code is written in such a way that the integration method is included from another file, so it would be easier to modify it. Even though I wasn't smart enough to save all I was doing, I figured I might make a comparison between the results that different integration methods produce (Euler-Cromer, RK4, Verlet (haven't tried this one yet, but willing to), Leapfrog...). That shouldn't even take up too much time, all I have to do is make the algorithm for integration and "call" it.
If this animation thing works out fine, it'd be easy to compare the integration methods.
The idea itself seems fun to me, and I'd definitely profit from practice. What do others think?

[Zamfir]

Yeah, that's the way to go. One thing about animations: for scientific problems, graphing is nearly always more useful than moving animations. If you have Matlab or something similar (or simply Excel) around, you should just use it read in your output files, and make plots of x, y, energy, momentum, rotational momentum vs time, and x vs y to see the orbits.

They often call this stage 'postprocessing': a system that makes the output of calculations quickly usable.

[PerchloricAcid]

One thing about animations: for scientific problems, graphing is nearly always more useful than moving animations. If you have Matlab or something similar (or simply Excel) around, you should just use it read in your output files, and make plots of x, y, energy, momentum, rotational momentum vs time, and x vs y to see the orbits.

I totally agree, and I'd rather do that. However, the thing is that this is homework for an exam, and it has to include *something* graphical, with buttons, forms, etc, so I figured that animations might be best.
Nevertheless, perhaps I'll just make a program that would show these plots! At first, I'd ask the user about initial conditions, then they would be given the choice click on a button to choose integration type, then click on a button to get plots of this or that, click on a button to see the orbits, etc. Perhaps it's not a bad idea at all. Thanks, dude

[Zamfir]

An animation doesn't have to be hard (not much harder than a workable plotting tool), many frameworks allow you to give commands like:


c= getCanvas()
c.clear()
c.drawCircle(x1, y1, bigRadius)
c.drawCricle(x2, y2, bigRadius)
c.drawCricle(xMoon, yMoon, moonRadius)

or something that sounds like those. The trouble starts afterwards. There's no reason why the stepsize of your animaton should be related to the stepsize of your calculation, so you need some interpolation. How do you rescale if the user puts in numbers that create an orbit far outside your standard window? They won't see anything otherwise.

Is there a fast-forward or replay? Do you want axes with numbers on them? After rescaling, won't those numbers be far too close to each other, or too far apart? Do you want to calculate first then play, or calculate on the fly as the user can fiddle around with the inputs? Do you want a line behind the moon that traces its orbit? If so, what happens if the moon crosses an existing trace line, will it cut out part of the line? Or do you remember the line and redraw it each frame?

If you like this programming stuff, then those are all very solvable problems and lots of fun, but each takes some time and it adds up fast. As a simple introduction to UI programming, you're better off making a box with a few buttons and input fields, that asks the user where to save the comma-separated output files. By scientific programming standards, that's already fancy. Quite some simulations work on the tried-and-true UI of putting some input variables at the top of the source code, then ask the user to recompile if they made changes. An advanced scientific UI means it can read variables from an input txt file.

If you really want to do animations, you might look at Processing. It's a animations and graphics language aimed at beginners, with lots of examples of simulations like yours. It also has a nearly perfect wrapper around javascript, so you can run your code in a browser at very little extra effort. Sadly, it doesn't do input forms very well.

[PerchloricAcid]

Oh my... there seem to be so many issues yet to be tackled! However, it seems kinda fun, too. I'm not that interested in UI programming actually, but as it's a must, I'll be glad to do my best and learn stuff while I'm at it.

Thanks for the detailed reply, you've given me great guidelines.

[sigsfried]

Hm, if doing this in C++ would be masochistic, I would do as you did.
My simulation already outputs the data to a file, and from what you wrote, it seems to me that it wouldn't be so time-consuming to do it in C#. [Other exams have yet to be studied for, and I won't really have much time for other than studying.]
A colleague of mine [with much more experience than me] suggested I use C#, too.


Btw, I came up with an idea that seems pretty neat to me, but opinions from others are always welcome. While I was coding, I fucked around with different integration methods, and the code is written in such a way that the integration method is included from another file, so it would be easier to modify it. Even though I wasn't smart enough to save all I was doing, I figured I might make a comparison between the results that different integration methods produce (Euler-Cromer, RK4, Verlet (haven't tried this one yet, but willing to), Leapfrog...). That shouldn't even take up too much time, all I have to do is make the algorithm for integration and "call" it.
If this animation thing works out fine, it'd be easy to compare the integration methods.
The idea itself seems fun to me, and I'd definitely profit from practice. What do others think?

This is certainly a good idea but be careful. If you are still solving a three body problem it is chaotic. This makes comparing integration schemes tricky as they will not converge.

[Yakk]

He has two fixed bodies. It isn't going to be all that bad.

[sigsfried]

It is still chaotic. How do you judge which of the solutions is correct, when they all have very different trajectories? Also notably you can't say something like the average trajectory is correct.

Also if you want something with a pretty result I highly recommend looking into the Duffing Oscillator.

[Zamfir]

Set one mass to zero, see how it looks.

[sigsfried]

That is fine, except some otherwise very poor solver, will solve perfectly that system.

I would recommend using the solvers on some other potential for example an anharmonic potential well.

[PerchloricAcid]

I completed the simulation, and I got max grade for it.
Now I'm willing to improve it and do at least some of the things I contemplated about in the thread.

Thanks for the feedback on comparing different integration methods. I might as well stick to the simplified three body problem, so as I see, it should all be fine? My code is written in such a manner that it will be rather easy to modify it so all of the bodies move.

[Sagekilla]

That is fine, except some otherwise very poor solver, will solve perfectly that system.

I would recommend using the solvers on some other potential for example an anharmonic potential well.

The simplest way is to analyze the two-body problem, or even the reduced one-body problem (central mass fixed).
Set it up with the initial conditions such that the orbiting body is in a circular orbit, and monitor the energy over time.
Measure the energy at t = 0 and every time thereafter. For symplectic integrators it should oscillate slightly but it will
be preserved to an extremely high degree. How poor the approximation is, is determined by the error w.r.t time and time step.

Alternatives: Subtract off the center of mass from each body, and subtract off the center of velocity from each body.
Basically, place it in the frame where the center of mass is at the origin and the center of momentum is at the origin.

Keep track of both of these quantities as you step the simulation in time. Both quantities will increase over time for
non-symplectic integrators, and you can determine how "good" the integration is by how fast these quantities grow
with respect to time and the time step.

[sigsfried]

I'm well aware of that, but there are schemes that can solve such systems perfectly, I was using one the other day I will look it up. They are mostly rubbish for more realistic systems.

[Sagekilla]

I'm well aware of that, but there are schemes that can solve such systems perfectly, I was using one the other day I will look it up. They are mostly rubbish for more realistic systems.

You may be but the OP certainly isn't -- that's why I posted that.

More info for the OP:

For systems which don't have periodic behavior, sometimes a high order dissipative integrator works better. You may be
familiar with the Runge-Kutta family of integrators. Practically all of them (I'm pretty sure every one of them is) are dissipative,
which means that when they're used on periodic orbits, the orbits will slowly decay over time. The symplectic integrator most
people are familiar with is the Leapfrog method, and most high order Runge-Kutta integrators (i.e. RK4, RKF, etc) will outperform
the Leapfrog method for non-periodic systems.

What I said about keeping track of the total energy (KE + PE for each body) and the center of mass / center of momentum still
holds true regardless of the system size. In the absence of any external forces (i.e. only gravity interacting), center of mass
and center of momentum should be preserved quantities in the center of momentum frame.

The better those quantities are preserved, the better your integrator is doing.

[PerchloricAcid]

Thanks
That will be very helpful to me

[sigsfried]

I may be wrong, but isn't the 3 body problem in a rotating reference frame? If so you have to be careful when dealing with conservation of energy.

You may be but the OP certainly isn't -- that's why I posted that.

I was only trying to point out to be cautious when checking for the quality of solvers, especially as such a generic code may well one day have far more complex solvers used.

[Sagekilla]

I may be wrong, but isn't the 3 body problem in a rotating reference frame? If so you have to be careful when dealing with conservation of energy.

You're probably right about. I believe it's a rotating frame. If you work in the center of momentum frame though, checking to make sure that center
of mass and center of momentum stay constant at the origin is still a valid way to make sure you're not getting a spurious solution.

It's quite easy for the OP to prove that using Newton's 2nd Law + 3rd Law in a system with only gravitational forces (what he's doing).

Keeping track of the magnitude of those two quantities is probably cheaper than calculating the energy of the system too.