Message Boards Message Boards

0
|
6617 Views
|
2 Replies
|
0 Total Likes
View groups...
Share
Share this post:
GROUPS:

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

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

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

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
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract