newtonian mechanics - Elliptic Orbit Solution based on initial conditions

$$\ddot{\bf{r}}=-\frac{\bf{r}}{|r|}\frac{k}{|r|^2}$$

$k$ here is a constant dependent on the gravitational constant, and the masses of the two objects. If I transform it into cartesian coordinates:

$$\ddot{X}(t)=-\frac{k X(t)}{\left(X(t)^2+Y(t)^2\right)^{3/2}}$$

$$\ddot{Y}(t)=-\frac{k Y(t)}{\left(X(t)^2+Y(t)^2\right)^{3/2}}$$

I can not solve this system of equations. Perhaps it would require some special functions like the elliptic function etc. Maybe I should get rid of the time dependency and represent it as an implicit function but I do not know how. I realize that solving an elliptic orbit is know, I am curious about how one would solve it in the fundamental f=ma kind of way, without any other assumption.

Answer

First of all, it is best to solve the equation in polar coordinates rather than Cartesian coordinates: this works well with the fact that the force vector - and therefore the acceleration vector - is always pointing in a radial direction.

While it is possible to express the equations of motions as differential equations in time, it is not possible to obtain closed form solutions of the orbital position as a function of time, and so these time differential equations cannot be solved as they as. The strategy is to reformulate the time differential equations into "positional" differential equations, i.e. the derivative terms are differentiated with respect to the orbital angular position. Details for this are provided below.

Formulation of polar equation of motion in time

$r(t)$ and $\theta(t)$ are the polar coordinates. It is useful to be able to relate the acceleration vector to the polar coordinates and their derivatives. The position vector can be represented as:

$$\mathbf{r}(t) = \begin{bmatrix} r\cos{\theta}\\ r\sin{\theta} \end{bmatrix}$$

The acceleration is obtain by double differentiation in time:

$$\ddot{\mathbf{r}}(t) = \begin{bmatrix} \ddot{r}\cos{\theta} - 2\dot{r}\dot{\theta}\sin{\theta} - r\dot{\theta}^2\cos{\theta} - r\ddot{\theta}\sin{\theta}\\ \ddot{r}\sin{\theta} + 2\dot{r}\dot{\theta}\cos{\theta} - r\dot{\theta}^2\sin{\theta} + r\ddot{\theta}\cos{\theta} \end{bmatrix} = \left(\ddot{r} - r\dot{\theta}^2\right) \begin{bmatrix} r\cos{\theta}\\ r\sin{\theta} \end{bmatrix} + \left(2\dot{r}\dot{\theta} + r\ddot{\theta}\right) \begin{bmatrix} -r\sin{\theta}\\ r\cos{\theta} \end{bmatrix}$$

Note how $\ddot{r} - r\dot{\theta}^2$ is the radial component of acceleration and $2\dot{r}\dot{\theta} + r\ddot{\theta}$ is the tangential component. However, we know that gravity only imparts a radial acceleration on the body, and so two simultaneous equations of motion are obtained:

$$\ddot{r} - r\dot{\theta}^2 = -\frac{k}{r^2}$$

$$2\dot{r}\dot{\theta} + r\ddot{\theta} = 0$$

These equations cannot be immediately solved as they are, since both equations contain $r$ and $\theta$ terms together. However, it can be shown that the second equation amounts to the conservation of angular momentum. Multiply the second equation by $r$ to yield

$$2 r \dot{r}\dot{\theta} + r^2\ddot{\theta} = 0$$

Noting that $\frac{\mathrm{d}}{\mathrm{d}t}r^2 = 2r\dot{r}$, the equation can be re-expressed as

$$\frac{\mathrm{d}}{\mathrm{d}t}\left(r^2\right)\dot{\theta} + r^2\frac{\mathrm{d}}{\mathrm{d}t}\left(\dot{\theta}\right) = 0 \longrightarrow \frac{\mathrm{d}}{\mathrm{d}t}\left(r^2 \dot{\theta}\right) = 0$$

This means that

$$r^2\dot{\theta} = \textrm{constant}$$

