Can Someone Help Me WIth My Code?

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

Answer 1

[Yao Li]

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]

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.

Answer 2

[Roger Stafford]

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.

[Craig]

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]

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.