Get the trajectory of an object in a gravitational field

Question

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[0] == a, y[0] == b, Derivative[1][x][0] == c, 
  Derivative[1][y][0] == 
   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. — Ulrich Neumann, Jun 11 '21 at 09: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 — Artes, Jun 11 '21 at 11:40
It is well known (Newton) that this problem can be solved in polarcoordinates with the additional constraint "Angular momentum constant" — Ulrich Neumann, Jun 11 '21 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. — Michael Seifert, Jun 11 '21 at 17:41

score 4 · Accepted Answer · edited Jun 11 '21 at 17:45

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[0] == a, y[0] == b, x'[0] == c, y'[0] == 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[0] == a, y[0] == b, x'[0] == c, y'[0] == 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}]

Get the trajectory of an object in a gravitational field

1 Answers1