I want to integrate a circular Keplerian orbit with GM = 1, a = 1 over 100 orbital periods. But the thing is that I can only make the Orbit change when I change a. And the Orbit remains the same when I change b..... Can someone see my mistake?

sol = NDSolve[{

x''[t] == -((Gm * x[t])/(x[t]^2 + y[t]^2)^(3/2)),
y''[t] == -((Gm * y[t])/(x[t]^2 + y[t]^2)^(3/2)),

x[0] == x0,
y[0] == y0,

x'[0] == Sqrt[Gm (2/Norm[{x[0],y[0]}] - 1/a)] * Normalize[{-a*Sin[ArcTan[x[0],y[0]]], b*Cos[ArcTan[x[0],y[0]]]}][[1]],
y'[0] == Sqrt[Gm (2/Norm[{x[0],y[0]}] - 1/a)] * Normalize[{-a*Sin[ArcTan[x[0],y[0]]], b*Cos[ArcTan[x[0],y[0]]]}][[2]]


{x[t],y[t]},{t,0, 100 * 2 Pi Sqrt[a^3/(Gm)]}   ];

ParametricPlot[Evaluate[{x[t], y[t]} /. sol], {t, 0, 100*2 Pi Sqrt[a^3/(Gm)]}]


a=1 and b=1
enter image description here

a=20 and b=1 enter image description here

a=1 and b=20 (no change)

enter image description here

  • $\begingroup$ You only show function definition but do not show the call itself. Better to show the call also so others do not have to guess how you called the function orbit $\endgroup$
    – Nasser
    Nov 25 '21 at 13:52
  • $\begingroup$ What are the parameters a and b supposed to physically represent? $\endgroup$ Nov 25 '21 at 15:48

"b" only appears in the initial conditions for y' and y'. And only in the Term "Normalize[..]". Now, look at this term:

Normalize[{-a*Sin[ArcTan[x[0], y[0]]], 
   b*Cos[ArcTan[x[0], y[0]]]}] // Simplify

enter image description here

You will not that this term is independent of "b".


I simplified your program to show the essential:

orbit[x0_, y0_, v0x_, v0y_] := Module[{sol, x, y, t},
  sol = NDSolve[{x''[t] == -((x[t])/(x[t]^2 + y[t]^2)^(3/2)), 
     y''[t] == -((y[t])/(x[t]^2 + y[t]^2)^(3/2)), x[0] == x0, 
     y[0] == y0,
     x'[0] == v0x, y'[0] == v0y},
    {x[t], y[t]},
    {t, 0, 100}];
  ParametricPlot[Evaluate[{x[t], y[t]} /. sol], {t, 0, 100 }]

Now if we choose x0=1, y0=0, v0x=0,x0y=v the path will have a major axis along the x-axis. {1,0} will either be the point of the trajectory farthest, or for larger velocities, nearest point to the origin. For still higher velocities, the body escapes on an hyperbola:

orbit[1, 0, 0, 0.5]

orbit[1, 0, 0, 1.3]

enter image description here

orbit[1, 0, 0, 1.5]

enter image description here

Note also, if the body comes close to the sun, the acceleration will be large and you get numerical problems. You need to decrease the step size and increase orecision.

  • $\begingroup$ Thank you. But how can I create an elliptic Orbit, where I can change the values for b? My first idea was to set x[t]= a cos (t) and y[t]= b sin (t ), but when I do this, my Plot becomes very odd. $\endgroup$
    – Nick
    Nov 25 '21 at 13:47
  • $\begingroup$ Thank you for the explanation. But I have the problem, if I change the initial conditions, so set y0 to 0 and x0 to 1 then I get only error messages out. And I do not understand why. $\endgroup$
    – Nick
    Nov 25 '21 at 17:17
  • $\begingroup$ I deleted my comment because it is for the case where the center of the ellipse is at the origin. But in your case, the sun is at the origin. Instead I added a part to my answer where I simplified your input to the essential. $\endgroup$ Nov 26 '21 at 9:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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