0
|
6800 Views
|
2 Replies
|
0 Total Likes
View groups...
Share
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
2 Replies
Sort By:
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
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