0
$\begingroup$

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:

enter image description here

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?

Improve this question
$\endgroup$
5
  • $\begingroup$ Change the second derivatives to ` x''[t] and y''[t]` . Still Mathematica isn't able to solve this odes. $\endgroup$
    – Ulrich Neumann
    Jun 11, 2021 at 9:05
  • $\begingroup$ 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 $\endgroup$
    – Artes
    Jun 11, 2021 at 11:40
  • $\begingroup$ It is well known (Newton) that this problem can be solved in polarcoordinates with the additional constraint "Angular momentum constant" $\endgroup$
    – Ulrich Neumann
    Jun 11, 2021 at 11:48
  • 1
    $\begingroup$ @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. $\endgroup$
    – Michael Seifert
    Jun 11, 2021 at 17:41
  • $\begingroup$ @MichaelSeifert You're right, thanks for your hint. $\endgroup$
    – Ulrich Neumann
    Jun 11, 2021 at 20:10

1 Answer 1

Reset to default
4
$\begingroup$

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}]

enter image description here

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}]

enter image description here

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.