Can Someone Help Me WIth My Code?

1 vue (au cours des 30 derniers jours)
Craig
Craig le 15 Mai 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

Réponses (2)

Yao Li
Yao Li le 15 Mai 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.
  1 commentaire
Craig
Craig le 15 Mai 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.

Connectez-vous pour commenter.


Roger Stafford
Roger Stafford le 15 Mai 2013
Modifié(e) : Roger Stafford le 15 Mai 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 commentaires
Craig
Craig le 15 Mai 2013
I did do the vectors and calculations by hand, however I wasn't sure how to put them into Matlab. I am still a beginner programer and don't really know how to do much of anything.
My equations of motion before putting them into matlab were:
a1 = G( m2 )[x1] + G( m3 )[x1]
(r1^3)[y1] (r3^3)[y1]
a2 = G( m1 )[x2] + G( m3 )[x2]
(r1^3)[y2] (r2^3)[y2]
a3 = G( m1 )[x3] + G( m2 )[x3]
(r3^3)[y3] (r2^3)[y3]
I may have messed them up when coding them.
As for the DE solver, I am not really sure how to use that either.
Roger Stafford
Roger Stafford le 15 Mai 2013
Modifié(e) : Roger Stafford le 15 Mai 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.

Connectez-vous pour commenter.

Catégories

En savoir plus sur Numerical Integration and Differential Equations dans Help Center et File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by