# Get the trajectory of an object in a gravitational field

Given the initial position and velocity, I want get the trajectory of an object in a gravitational field. Instead of using Kepler's Laws, I want to solve the following differential equations: DSolve[{x == a, y == b, Derivative[x] == c,
Derivative[y] ==
d, (x^\[Prime]\[Prime])[t] == -((G M x[t])/(x[t]^2 + y[t]^2)^(
3/2)), (y^\[Prime]\[Prime])[t] == -((G M y[t])/(x[t]^2 + y[t]^2)^(
3/2))}, {x, y}, t]


But instead of solving it, Mathematica just returns the expression itself. Any help?

• Change the second derivatives to  x''[t] and y''[t] . Still Mathematica isn't able to solve this odes. Jun 11, 2021 at 9:05
• This problem can be solved with a more systematic approach transforming slightly given differential equations. Take a look at analogous problem solved exactly within general relativity: The time-like geodesics (orbits) in the Schwarzschild spacetime Jun 11, 2021 at 11:40
• It is well known (Newton) that this problem can be solved in polarcoordinates with the additional constraint "Angular momentum constant" Jun 11, 2021 at 11:48
• @UlrichNeumann: To the best of my knowledge you can only get $r(\theta)$ that way. The OP appears to want $r(t)$ and $\theta(t)$ (or $x(t)$ and $y(t)$), and to the best of my knowledge there is not a known closed-form solution for that. I suspect they will need to resort to using NDSolve. Jun 11, 2021 at 17:41
• @MichaelSeifert You're right, thanks for your hint. Jun 11, 2021 at 20:10

After correcting several syntax errors, you may get a numerical solution by e.g. with arbitrary values for G and M:

G = 1;
M = 1;
a = 1; b = 0;
c = 0; d = 0.5;
tmax = 3;
sol = NDSolve[{x == a, y == b, x' == c, y' == d,
x''[t] == -((G M x[t])/(x[t]^2 + y[t]^2)^(3/2)),
y''[t] == -((G M y[t])/(x[t]^2 + y[t]^2)^(3/2))}, {x, y}, {t, 0,
tmax}];

ParametricPlot[Evaluate[{x[t], y[t]} /. sol], {t, 0, tmax}] Or an ubounded orbit with a bit more initial velocity:

G = 1;
M = 1;
a = 1; b = 0;
c = 0; d = 1.5;
tmax = 10;
sol = NDSolve[{x == a, y == b, x' == c, y' == d,
x''[t] == -((G M x[t])/(x[t]^2 + y[t]^2)^(3/2)),
y''[t] == -((G M y[t])/(x[t]^2 + y[t]^2)^(3/2))}, {x, y}, {t, 0,
tmax}];

ParametricPlot[Evaluate[{x[t], y[t]} /. sol], {t, 0, tmax}] 