Message Boards

WOLFRAM COMMUNITY

6452 Views

2 Replies

0 Total Likes

View groups...

Follow this post

Share this post:

GROUPS:

Astronomy

What's wrong with my code? (Modelling planetary motion)

Miriam Ahmed

Posted 9 years ago

Differential equations for modelling planetary motion (3 bodies), getting errors saying that it is 'null'- anyone point out where I've gone wrong? sols = NDSolve[ {x1''[t] == (G m2 (x2 - x1[t]))/((x2 - x1[t])^2 - (y2 - y1)^2)^(3/ 2) + (G m3 (x3 - x1[t]))/((x3 - x1[t])^2 - (y3 - y1)^2)^(3/2), y1''[t] == (G m2 (y2 - y1[t]))/((x2 - x1)^2 - (y2 - y1[t])^2)^(3/ 2) + (G m3 (y3 - y1[t]))/((x3 - x1)^2 - (y3 - y1[t])^2)^(3/2), x2''[t] == (G m1 (x1 - x2[t]))/((x1 - x2[t])^2 - (y1 - y2)^2)^(3/ 2) + (G m3 (x3 - x2[t]))/((x3 - x2[t])^2 - (y3 - y2)^2)^(3/2), y2''[t] == (G m1 (y1 - y2[t]))/((x1 - x2)^2 - (y1 - y2[t])^2)^(3/ 2) + (G m3 (y3 - y2[t]))/((x3 - x2)^2 - (y3 - y2[t])^2)^(3/2), x3''[t] == (G m1 (x1 - x3[t]))/((x1 - x3[t])^2 - (y1 - y3)^2)^(3/ 2) + (G m2 (x2 - x3[t]))/((x2 - x3[t])^2 - (y2 - y3)^2)^(3/2), y3''[t] == (G m1 (y1 - y3[t]))/((x1 - x3)^2 - (y1 - y3[t])^2)^(3/ 2) + (G m2 (y2 - y3[t]))/((x2 - x3)^2 - (y2 - y3[t])^2)^(3/2), x1[0] == d, y1[0] == 0, x1'[0] == -0.8, y1'[0] == v, x2[0] == d, y2[0] == 0, x2'[0] == -0.8, y2'[0] == v, x3[0] == d, y3[0] == 0, x3'[0] == -0.8, y3'[0] == v} /. {G -> 1, m1 -> 1, m2 -> 1, m3 -> 1, d -> 1, v -> 1,}, {x1[t], y[1 t], x2[t], y2[t], x3[t], y3[t]}, {t, 0, 100}]; (if you paste it into mathematica hopefully it'll look nicer), also I'm a sixteen year old girl so please go easy

POSTED BY: Miriam Ahmed

2 Replies

Sort By:

Marco Thiel

Marco Thiel, University of Aberdeen - Dept. of Physics/Mathematics

Posted 9 years ago