We can represent this constant quantity with $h = r^2\dot{\theta}$, and $h$ is the specific angular momentum. Now we can substituted $\dot{\theta} = h/r^2$ into the first equation of motion, yielding

$$\ddot{r} - \frac{h^2}{r^3} + \frac{k}{r^2} = 0$$

This is the polar equation of motion in time, but this form is still not sufficient in determining the shape of an orbit.

Obtaining the "positional" equation of motion

A useful step to perform is to defined a new quantity $q = \frac{1}{r}$, and substitute this into the equation of motion. In the process, we can convert time derivatives of $r(t)$ into positional derivatives of $q(\theta)$

Note how $\dot{r}$ can be expanded using the chain rule:

$$\dot{r} = \frac{\mathrm{d}r}{\mathrm{d}t} = \frac{\mathrm{d}r}{\mathrm{d}\theta}\frac{\mathrm{d}\theta}{\mathrm{d}t} = \frac{\mathrm{d}}{\mathrm{d}\theta}\left(\frac{1}{q}\right)\dot{\theta} = -\frac{1}{q^2}\frac{\mathrm{d}q}{\mathrm{d}\theta}hq^2 = -h\frac{\mathrm{d}q}{\mathrm{d}\theta}$$

Likewise, $\ddot{r}$ becomes

$$\ddot{r} = \frac{\mathrm{d}\dot{r}}{\mathrm{d}t} = \frac{\mathrm{d}\dot{r}}{\mathrm{d}\theta}\frac{\mathrm{d}\theta}{\mathrm{d}t} = \frac{\mathrm{d}}{\mathrm{d}\theta}\left(-h\frac{\mathrm{d}q}{\mathrm{d}\theta}\right)\dot{\theta} = -h^2 q^2\frac{\mathrm{d}^2 q}{\mathrm{d}\theta^2}$$

Therefore, the equation of motion can be re-expressed as:

$$-h^2 q^2 \frac{\mathrm{d}^2 q}{\mathrm{d}\theta^2} - h^2 q^3 + k q^2 = 0$$

or

$$\frac{\mathrm{d}^2 q}{\mathrm{d}\theta^2} + q = \frac{k}{h^2}$$

This equation can be solved, and is in fact a linear differential equation with a general solution:

$$q(\theta) = A\cos{\theta} + B\sin{\theta} + \frac{k}{h^2}$$

Solving the equation of motion

$A$ and $B$ are yet to be known, but can be determined with two boundary conditions. $h$ is also unknown, and can be determined with a third boundary condition.

Let's say that $r_0$ is the initial distance of the body from the orbited planet, $u_0$ is the initial (outward) radial velocity, and $v_0$ is the initial (anticlockwise) tangential velocity. With no loss of generality, we can assume the orbiting body is initially at $\theta = 0$.

The first boundary condition is that $r = r_0$ at $t=0$. In terms of $q$, this becomes $q = \frac{1}{r_0}$ at $\theta = 0$.

The second boundary condition involves the radial velocity: $\dot{r} = u_0$ at $t = 0$. This becomes $\mathrm{d}q/\mathrm{d}\theta = -u_0/h$ at $\theta = 0$.

As for the final boundary condition, note that $h = r^2\dot{\theta}$. Then $h$ can be determined from the initial conditions, such that $h = r_0 v_0$.

With these conditions, the expression for $q$ becomes:

$$q(\theta) = \left(\frac{1}{r_0} - \frac{k}{r_0^2 v_0^2}\right)\cos{\theta} - \frac{u_0}{r_0 v_0}\sin{\theta} + \frac{k}{r_0^2 v_0^2}$$

Then, since $r = 1/q$:

$$r(\theta) = \frac{r_0}{\left(1 - \frac{k}{r_0 v_0^2}\right)\cos{\theta} - \frac{u_0}{v_0}\sin{\theta} + \frac{k}{r_0 v_0^2}}$$

This equation is the polar equation of the orbit, and is the equation of an ellipse for certain choices of initial conditions $r_0$, $u_0$ and $v_0$. This equation will also account for parabolic and hyperbolic "orbits".