# Orbit integrator / Keplerian orbit

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?

orbit[x0_,y0_,a_,b_]:=Module[{sol,x,y,Gm,t},
Gm=1;
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

a=20 and b=1

a=1 and b=20 (no change)

• 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 Nov 25 at 13:52
• What are the parameters a and b supposed to physically represent? Nov 25 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


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]


orbit[1, 0, 0, 1.5]


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.

• 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.
– Nick
Nov 25 at 13:47
• 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.
– Nick
Nov 25 at 17:17
• 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. Nov 26 at 9:02