# Can Someone Help Me WIth My Code?

5 views (last 30 days)
Craig on 15 May 2013
I am trying to code an approximation for the three body problem, and I am having a lot of trouble.
I solved for the three equations of motion by hand and got
a1 = G*(m1*m2)/(r1)^3 + G*(m1*m3)/(r3)^3; a2 = G*(m1*m2)/(r1)^3 + G*(m2*m3)/(r2)^3; a3 = G*(m1*m3)/(r3)^3 + G*(m2*m3)/(r2)^3;
r1, r2, r3, being the distances between the three different bodies.
I want to add an origin (arbitrary or not, it doesn't matter because the bodies will move, they just need a starting reference point. Right?), but I am not sure how.
I would also like to plot each bodies path on the same graph, and if possible make an animation to show the bodies moving at each timestep (possibly the comet function?)
But that is only if I can get it to work first. I really need help with the basic plot. I don't think I did it right.
Any help would be really appreciated.
Here is my code:
function y = threebody(m1,x1,y1,v1,m2,x2,y2,v2,m3,x3,y3,v3)
r3 = ((x1-x3)^2+(y1+y3)^2)^(1/2);
r2 = ((x2-x3)^2+(y2+y3)^2)^(1/2);
r1 = ((x2-x1)^2+(y1-y2)^2)^(1/2);
G = 6.67384*10^(-11);
a1(1) = G*(m1*m2)/(r1(1))^3 + G*(m1*m3)/(r3(1))^3;
a2(1) = G*(m1*m2)/(r1(1))^3 + G*(m2*m3)/(r2(1))^3;
a3(1) = G*(m1*m3)/(r3(1))^3 + G*(m2*m3)/(r2(1))^3;
steps = 50;
%t = 0.0;
dt = 50/steps;
%d1 = v1*t+(.5)*a1*t^2;
%d2 = v2*t+(.5)*a2*t^2;
%d3 = v3*t+(.5)*a3*t^2;
for i = 1:steps
r1(i+1) = r1(i) - r1(i)*dt;
r2(i+1) = r2(i) - r2(i)*dt;
r3(i+1) = r3(i) - r3(i)*dt;
a1(i+1) = G*(m1*m2)/(r1(i+1))^3 + G*(m1*m3)/(r3(i+1))^3;
a2(i+1) = G*(m1*m2)/(r1(i+1))^3 + G*(m2*m3)/(r2(i+1))^3;
a3(i+1) = G*(m1*m3)/(r3(i+1))^3 + G*(m2*m3)/(r2(i+1))^3;
%d1(i+1) = v1*t+(.5)*a1(i+1)*t^2;
%d2(i+1) = v2*t+(.5)*a2(i+1)*t^2;
%d3(i+1) = v3*t+(.5)*a3(i+1)*t^2;
end
plot(r1,'-b');
hold on
plot(r2,'-g');
plot(r3,'-r');
end
Thanks again

Yao Li on 15 May 2013
You haven't defined the output y for the function. If you don't need an output, just remove 'y=' in the first line.
If it's possible, would u pls. give me a set of initial values of the inputs?
Also, if you have the simMechanics toolbox in simulink, it's easier to build the physical model with simMechanics.
Craig on 15 May 2013
Hello Yao Li,
I used (1.37e+22,1000,250,1023,5.972e+24,1050,300,29722,1.989e+30,0,0,0) for Initial conditions.
1.37e+22 being the mass of the moon and 1023 m/s being the velocity,
5.972e+24 being the mass of the Earth and 29722 being the velocity,
1.989e+30 being the mass of the Sun and 0 being its velocity.
The rest were just chosen at random because I am having trouble with the origin and setting the initial positions.
I know the position, velocity, and acceleration are changing for every time step, but I don't know how to model those changes. I know the for loop can update them with every iteration, but I am having trouble doing so.
I am still a beginner at Matlab and I don't know what simulink is sorry.

Roger Stafford on 15 May 2013
Edited: Roger Stafford on 15 May 2013
Those aren't valid equations for the accelerations of your three bodies under mutual gravitational attraction. Acceleration for your two-dimensional problem should be two-element vectors containing the acceleration components. Remember, acceleration is a vector.
a1 = G*m2/((x2-x1)^2+(y2-y1)^2)^(3/2)*[x2-x1,y2-y1] + ...
G*m3/((x3-x1)^2+(y3-y1)^2)^(3/2)*[x3-x1,y3-y1];
What you have in your code are scalar quantities - old Sir Issac would turn over in his grave.
Also to do the computation you should be using matlab's differential equation solver, ode45 or the like, for the integration to produce curves.
##### 2 CommentsShowHide 1 older comment
Roger Stafford on 15 May 2013
Edited: Roger Stafford on 15 May 2013
Those are closer to being right. However the vector part is wrong. It should be:
a1 = G*m2/r21^3*[x2-x1,y2-y1] + G*m3/r31^3*[x3-x1,y3-y1];
where
r21 = sqrt((x2-x1)^2+(y2-y1)^2)
r31 = sqrt((x3-x1)^2+(y3-y1)^2)
This is because the first term of the acceleration points in the direction of the force from point (x1,y1) toward point (x2,y2) and is inversely proportional to the square of the distance between them. Note that m1 is not present because the force is proportional to the mass m1, but force equals mass times acceleration so m1 drops out of the equation.