Hi, the following equations will produce some output. sols = NDSolve[{x1''[ t] == (G m2 (x2[t] - x1[t]))/((x2[t] - x1[t])^2 - (y2[t] - y1[t])^2)^(3/ 2) + (G m3 (x3[t] - x1[t]))/((x3[t] - x1[t])^2 - (y3[t] - y1[t])^2)^(3/2), y1''[t] == (G m2 (y2[t] - y1[t]))/((x2[t] - x1[t])^2 - (y2[t] - y1[t])^2)^(3/ 2) + (G m3 (y3[t] - y1[t]))/((x3[t] - x1[t])^2 - (y3[t] - y1[t])^2)^(3/2), x2''[t] == (G m1 (x1[t] - x2[t]))/((x1[t] - x2[t])^2 - (y1[t] - y2[t])^2)^(3/ 2) + (G m3 (x3[t] - x2[t]))/((x3[t] - x2[t])^2 - (y3[t] - y2[t])^2)^(3/2), y2''[t] == (G m1 (y1[t] - y2[t]))/((x1[t] - x2[t])^2 - (y1[t] - y2[t])^2)^(3/ 2) + (G m3 (y3[t] - y2[t]))/((x3[t] - x2[t])^2 - (y3[t] - y2[t])^2)^(3/2), x3''[t] == (G m1 (x1[t] - x3[t]))/((x1[t] - x3[t])^2 - (y1[t] - y3[t])^2)^(3/ 2) + (G m2 (x2[t] - x3[t]))/((x2[t] - x3[t])^2 - (y2[t] - y3[t])^2)^(3/2), y3''[t] == (G m1 (y1[t] - y3[t]))/((x1[t] - x3[t])^2 - (y1[t] - y3[t])^2)^(3/ 2) + (G m2 (y2[t] - y3[t]))/((x2[t] - x3[t])^2 - (y2[t] - y3[t])^2)^(3/2), x1[0] == d1, y1[0] == 0, x1'[0] == -0.8, y1'[0] == v1, x2[0] == d2, y2[0] == 0, x2'[0] == -0.2, y2'[0] == v2, x3[0] == d3, y3[0] == 0, x3'[0] == -0.1, y3'[0] == v3} /. {G -> 0.1, m1 -> 1, m2 -> 1, m3 -> 1, d1 -> 10, v1 -> 1, d2 -> 20, v2 -> 2, d3 -> 30, v3 -> 1}, {x1[t], y1[ t], x2[t], y2[t], x3[t], y3[t]}, {t, 0, 10}]; You can plot functions like this: Plot[{x1[t], x2[t]} /. sols, {t, 0, 10}] You can also use ParametricPlot: ParametricPlot[{{x1[t], y1[t]} /. sols, {x2[t], y2[t]} /. sols, {x3[t], y3[t]} /. sols}, {t, 0, 10}, PlotRange -> {{-30, 30}, {-30, 30}}] The results are still not what you would expect. This is because your equations for the gravitational force are not correct. Where you divide by the squared distances in the x and y direction you should actually add not subtract: sols = NDSolve[{x1''[ t] == (G m2 (x2[t] - x1[t]))/((x2[t] - x1[t])^2 + (y2[t] - y1[t])^2)^(3/ 2) + (G m3 (x3[t] - x1[t]))/((x3[t] - x1[t])^2 + (y3[t] - y1[t])^2)^(3/2), y1''[t] == (G m2 (y2[t] - y1[t]))/((x2[t] - x1[t])^2 + (y2[t] - y1[t])^2)^(3/ 2) + (G m3 (y3[t] - y1[t]))/((x3[t] - x1[t])^2 + (y3[t] - y1[t])^2)^(3/2), x2''[t] == (G m1 (x1[t] - x2[t]))/((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2)^(3/ 2) + (G m3 (x3[t] - x2[t]))/((x3[t] - x2[t])^2 + (y3[t] - y2[t])^2)^(3/2), y2''[t] == (G m1 (y1[t] - y2[t]))/((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2)^(3/ 2) + (G m3 (y3[t] - y2[t]))/((x3[t] - x2[t])^2 + (y3[t] - y2[t])^2)^(3/2), x3''[t] == (G m1 (x1[t] - x3[t]))/((x1[t] - x3[t])^2 + (y1[t] - y3[t])^2)^(3/ 2) + (G m2 (x2[t] - x3[t]))/((x2[t] - x3[t])^2 + (y2[t] - y3[t])^2)^(3/2), y3''[t] == (G m1 (y1[t] - y3[t]))/((x1[t] - x3[t])^2 + (y1[t] - y3[t])^2)^(3/ 2) + (G m2 (y2[t] - y3[t]))/((x2[t] - x3[t])^2 + (y2[t] - y3[t])^2)^(3/2), x1[0] == d1, y1[0] == 0, x1'[0] == -0.8, y1'[0] == v1, x2[0] == d2, y2[0] == 0, x2'[0] == -0.2, y2'[0] == v2, x3[0] == d3, y3[0] == 0, x3'[0] == -0.1, y3'[0] == v3} /. {G -> 2, m1 -> 1, m2 -> 1, m3 -> 1, d1 -> 1, v1 -> 0.1, d2 -> 2, v2 -> 0.2, d3 -> 3, v3 -> 0.4}, {x1[t], y1[ t], x2[t], y2[t], x3[t], y3[t]}, {t, 0, 100}]; You will notice that this is probably still not what you want, or is it? Cheers, Marco

Hi,

the following equations will produce some output.

sols = NDSolve[{x1''[
       t] == (G m2 (x2[t] - 
            x1[t]))/((x2[t] - x1[t])^2 - (y2[t] - y1[t])^2)^(3/
           2) + (G m3 (x3[t] - 
            x1[t]))/((x3[t] - x1[t])^2 - (y3[t] - y1[t])^2)^(3/2), 
     y1''[t] == (G m2 (y2[t] - 
            y1[t]))/((x2[t] - x1[t])^2 - (y2[t] - y1[t])^2)^(3/
           2) + (G m3 (y3[t] - 
            y1[t]))/((x3[t] - x1[t])^2 - (y3[t] - y1[t])^2)^(3/2), 
     x2''[t] == (G m1 (x1[t] - 
            x2[t]))/((x1[t] - x2[t])^2 - (y1[t] - y2[t])^2)^(3/
           2) + (G m3 (x3[t] - 
            x2[t]))/((x3[t] - x2[t])^2 - (y3[t] - y2[t])^2)^(3/2), 
     y2''[t] == (G m1 (y1[t] - 
            y2[t]))/((x1[t] - x2[t])^2 - (y1[t] - y2[t])^2)^(3/
           2) + (G m3 (y3[t] - 
            y2[t]))/((x3[t] - x2[t])^2 - (y3[t] - y2[t])^2)^(3/2), 
     x3''[t] == (G m1 (x1[t] - 
            x3[t]))/((x1[t] - x3[t])^2 - (y1[t] - y3[t])^2)^(3/
           2) + (G m2 (x2[t] - 
            x3[t]))/((x2[t] - x3[t])^2 - (y2[t] - y3[t])^2)^(3/2), 
     y3''[t] == (G m1 (y1[t] - 
            y3[t]))/((x1[t] - x3[t])^2 - (y1[t] - y3[t])^2)^(3/
           2) + (G m2 (y2[t] - 
            y3[t]))/((x2[t] - x3[t])^2 - (y2[t] - y3[t])^2)^(3/2), 
     x1[0] == d1, y1[0] == 0, x1'[0] == -0.8, y1'[0] == v1, 
     x2[0] == d2, y2[0] == 0, x2'[0] == -0.2, y2'[0] == v2, 
     x3[0] == d3, y3[0] == 0, x3'[0] == -0.1, 
     y3'[0] == v3} /. {G -> 0.1, m1 -> 1, m2 -> 1, m3 -> 1, d1 -> 10, 
     v1 -> 1, d2 -> 20, v2 -> 2, d3 -> 30, v3 -> 1}, {x1[t], y1[ t], 
    x2[t], y2[t], x3[t], y3[t]}, {t, 0, 10}];

You can plot functions like this:

Plot[{x1[t], x2[t]} /. sols, {t, 0, 10}]

enter image description here

You can also use ParametricPlot:

ParametricPlot[{{x1[t], y1[t]} /. sols, {x2[t], y2[t]} /. sols, {x3[t], y3[t]} /. sols}, {t, 0, 10}, PlotRange -> {{-30, 30}, {-30, 30}}]

enter image description here

The results are still not what you would expect. This is because your equations for the gravitational force are not correct. Where you divide by the squared distances in the x and y direction you should actually add not subtract:

sols = NDSolve[{x1''[
       t] == (G m2 (x2[t] - 
            x1[t]))/((x2[t] - x1[t])^2 + (y2[t] - y1[t])^2)^(3/
           2) + (G m3 (x3[t] - 
            x1[t]))/((x3[t] - x1[t])^2 + (y3[t] - y1[t])^2)^(3/2), 
     y1''[t] == (G m2 (y2[t] - 
            y1[t]))/((x2[t] - x1[t])^2 + (y2[t] - y1[t])^2)^(3/
           2) + (G m3 (y3[t] - 
            y1[t]))/((x3[t] - x1[t])^2 + (y3[t] - y1[t])^2)^(3/2), 
     x2''[t] == (G m1 (x1[t] - 
            x2[t]))/((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2)^(3/
           2) + (G m3 (x3[t] - 
            x2[t]))/((x3[t] - x2[t])^2 + (y3[t] - y2[t])^2)^(3/2), 
     y2''[t] == (G m1 (y1[t] - 
            y2[t]))/((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2)^(3/
           2) + (G m3 (y3[t] - 
            y2[t]))/((x3[t] - x2[t])^2 + (y3[t] - y2[t])^2)^(3/2), 
     x3''[t] == (G m1 (x1[t] - 
            x3[t]))/((x1[t] - x3[t])^2 + (y1[t] - y3[t])^2)^(3/
           2) + (G m2 (x2[t] - 
            x3[t]))/((x2[t] - x3[t])^2 + (y2[t] - y3[t])^2)^(3/2), 
     y3''[t] == (G m1 (y1[t] - 
            y3[t]))/((x1[t] - x3[t])^2 + (y1[t] - y3[t])^2)^(3/
           2) + (G m2 (y2[t] - 
            y3[t]))/((x2[t] - x3[t])^2 + (y2[t] - y3[t])^2)^(3/2), 
     x1[0] == d1, y1[0] == 0, x1'[0] == -0.8, y1'[0] == v1, 
     x2[0] == d2, y2[0] == 0, x2'[0] == -0.2, y2'[0] == v2, 
     x3[0] == d3, y3[0] == 0, x3'[0] == -0.1, 
     y3'[0] == v3} /. {G -> 2, m1 -> 1, m2 -> 1, m3 -> 1, d1 -> 1, 
     v1 -> 0.1, d2 -> 2, v2 -> 0.2, d3 -> 3, v3 -> 0.4}, {x1[t], 
    y1[ t], x2[t], y2[t], x3[t], y3[t]}, {t, 0, 100}];

You will notice that this is probably still not what you want, or is it?

Cheers,

Marco

POSTED BY: Marco Thiel

Marco Thiel

Marco Thiel, University of Aberdeen - Dept. of Physics/Mathematics

Posted 9 years ago

Hi Miriam, it would be helpful tonpost the code in one of the "code-windows". There is a button when you compose your message. Glancing over the code it seems that there are still x and y variables that do not have the time dependence. So x1 should be x1[t]. Cheers, Marco

POSTED BY: Marco Thiel

Reply to this discussion

Reply Preview

Attachments

Remove Add a file to this post

Follow this discussion

or Discard

Group Abstract

Feedback