# What's wrong with my code? (Modelling planetary motion) **URL:** \ **Category:** Questions\ **Tags:** astronomy\ **Created:** [July 7, 2015, 2:44pm UTC](https://community.wolfram.com/t/whats-wrong-with-my-code-modelling-planetary-motion/6107 "2015-07-07T14:44:29Z")\ **Posts on this page:** 3\ **Page:** 1
**Author:** ![miriam98ahmed](https://avatars.discourse-cdn.com/v4/letter/m/b487fb/32.png) [@miriam98ahmed](https://community.wolfram.com/u/miriam98ahmed)\ **Post date:** [July 7, 2015, 2:44pm UTC](https://community.wolfram.com/t/whats-wrong-with-my-code-modelling-planetary-motion/6107/1 "2015-07-07T14:44:29Z")
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 ---
**Author:** ![mthiel](https://yyz1.discourse-cdn.com/wolfram/user_avatar/community.wolfram.com/mthiel/32/42877_2.png) [@mthiel](https://community.wolfram.com/u/mthiel)\ **Post date:** [July 7, 2015, 3:11pm UTC](https://community.wolfram.com/t/whats-wrong-with-my-code-modelling-planetary-motion/6107/2 "2015-07-07T15:11:29Z")
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 ---
**Author:** ![mthiel](https://yyz1.discourse-cdn.com/wolfram/user_avatar/community.wolfram.com/mthiel/32/42877_2.png) [@mthiel](https://community.wolfram.com/u/mthiel)\ **Post date:** [July 7, 2015, 4:24pm UTC](https://community.wolfram.com/t/whats-wrong-with-my-code-modelling-planetary-motion/6107/3 "2015-07-07T16:24:39Z")
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}] ``` ![ScreenShot2015-07-07at17.09.20.png](https://global.discourse-cdn.com/wolfram/original/3X/2/1/21c5a72f189f189a7448c5ef082bd3e97fb00cce.png) 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}}] ``` ![ScreenShot2015-07-07at17.11.02.png](https://global.discourse-cdn.com/wolfram/original/3X/d/b/dbba641b859bb060c37a0963ed2aa4232dfcfc08.png) 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