MATLAB Where Should I Include Particle Mass in My Hamiltonian CR3BP Code?

[Deadstar]

Hey folks I'm trying to code the Hamiltonian form of the circular restricted three body problem and would like a bit of help here. I've managed to code the non-Hamiltonian form and have gotten it to work fine but came run into trouble here.

I'm using the synodic coordinate system as described here http://en.wikipedia.org/wiki/Jacobi_integral#Synodic_system and have found the Hamiltonian form of the system to be...

H = \frac{1}{2}(p_x^2 + p_y^2) - xp_y + yp_x - \Omega

where p is the particles momentum and \Omega = \frac{1-mu}{r_1} + \frac{\mu}{r_2} . Note that I'm not involving the z-component yet until I get the basics sorted out.

Now this can be written as the following set of first order differential equations...

\dot{x} = p_x + y
\dot{p_x} = p_y - \frac{\partial \Omega}{\partial x}
\dot{y} = p_y - x
\dot{p_y} = -p_x - \frac{\partial \Omega}{\partial y}

(I prefer to write it like this, keeping x parts together and so on rather than keeping position parts together etc...)

Now this is where I get a bit stuck as I'm a bit confused about how to deal the the particles momentum. This is the code I wrote for the function...

function xdot = symCR3BP(t,x,mu)
xdot = zeros(4,1);
xdot(1) = x(2) + x(3);
xdot(2) = x(4) - (1-mu)*(x(1) + mu)/(((x(1) + mu)^2 + x(3)^2)^1.5) - mu*(x(1)-1+mu)/(((x(1) - (1 - mu))^2 + x(3)^2)^1.5);
xdot(3) = x(4) - x(1);
xdot(4) = -x(2) - (1-mu)*x(3)/(((x(1) + mu)^2 + x(3)^2)^1.5) - mu*x(3)/(((x(1) - (1 - mu))^2 + x(3)^2)^1.5);
end


Now notice I have never defined the particles mass and I'm not quite sure where to put it! An integration of this with a symplectic integrator just produces rubbish and any attempt to include mass in this yields momentums that are so small that they don't contribute to anything.

So where do I put the particles mass? I have tried to include it initial conditions but that did not work either. My usual test for integrators and this problem is to start near the L4 or L5 points and doing that with the above quickly sends the particle shooting out in an ever growing spiral which does not follow my other results. (my initial condition in the rotating frame is [0.5;0;0.9;0] which produces some form of horseshoe/tadpole type orbit for the standard 3 body problem function using any of my integrators.)

If anyones needs any more information just ask.
Thanks.

 

[mikeph]

I studied this a while back, but my code wasn't anything like this so I don't know how much help I can be.

The first thing I wonder is why you're even using momentum? In the circular restricted 3 body problem, the mass of the 3rd body is negligible, so it doesn't affect the paths of the massive bodies. Whatever its mass actually is doesn't matter because as long as the inertial mass = the gravitational mass, the acceleration of the body is the same.

I don't see any reference to momentum in your link so it confuses me. What is your overall objective? TO demonstrate the existence of stationary solutions? Or the stability of these solutions? This is all possible analytically. If you want to make simulations of unstable orbits you can just use some basic (or complex) predictions, eg. new_position = old_position + speed*timestep + 0.5*acceleration*timestep^2, or something like that (I've forgotten the details). A suitably complex version of this code will produce thousands of orbits before there is any observable divergence.

 

[Deadstar]

Hey thanks for the reply sorry for the delay in replying back!

My goal was to create a symplectic integrator in Matlab which I have done. It has worked fine with every other function I have used.

Then I turned my attentions to the three body problem and wondered if I could do a symplectic integration of that. Now perhaps I just plainly don't understand something but I'm under the impression that I should convert the 3BP into a set of Hamiltonian equations (shown in first post). Then I can use them in my symplectic integrator such that my particles trajectory conserves the Jacobi constant.

I have been reading though that because of the coordinate system... The Hamiltonian of the 3BP does not conserve energy so perhaps this is all fruitless...


I have already coded the 3BP in the usual sense and have gotten it to work in that it produces the expected particle trajectories (horseshoe orbits, orbits around L4/5 etc...) however the Jacobi constant always ends up increasing and the particle starts shooting of to infinity after some time.

My goal is...
Perform a symplectic integration of the circular restricted three body problem such that this does not happen (or at least not within large integration limits).

My problems are...
How to involve momentum.
Can it be done using an implicit Runge Kutta.

 

[mikeph]

I'm not surprised it flies off. Some orbits will really appear stable then after a while, fly off. Eg. how did the Trojan asteroids get there in the first place? Many probably were passing by and became caught in that orbit by the immense gravity of Jupiter/sun. If this can happen in dynamics, the time-reversal is also possible - orbits spending a long time at L4/5 and then flying off.

The alternative might be because of the chaotic nature of the system- it is unpredictable, though still deterministic, so your calculation will always diverge to some degree. My simulations never perfectly conserved energy either.I don't know how I can help with symplectic integrators, I never studied them, so that might explain it, but I can't understand why you're trying to involve momentum. The 3rd body has effectively zero mass, so it sounds like you are trying to force a particular solution (using Hamiltonians for problems involving momentum) on a certain problem. In fact this is explicitly stated in your goal: "solve X using Y". To me at least, I cannot see how Y will help you solve X.