Confirming Einstein’s Theory of General Relativity With Calculus, Part 6c: Solving New Differential Equation with Variation of Parameters

In this series, I’m discussing how ideas from calculus and precalculus (with a touch of differential equations) can predict the precession in Mercury’s orbit and thus confirm Einstein’s theory of general relativity. The origins of this series came from a class project that I assigned to my Differential Equations students maybe 20 years ago.

Under general relativity, the motion of a planet around the Sun —in polar coordinates $(r,\theta)$ , with the Sun at the origin — satisfies the initial-value problem

$u''(\theta) + u(\theta) = \displaystyle \frac{1}{\alpha} + \delta \left( \frac{1 + \epsilon \cos \theta}{\alpha} \right)^2$ ,

$u(0) = \displaystyle \frac{1}{P}$ ,

$u'(0) = 0$ ,

where $u = \displaystyle \frac{1}{r}$ , $\displaystyle \frac{1}{\alpha} = \frac{GMm^2}{\ell^2}$ , $\delta = \displaystyle \frac{3GM}{c^2}$ , $G$ is the gravitational constant of the universe, $m$ is the mass of the planet, $M$ is the mass of the Sun, $\ell$ is the constant angular momentum of the planet, $c$ is the speed of light, and $P$ is the smallest distance of the planet from the Sun during its orbit (i.e., at perihelion).

We now take the perspective of a student who is taking a first-semester course in differential equations. There are two standard techniques for solving a second-order non-homogeneous differential equations with constant coefficients. One of these is the method of variation of parameters. First, we solve the associated homogeneous differential equation

$u''(\theta) + u(\theta) = 0$ .

The characteristic equation of this differential equation is $r^2 + 1 = 0$ , which clearly has the two imaginary roots $r = \pm i$ . Therefore, two linearly independent solutions of the associated homogeneous equation are $u_1(\theta) = \cos \theta$ and $u_2(\theta) = \sin \theta$ .

(As an aside, this is one answer to the common question, “What are complex numbers good for?” The answer is naturally above the heads of Algebra II students when they first encounter the mysterious number $i$ , but complex numbers provide a way of solving the differential equations that model multiple problems in statics and dynamics.)

According to the method of variation of parameters, the general solution of the original nonhomogeneous differential equation

$u''(\theta) + u(\theta) = g(\theta)$

$u(\theta) = f_1(\theta) u_1(\theta) + f_2(\theta) u_2(\theta)$ ,

where

$f_1(\theta) = -\displaystyle \int \frac{u_2(\theta) g(\theta)}{W(\theta)} d\theta$ ,

$f_2(\theta) = \displaystyle \int \frac{u_1(\theta) g(\theta)}{W(\theta)} d\theta$ ,

and $W(\theta)$ is the Wronskian of $u_1(\theta)$ and $u_2(\theta)$ , defined by the determinant

$W(\theta) = \displaystyle \begin{vmatrix} u_1(\theta) & u_2(\theta) \\ u_1'(\theta) & u_2'(\theta) \end{vmatrix} = u_1(\theta) u_2'(\theta) - u_1'(\theta) u_2(\theta)$ .

Well, that’s a mouthful.

The only good news is that $W$ is easy to compute. Since $u_1(\theta) = \cos \theta$ and $u_2(\theta) = \sin \theta$ , we have

$W(\theta) = (\cos \theta)(\cos \theta) - (\sin \theta)(-\sin \theta) = \cos^2 \theta + \sin^2 \theta = 1$

from the usual Pythagorean trigonometric identity. Therefore, the denominators in the integrals for $f_1(\theta)$ and $f_2(\theta)$ essentially disappear.

Unfortunately, computing $f_1(\theta)$ and $f_2(\theta)$ , using

$g(\theta) = \displaystyle \frac{1}{\alpha} + \delta \left( \frac{1 + \epsilon \cos \theta}{\alpha} \right)^2$ ,

is a beast, requiring the creative use of multiple trigonometric identities. We begin with $f_1(\theta)$ , using the substitution $t = \cos \theta$ :

$f_1(\theta) = -\displaystyle \int \left[ \frac{1}{\alpha} + \delta \left( \frac{1 + \epsilon \cos \theta}{\alpha} \right)^2 \right] \sin \theta \, d\theta$

$= \displaystyle \int \left[ \frac{1}{\alpha} + \delta \left( \frac{1 + \epsilon t}{\alpha} \right)^2 \right] \, dt$

$= \displaystyle \frac{t}{\alpha} + \frac{\alpha\delta}{3\epsilon} \left( \frac{1 + \epsilon t}{\alpha} \right)^3 + a$

$= \displaystyle \frac{\cos \theta}{\alpha} + \frac{\alpha\delta}{3\epsilon} \left( \frac{1 + \epsilon \cos \theta}{\alpha} \right)^3 + a$ ,

where we use $+a$ for the constant of integration instead of the usual $+C$ . Second,

$f_2(\theta) = \displaystyle \int \left[ \frac{1}{\alpha} + \delta \left( \frac{1 + \epsilon \cos \theta}{\alpha} \right)^2 \right] \cos\theta \, d\theta$ .

Unfortunately, this is not easily simplified with a substitution, so we have to expand the integrand:

$f_2(\theta) = \displaystyle \int \left[ \frac{\cos \theta}{\alpha} + \frac{\delta \cos \theta}{\alpha^2} + \frac{2 \delta \epsilon \cos^2 \theta}{\alpha^2} + \frac{\delta \epsilon^2 \cos^3 \theta}{\alpha^2} \right] \, d\theta$

$= \displaystyle \int \left[ \frac{\cos \theta}{\alpha} + \frac{\delta \cos \theta}{\alpha^2} + \frac{\delta \epsilon (1 + \cos 2 \theta)}{\alpha^2} + \frac{\delta \epsilon^2 \cos \theta \cos^2 \theta}{\alpha^2} \right] \, d\theta$

$= \displaystyle \int \left[ \frac{\cos \theta}{\alpha} + \frac{\delta \cos \theta}{\alpha^2} + \frac{\delta \epsilon}{\alpha^2} + \frac{\delta \epsilon \cos 2 \theta}{\alpha^2}+ \frac{\delta \epsilon^2 \cos \theta (1- \sin^2 \theta)}{\alpha^2} \right] \, d\theta$

$= \displaystyle \int \left[ \frac{\cos \theta}{\alpha} + \frac{\delta (1 + \epsilon^2) \cos \theta}{\alpha^2} + \frac{\delta \epsilon}{\alpha^2} + \frac{\delta \epsilon \cos 2 \theta}{\alpha^2} - \frac{\delta \epsilon^2 \cos \theta \sin^2 \theta}{\alpha^2} \right] \, d\theta$

$= \displaystyle \frac{\sin\theta}{\alpha} + \frac{\delta (1 + \epsilon^2) \sin\theta}{\alpha^2} + \frac{\delta \epsilon \theta}{\alpha^2} + \frac{\delta \epsilon \sin 2 \theta}{2\alpha^2} - \frac{\delta \epsilon^2 \sin^3 \theta}{3\alpha^2} + b$ ,

using $+b$ for the constant of integration. Therefore, by variation of parameters, the general solution of the nonhomogeneous differential equation is

$u(\theta) = f_1(\theta) u_1(\theta) + f_2(\theta) u_2(\theta)$

$= \displaystyle \left( \frac{\cos \theta}{\alpha} + \frac{\delta}{3\alpha^2\epsilon} \left( 1 + \epsilon \cos \theta \right)^3 + a \right) \cos \theta$

$+ \displaystyle \left(\frac{\sin\theta}{\alpha} + \frac{\delta (1 + \epsilon^2) \sin\theta}{\alpha^2} + \frac{\delta \epsilon \theta}{\alpha^2} + \frac{\delta \epsilon \sin 2 \theta}{2\alpha^2} - \frac{\delta \epsilon^2 \sin^3 \theta}{3\alpha^2} + b \right) \sin \theta$

$= a\cos \theta + b \sin \theta + \displaystyle \frac{\cos^2 \theta + \sin^2 \theta}{\alpha} + \frac{\delta \cos \theta}{3\alpha^2\epsilon} + \frac{\delta \cos^2 \theta}{\alpha^2} + \frac{\delta \epsilon \cos^3 \theta}{\alpha^2}+ \frac{\delta\epsilon^2 \cos^4 \theta}{3\alpha^2}$

$+ \displaystyle \frac{\delta \sin^2 \theta}{\alpha^2} + \frac{\delta \epsilon^2 \sin^2\theta}{\alpha^2} + \frac{\delta \epsilon \theta \sin \theta}{\alpha^2} + \frac{\delta \epsilon \sin \theta\sin 2 \theta}{2\alpha^2} - \frac{\delta \epsilon^2 \sin^4 \theta}{3\alpha^2}$

$= a\cos \theta + \displaystyle \frac{\delta \cos \theta}{3\alpha^2\epsilon} + b \sin \theta + \displaystyle \frac{1}{\alpha} + \frac{\delta (\cos^2 \theta+\sin^2 \theta)}{\alpha^2} + \frac{\delta \epsilon \cos^3 \theta}{\alpha^2}+ \frac{\delta\epsilon^2 \cos^4 \theta}{3\alpha^2}$

$\displaystyle + \frac{\delta \epsilon^2 \sin^2\theta}{\alpha^2} + \frac{\delta \epsilon \theta \sin \theta}{\alpha^2} + \frac{\delta \epsilon \sin \theta\sin 2 \theta}{2\alpha^2} - \frac{\delta \epsilon^2 \sin^4 \theta}{3\alpha^2}$

$= a\cos \theta+ \displaystyle \frac{\delta \cos \theta}{3\alpha^2\epsilon} + b \sin \theta + \displaystyle \frac{1}{\alpha} + \frac{\delta}{\alpha^2} + \frac{\delta \epsilon \theta \sin \theta}{\alpha^2} + \frac{\delta \epsilon \cos^3 \theta}{\alpha^2} + \frac{\delta \epsilon \sin \theta\sin 2 \theta}{2\alpha^2}$

$\displaystyle + \frac{\delta \epsilon^2 \sin^2\theta}{\alpha^2} + \frac{\delta \epsilon^2 (\cos^4 \theta -\sin^4 \theta)}{3\alpha^2}$

$= a\cos \theta + \displaystyle \frac{\delta \cos \theta}{3\alpha^2\epsilon} + b \sin \theta + \displaystyle \frac{1}{\alpha} + \frac{\delta}{\alpha^2} + \frac{\delta \epsilon \theta \sin \theta}{\alpha^2} + \frac{\delta \epsilon \cos^3 \theta}{\alpha^2} + \frac{\delta \epsilon \sin \theta \cdot 2 \sin \theta \cos \theta}{2\alpha^2}$

$\displaystyle + \frac{\delta \epsilon^2 (1-\cos 2\theta)}{2\alpha^2} + \frac{\delta \epsilon^2 (\cos^2 \theta +\sin^2\theta)(\cos^2 \theta - \sin^2 \theta)}{3\alpha^2}$

$\displaystyle + \frac{\delta \epsilon^2}{2\alpha^2} -\frac{\delta \epsilon^2 \cos 2\theta}{2\alpha^2} + \frac{\delta \epsilon^2 \cos 2 \theta}{3\alpha^2}$

$= a\cos \theta + \displaystyle \frac{\delta \cos \theta}{3\alpha^2\epsilon}+ b \sin \theta + \displaystyle \frac{1}{\alpha} + \frac{\delta}{\alpha^2} + \frac{\delta \epsilon^2}{2\alpha^2} + \frac{\delta \epsilon \theta \sin \theta}{\alpha^2} -\frac{\delta \epsilon^2 \cos 2\theta}{6\alpha^2}$

$\displaystyle + \frac{\delta \epsilon \cos^3 \theta}{\alpha^2} + \frac{\delta \epsilon \sin^2 \theta\cos \theta}{\alpha^2}$

$= a\cos \theta + \displaystyle \frac{\delta \cos \theta}{3\alpha^2\epsilon}+ b \sin \theta + \displaystyle \frac{1}{\alpha} + \frac{\delta}{\alpha^2} + \frac{\delta \epsilon^2}{2\alpha^2} + \frac{\delta \epsilon \theta \sin \theta}{\alpha^2} -\frac{\delta \epsilon^2 \cos 2\theta}{6\alpha^2}$

$\displaystyle + \frac{\delta \epsilon \cos \theta (\cos^2 \theta + \sin^2 \theta)}{\alpha^2}$

$= a\cos \theta + \displaystyle \frac{\delta \cos \theta}{3\alpha^2\epsilon}+ b \sin \theta + \displaystyle \frac{1}{\alpha} + \frac{\delta}{\alpha^2} + \frac{\delta \epsilon^2}{2\alpha^2} + \frac{\delta \epsilon \theta \sin \theta}{\alpha^2} -\frac{\delta \epsilon^2 \cos 2\theta}{6\alpha^2} + \frac{\delta \epsilon \cos \theta}{\alpha^2}$

$= \displaystyle \left(a + \displaystyle \frac{\delta}{3\alpha^2\epsilon} + \frac{\delta \epsilon}{\alpha^2} \right) \cos \theta + b \sin \theta + \displaystyle \frac{1}{\alpha} + \frac{\delta}{\alpha^2} + \frac{\delta \epsilon^2}{2\alpha^2} + \frac{\delta \epsilon \theta \sin \theta}{\alpha^2} -\frac{\delta \epsilon^2 \cos 2\theta}{6\alpha^2}$

$= \displaystyle A \cos \theta + b \sin \theta + \displaystyle \frac{1}{\alpha} + \frac{\delta}{\alpha^2} + \frac{\delta \epsilon^2}{2\alpha^2} + \frac{\delta \epsilon \theta \sin \theta}{\alpha^2} -\frac{\delta \epsilon^2 \cos 2\theta}{6\alpha^2}$ ,

where $A$ is another arbitrary constant.

Next, we use the initial conditions to find the constants $A$ and $b$ . From the initial condition $u(0) = \displaystyle \frac{1}{P} = \frac{1+\epsilon}{\alpha}$ , we obtain

$u(0) = \displaystyle A \cos 0 + b \sin 0 + \displaystyle \frac{1}{\alpha} + \frac{\delta}{\alpha^2} + \frac{\delta \epsilon^2}{2\alpha^2} + \frac{\delta \epsilon \cdot 0 \cdot \sin 0}{\alpha^2} -\frac{\delta \epsilon^2 \cos 0}{6\alpha^2}$

$\displaystyle \frac{1+\epsilon}{\alpha} = A + \displaystyle \frac{1}{\alpha} + \frac{\delta}{\alpha^2} + \frac{\delta \epsilon^2}{2\alpha^2} -\frac{\delta \epsilon^2}{6\alpha^2}$

$\displaystyle \frac{\epsilon}{\alpha} = A + \displaystyle \frac{3\delta +\delta \epsilon^2}{3\alpha^2}$ ,

so that

$A = \displaystyle \frac{\epsilon}{\alpha} - \frac{\delta(3 + \epsilon^2)}{3\alpha^2}$ .

Next, we compute $u'(\theta)$ and use the initial condition $u'(0) = 0$ :

$u'(\theta) = \displaystyle -A \sin \theta + b \cos \theta + \frac{\delta \epsilon}{\alpha^2} (\sin \theta + \theta \cos \theta) + \frac{\delta \epsilon^2 \sin 2\theta}{3\alpha^2}$

$u'(0) = \displaystyle -A \sin 0 + b \cos 0 + \frac{\delta \epsilon}{\alpha^2} (\sin 0 + 0 \cos 0) + \frac{\delta \epsilon^2 \sin 0}{3\alpha^2}$

$0 = b$ .

Substituting these values for $A$ and $b$ , we finally arrive at the solution

$u(\theta) = \displaystyle \left(\frac{\epsilon}{\alpha} - \frac{\delta(3 + \epsilon^2)}{3\alpha^2} \right) \cos \theta + \displaystyle \frac{1}{\alpha} + \frac{\delta}{\alpha^2} + \frac{\delta \epsilon^2}{2\alpha^2} + \frac{\delta \epsilon \theta \sin \theta}{\alpha^2} -\frac{\delta \epsilon^2 \cos 2\theta}{6\alpha^2}$

$= \displaystyle \frac{1 + \epsilon \cos \theta}{\alpha} + \frac{\delta}{\alpha^2} + \frac{\delta \epsilon^2}{2\alpha^2} + \frac{ \delta\epsilon}{\alpha^2} \theta \sin \theta - \frac{ \delta \epsilon^2}{6\alpha^2} \cos 2\theta - \frac{\delta(3+\epsilon^2)}{3\alpha^2} \cos \theta$ .

Confirming Einstein’s Theory of General Relativity With Calculus, Part 6b: Checking Solution of New Differential Equation with Calculus

In the last post, we showed that if the motion of a planet around the Sun is expressed in polar coordinates $(r,\theta)$ , with the Sun at the origin, then under general relativity the motion of the planet follows the initial-value problem

$u''(\theta) + u(\theta) = \displaystyle \frac{1}{\alpha} + \delta \left( \frac{1 + \epsilon \cos \theta}{\alpha} \right)^2$ ,

$u(0) = \displaystyle \frac{1}{P}$ ,

$u'(0) = 0$ ,

I won’t sugar-coat it; the solution is a big mess:

$u(\theta) = \displaystyle \frac{1 + \epsilon \cos \theta}{\alpha} + \frac{\delta}{\alpha^2} + \frac{\delta \epsilon^2}{2\alpha^2} + \frac{\epsilon \delta}{\alpha^2} \theta \sin \theta - \frac{\epsilon^2 \delta}{6\alpha^2} \cos 2\theta - \frac{\delta(3+\epsilon^2)}{3\alpha^2} \cos \theta$ .

That said, it is an elementary, if complicated, exercise in calculus to confirm that this satisfies all three equations above. We’ll start with the second one:

$u(0) = \displaystyle \frac{1 + \epsilon \cos 0}{\alpha} + \frac{\delta}{\alpha^2} + \frac{\delta \epsilon^2}{2\alpha^2} + \frac{\epsilon \delta}{\alpha^2} \cdot 0 \sin 0 - \frac{\epsilon^2 \delta}{6\alpha^2} \cos 0 - \frac{\delta(3+\epsilon^2)}{3\alpha^2} \cos 0$

$= \displaystyle \frac{1 + \epsilon}{\alpha} + \frac{\delta}{\alpha^2} + \frac{\delta \epsilon^2}{2\alpha^2} - \frac{\epsilon^2 \delta}{6\alpha^2} - \frac{\delta(3+\epsilon^2)}{3\alpha^2}$

$= \displaystyle \frac{ 6\alpha(1+\epsilon) + 6\delta + 3 \delta \epsilon^2 - \delta \epsilon^2 - 2\delta (3+\epsilon^2)}{6\alpha^2}$

$= \displaystyle \frac{ 6\alpha(1+\epsilon) + 6\delta + 3 \delta \epsilon^2 - \delta \epsilon^2 - 6\delta - 2\delta \epsilon^2}{6\alpha^2}$

$= \displaystyle \frac{ 6\alpha(1+\epsilon)}{6\alpha^2}$

$= \displaystyle \frac{ 1+\epsilon}{\alpha}$

$= \displaystyle \frac{1}{P}$ ,

where in the last step we used the equation $P = \displaystyle \frac{\alpha}{1 + \epsilon}$ that was obtained earlier in this series.

Next, to check the initial condition $u'(0) = 0$ , we differentiate:

$u'(\theta) = \displaystyle -\frac{\epsilon \sin \theta}{\alpha} + \frac{\epsilon \delta}{\alpha^2} (\sin \theta + \theta \cos \theta) + \frac{\epsilon^2 \delta}{3\alpha^2} \sin 2\theta + \frac{\delta(3+\epsilon^2)}{3\alpha^2} \sin\theta$

$u'(0) = \displaystyle -\frac{\epsilon \sin 0}{\alpha} + \frac{\epsilon \delta}{\alpha^2} (\sin 0 + 0 \cdot \cos 0) + \frac{\epsilon^2 \delta}{3\alpha^2} \sin 0 + \frac{\delta(3+\epsilon^2)}{3\alpha^2} \sin 0 = 0$ .

Finally, to check the differential equation itself, we compute the second derivative:

$u''(\theta) = \displaystyle -\frac{\epsilon \cos \theta}{\alpha} + \frac{\epsilon \delta}{\alpha^2} (2 \cos \theta - \theta \sin \theta) + \frac{2\epsilon^2 \delta}{3\alpha^2} \cos 2\theta + \frac{\delta(3+\epsilon^2)}{3\alpha^2} \cos\theta$ .

Adding $u''(\theta)$ and $u(\theta)$ , we find

$u''(\theta) + u(\theta) = \displaystyle \frac{1 + \epsilon \cos \theta - \epsilon \cos \theta}{\alpha} + \frac{\delta}{\alpha^2} + \frac{\delta \epsilon^2}{2\alpha^2} + \frac{\epsilon \delta}{\alpha^2} (\theta \sin \theta + 2 \cos \theta - \theta \sin \theta)$

$\displaystyle - \frac{\epsilon^2 \delta}{6\alpha^2} \cos 2\theta + \frac{2\epsilon^2 \delta}{3\alpha^2} \cos 2\theta - \frac{\delta(3+\epsilon^2)}{3\alpha^2} \cos \theta + \frac{\delta(3+\epsilon^2)}{3\alpha^2} \cos\theta$ ,

which simplifies considerably:

$u''(\theta) + u(\theta) = \displaystyle \frac{1}{\alpha} + \frac{\delta}{\alpha^2} + \frac{\delta \epsilon^2}{2\alpha^2} + \frac{2 \epsilon \delta}{\alpha^2} \cos \theta + \frac{\epsilon^2 \delta}{2\alpha^2} \cos 2\theta$

$= \displaystyle \frac{1}{\alpha} + \frac{\delta}{\alpha^2} \left( 1 + \frac{\epsilon^2}{2} + 2 \epsilon \cos \theta + \frac{\epsilon^2}{2} \cos 2\theta \right)$

$= \displaystyle \frac{1}{\alpha} + \frac{\delta}{\alpha^2} \left( 1 + 2 \epsilon \cos \theta + \epsilon^2 \frac{1+\cos 2\theta}{2} \right)$

$= \displaystyle \frac{1}{\alpha} + \frac{\delta}{\alpha^2} \left( 1 + 2 \epsilon \cos \theta + \epsilon^2 \cos^2 \theta \right)$

$= \displaystyle \frac{1}{\alpha} + \frac{\delta}{\alpha^2} (1 + \epsilon \cos \theta)^2$

$= \displaystyle \frac{1}{\alpha} + \delta \left( \frac{1+\epsilon \cos \theta}{\alpha} \right)^2$ ,

where we used the power-reduction trigonometric identity

$\cos^2 \theta = \displaystyle \frac{1 + \cos 2\theta}{2}$

on the second-to-last step.

While we have verified the proposed solution of the initial-value problem, and the steps for doing so lie completely within the grasp of a good calculus student, I’ll be the first to say that this solution is somewhat unsatisfying: the solution appeared seemingly out of thin air, and we just checked to see if this mysterious solution actually works. In the next few posts, I’ll discuss how this solution can be derived using standard techniques from first-semester differential equations.

Confirming Einstein’s Theory of General Relativity With Calculus, Part 6a: New Differential Equation

We previously showed that if the motion of a planet around the Sun is expressed in polar coordinates $(r,\theta)$ , with the Sun at the origin, then under Newtonian mechanics (i.e., without general relativity) the motion of the planet follows the differential equation

$u_0''(\theta) + u_0(\theta) = \displaystyle \frac{Gm^2 M}{\ell^2}$ ,

where $u_0 = 1/r$ , $G$ is the gravitational constant of the universe, $m$ is the mass of the planet, $M$ is the mass of the Sun, and $\ell$ is the constant angular momentum of the planet. For simplicity, we wrote $\displaystyle \frac{1}{\alpha} = \displaystyle \frac{Gm^2 M}{\ell^2}$ .

We will also impose the initial condition that the planet is at perihelion (i.e., is closest to the sun), at a distance of $P$ , when $\theta = 0$ . This means that $u$ obtains its maximum value of $1/P$ when $\theta = 0$ . This leads to the two initial conditions

$u_0(0) = \displaystyle \frac{1}{P} \qquad \hbox{and} \qquad u_0'(0) = 0$ ;

the second equation arises since $u$ has a local extremum at $\theta = 0$ .

In previous posts, we showed that the solution of this initial-value problem is

$u_0(\theta) = \displaystyle \frac{1 + \epsilon \cos \theta}{\alpha}$ ,

where $\epsilon = \displaystyle \frac{\alpha - P}{P}$ . Since $r = 1/u_0$ , we see that the planet’s orbit satisfies

$r = \displaystyle \frac{\alpha}{1 + \epsilon \cos \theta}$ ,

so that, as shown earlier in this series, the orbit is an ellipse with eccentricity $\epsilon$ .

At long last, we’re now ready to see what happens under general relativity. According to the theory of general relativity, the governing differential equation for the orbit of a planet should be changed ever so slightly to

$u''(\theta) + u(\theta) = \displaystyle \frac{Gm^2 M}{\ell^2} + \frac{3GM}{c^2} [u(\theta)]^2$ ,

where a second term was added to the right-hand side. (I will make no attempt here to justify the physics for this second term.) In this second term, $c$ represents the speed of light. Since $c$ is very large, this second term is very, very small. Using $\displaystyle \frac{1}{\alpha} = \displaystyle \frac{Gm^2 M}{\ell^2}$ as before and defining $\delta = \displaystyle \frac{3GM}{c^2}$ , this may be simplified to

$u''(\theta) + u(\theta) = \displaystyle \frac{1}{\alpha} + \delta [u(\theta)]^2$ ,

Finding an exact solution of this new differential equation is hopeless. While the previous differential equation was linear with constant coefficients, this new differential is decidedly nonlinear because of the new term $[u(\theta)]^2$ .

So, instead of attempting to find an exact solution, we will use the method of successive approximations, mentioned earlier in this series, to find an approximate solution that is very, very close to the exact solution. Since $\delta$ is very, very small, the solution of

$u''(\theta) + u(\theta) = \displaystyle \frac{1}{\alpha}$

will be approximately equal to the solution of the “real” differential equation under general relativity. We have already shown that the solution of this simpler differential equation is

$u_0(\theta) = \displaystyle \frac{1 + \epsilon \cos \theta}{\alpha}$ .

Therefore, the method of successive approximations suggests that a better approximation to the true orbit will be the solution of

$u''(\theta) + u(\theta) = \displaystyle \frac{1}{\alpha} + \delta [u_0(\theta)]^2$ ,

$u''(\theta) + u(\theta) = \displaystyle \frac{1}{\alpha} + \delta \left( \frac{1 + \epsilon \cos \theta}{\alpha} \right)^2$ .

While significantly messier than the differential equation under Newtonian mechanics, this is now a linear differential equation that can be solved using standard techniques from differential equations. In the next few posts, we will solve this new differential equation, thus finding the predicted orbit of a planet under general relativity.

Overlapping Circles

Posted in honor of the upcoming North American solar eclipse.

Source: https://xkcd.com/2769

Confirming Einstein’s Theory of General Relativity With Calculus, Part 5d: Deriving Orbits under Newtonian Mechanics Using Variation of Parameters

We previously showed that if the motion of a planet around the Sun is expressed in polar coordinates $(r,theta)$ , with the Sun at the origin, then under Newtonian mechanics (i.e., without general relativity) the motion of the planet follows the differential equation

$u''(\theta) + u(\theta) = \displaystyle \frac{1}{\alpha}$ ,

where $u = 1/r$ and $\alpha$ is a certain constant. We will also impose the initial condition that the planet is at perihelion (i.e., is closest to the sun), at a distance of $P$ , when $\theta = 0$ . This means that $u$ obtains its maximum value of $1/P$ when $\theta = 0$ . This leads to the two initial conditions

$u(0) = \displaystyle \frac{1}{P} \qquad \hbox{and} \qquad u'(0) = 0$ ;

the second equation arises since $u$ has a local extremum at $\theta = 0$ .

$u''(\theta) + u(\theta) = 0$ .

According to the method of variation of parameters, the general solution of the original nonhomogeneous differential equation

$u''(\theta) + u(\theta) = g(\theta)$

$u(\theta) = f_1(\theta) u_1(\theta) + f_2(\theta) u_2(\theta)$ ,

where

$f_1(\theta) = -\displaystyle \int \frac{u_2(\theta) g(\theta)}{W(\theta)} d\theta$ ,

$f_2(\theta) = \displaystyle \int \frac{u_1(\theta) g(\theta)}{W(\theta)} d\theta$ ,

and $W(\theta)$ is the Wronskian of $u_1(\theta)$ and $u_2(\theta)$ , defined by the determinant

$W(\theta) = \displaystyle \begin{vmatrix} u_1(\theta) & u_2(\theta) \\ u_1'(\theta) & u_2'(\theta) \end{vmatrix} = u_1(\theta) u_2'(\theta) - u_1'(\theta) u_2(\theta)$ .

Well, that’s a mouthful.

Fortunately, for the example at hand, these computations are pretty easy. First, since $u_1(\theta) = \cos \theta$ and $u_2(\theta) = \sin \theta$ , we have

$W(\theta) = (\cos \theta)(\cos \theta) - (\sin \theta)(-\sin \theta) = \cos^2 \theta + \sin^2 \theta = 1$

from the usual Pythagorean trigonometric identity. Therefore, the denominators in the integrals for $f_1(\theta)$ and $f_2(\theta)$ essentially disappear.

Since $g(\theta) = \displaystyle \frac{1}{\alpha}$ , the integrals for $f_1(\theta)$ and $f_2(\theta)$ are straightforward to compute:

$f_1(\theta) = -\displaystyle \int u_2(\theta) \frac{1}{\alpha} d\theta = -\displaystyle \frac{1}{\alpha} \int \sin \theta \, d\theta = \displaystyle \frac{1}{\alpha}\cos \theta + a$ ,

where we use $+a$ for the constant of integration instead of the usual $+C$ . Second,

$f_2(\theta) = \displaystyle \int u_1(\theta) \frac{1}{\alpha} d\theta = \displaystyle \frac{1}{\alpha} \int \cos \theta \, d\theta = \displaystyle \frac{1}{\alpha}\sin \theta + b$ ,

using $+b$ for the constant of integration. Therefore, by variation of parameters, the general solution of the nonhomogeneous differential equation is

$u(\theta) = f_1(\theta) u_1(\theta) + f_2(\theta) u_2(\theta)$

$= \left( \displaystyle \frac{1}{\alpha}\cos \theta + a \right) \cos \theta + \left( \displaystyle \frac{1}{\alpha}\sin\theta + b \right) \sin \theta$

$= a \cos \theta + b\sin \theta + \displaystyle \frac{\cos^2 \theta + \sin^2 \theta}{\alpha}$

$= a \cos \theta + b \sin \theta + \displaystyle \frac{1}{\alpha}$ .

Unsurprisingly, this matches the answer in the previous post that was found by the method of undetermined coefficients.

For the sake of completeness, I repeat the argument used in the previous two posts to determine $a$ and $b$ . This is require using the initial conditions $u(0) = \displaystyle \frac{1}{P}$ and $u'(0) = 0$ . From the first initial condition,

$u(0) = a \cos 0 + b \sin 0 + \displaystyle \frac{1}{\alpha}$

$\displaystyle \frac{1}{P} = a + \frac{1}{\alpha}$

$\displaystyle \frac{1}{P} - \frac{1}{\alpha} = a$

$\displaystyle \frac{\alpha - P}{\alpha P} = a$

From the second initial condition,

$u'(\theta) = -a \sin \theta + b \cos \theta$

$u'(0) = -a \sin 0 + b \cos 0$

$0 = b$ .

From these two constants, we obtain

$u(\theta) = \displaystyle \frac{\alpha - P}{\alpha P} \cos \theta + 0 \sin \theta + \displaystyle \frac{1}{\alpha}$

$= \displaystyle \frac{1}{\alpha} \left( 1 + \frac{\alpha-P}{P} \cos \theta \right)$

$= \displaystyle \frac{1 + \epsilon \cos \theta}{\alpha}$ ,

where $\epsilon = \displaystyle \frac{\alpha - P}{P}$ .

Finally, since $r = 1/u$ , we see that the planet’s orbit satisfies

$r = \displaystyle \frac{\alpha}{1 + \epsilon \cos \theta}$ ,

so that, as shown earlier in this series, the orbit is an ellipse with eccentricity $\epsilon$ .

Confirming Einstein’s Theory of General Relativity With Calculus, Part 5c: Deriving Orbits under Newtonian Mechanics with the Method of Undetermined Coefficients

$u''(\theta) + u(\theta) = \displaystyle \frac{1}{\alpha}$ ,

$u(0) = \displaystyle \frac{1}{P} \qquad \hbox{and} \qquad u'(0) = 0$ ;

the second equation arises since $u$ has a local extremum at $\theta = 0$ .

We now take the perspective of a student who is taking a first-semester course in differential equations. There are two standard techniques for solving a second-order non-homogeneous differential equations with constant coefficients. One of these is the method of undetermined coefficients. First, we solve the associated homogeneous differential equation

$u''(\theta) + u(\theta) = 0$ .

The characteristic equation of this differential equation is $r^2 + 1 = 0$ , which clearly has the two imaginary roots $r = \pm i$ . Therefore, the solution of the associated homogeneous equation is $u_h(\theta) = a \cos \theta + b \sin \theta$ .

Next, we find a particular solution to the original differential equation. Since the right-hand side is a constant and $r=0$ is not a solution of the characteristic equation, this leads to trying something of the form $U(\theta) = A$ as a solution, where $A$ is a soon-to-be-determined constant. (Guessing this form of $U(\theta)$ is a standard technique from differential equations; later in this series, I’ll give some justification for this guess.)

Clearly, $U'(\theta) = 0$ and $U''(\theta) = 0$ . Substituting, we find

$U''(\theta) + U(\theta) = \displaystyle \frac{1}{\alpha}$

$0 + A = \displaystyle \frac{1}{\alpha}$

$A = \displaystyle \frac{1}{\alpha}$

Therefore, $U(\theta) = \displaystyle \frac{1}{\alpha}$ is a particular solution of the nonhomogeneous differential equation.

Next, the general solution of the nonhomogeneous differential equation is found by adding the general solution to the associated homogeneous differential equation and the particular solution:

$u(\theta) = u_h(\theta) + U(\theta) = a \cos \theta + b\sin \theta + \displaystyle \frac{1}{\alpha}$ .

Finally, to determine $a$ and $b$ , we use the initial conditions $u(0) = \displaystyle \frac{1}{P}$ and $u'(0) = 0$ . From the first initial condition,

$u(0) = a \cos 0 + b \sin 0 + \displaystyle \frac{1}{\alpha}$

$\displaystyle \frac{1}{P} = a + \frac{1}{\alpha}$

$\displaystyle \frac{1}{P} - \frac{1}{\alpha} = a$

$\displaystyle \frac{\alpha - P}{\alpha P} = a$

From the second initial condition,

$u'(\theta) = -a \sin \theta + b \cos \theta$

$u'(0) = -a \sin 0 + b \cos 0$

$0 = b$ .

From these two constants, we obtain

$u(\theta) = \displaystyle \frac{\alpha - P}{\alpha P} \cos \theta + 0 \sin \theta + \displaystyle \frac{1}{\alpha}$

$= \displaystyle \frac{1}{\alpha} \left( 1 + \frac{\alpha-P}{P} \cos \theta \right)$

$= \displaystyle \frac{1 + \epsilon \cos \theta}{\alpha}$ ,

where $\epsilon = \displaystyle \frac{\alpha - P}{P}$ .

Finally, since $r = 1/u$ , we see that the planet’s orbit satisfies

$r = \displaystyle \frac{\alpha}{1 + \epsilon \cos \theta}$ ,

so that, as shown earlier in this series, the orbit is an ellipse with eccentricity $\epsilon$ .

The reader will notice that this solution is pretty much a carbon-copy of the previous post. The difference is that calculus students wouldn’t necessarily be able to independently generate the solutions of the associated homogeneous differential equation and a particular solution of the nonhomogeneous differential equation, and so the line of questioning is designed to steer students toward the answer.

Confirming Einstein’s Theory of General Relativity With Calculus, Part 5b: Deriving Orbits under Newtonian Mechanics with Calculus

$u''(\theta) + u(\theta) = \displaystyle \frac{1}{\alpha}$ ,

$u(0) = \displaystyle \frac{1}{P} \qquad \hbox{and} \qquad u'(0) = 0$ ;

the second equation arises since $u$ has a local extremum at $\theta = 0$ .

In the previous post, we confirmed that

$u(\theta) = \displaystyle \frac{1 + \epsilon \cos \theta}{\alpha}$

solved this initial-value problem. However, the solution was unsatisfying because it gave no indication of where this guess might have come from. In this post, I suggest a series of questions that good calculus students could be asked that would hopefully lead them quite naturally to this solution.

Step 1. Let’s make the differential equation simpler, for now, by replacing the right-hand side with 0:

$u''(\theta) + u(\theta) = 0$ ,

or

$u''(\theta) = -u(\theta)$ .

Can you think of a function or two that, when you differentiate twice, you get the original function back, except with a minus sign in front?

Answer to Step 1. With a little thought, hopefully students can come up with the standard answers of $u(\theta) = \cos \theta$ and $u(\theta) = \sin \theta$ .

Step 2. Using these two answers, can you think of a third function that works?

Answer to Step 2. This is usually the step that students struggle with the most, as they usually try to think of something completely different that works. This won’t work, but that’s OK… we all learn from our failures. If they can’t figure out, I’ll give a big hint: “Try multiplying one of these two answers by something.” In time, they’ll see that answers like $u(\theta) = 2\cos \theta$ and $u(\theta) = 3\sin \theta$ work. Once that conceptual barrier is broken, they’ll usually produce the solutions $u(\theta) = a \cos \theta$ and $u(\theta) = b \sin \theta$ .

Step 3. Using these two answers, can you think of anything else that works?

Answer to Step 3. Again, students might struggle as they imagine something else that works. If this goes on for too long, I’ll give a big hint: “Try combining them.” Eventually, we hopefully get to the point that they’ll see that the linear combination $u(\theta) = a \cos \theta + b \sin \theta$ also solves the associated homogeneous differential equation.

Step 4. Let’s now switch back to the original differential equation $u''(\theta) + u(\theta) = \displaystyle \frac{1}{\alpha}$ . Let’s start simple: $u''(\theta) + u(\theta) = 5$ . Can you think of an easy function that’s a solution?

Answer to Step 4. This might take some experimentation, and students will probably try unnecessarily complicated guesses first. If this goes on for too long, I’ll give a big hint: “Try a constant.” Eventually, they hopefully determine that if $u(\theta) = 5$ is a constant function, then clearly $u'(\theta) = 0$ and $u''(\theta) = 0$ , so that $u''(\theta) + u(\theta) = 5$ .

Step 5. Let’s return to $u''(\theta) + u(\theta) = \displaystyle \frac{1}{\alpha}$ . Any guesses on an answer to this one?

Answer to Step 5. Hopefully, students quickly realize that the constant function $u(\theta) = \displaystyle \frac{1}{\alpha}$ works.

Step 6. Let’s review. We’ve shown that anything of the form $u(\theta) = a\cos \theta + b \sin \theta$ is a solution of $u''(\theta) + u(\theta) = 0$ . We’ve also shown that $u(\theta) = \displaystyle\frac{1}{\alpha}$ is a solution of $u''(\theta) + u(\theta) = \displaystyle \frac{1}{\alpha}$ . Can you think use these two answers to find something else that works?

Answer to Step 6. Hopefully, with the experience learned from Step 3, students will guess that $u(\theta) = a\cos \theta + b\sin \theta + \displaystyle \frac{1}{\alpha}$ will work.

Step 7. OK, that solves the differential equation. Any thoughts on how to find the values of $a$ and $b$ so that $u(0) = \displaystyle \frac{1}{P}$ and $u'(0) = 0$ ?

Answer to Step 7. Hopefully, students will see that we should just plug into $u(\theta)$ :

$u(0) = a \cos 0 + b \sin 0 + \displaystyle \frac{1}{\alpha}$

$\displaystyle \frac{1}{P} = a + \frac{1}{\alpha}$

$\displaystyle \frac{1}{P} - \frac{1}{\alpha} = a$

$\displaystyle \frac{\alpha - P}{\alpha P} = a$

To find $b$ , we first find $u'(\theta)$ and then substitute $\theta = 0$ :

$u'(\theta) = -a \sin \theta + b \cos \theta$

$u'(0) = -a \sin 0 + b \cos 0$

$0 = b$ .

From these two constants, we obtain

$u(\theta) = \displaystyle \frac{\alpha - P}{\alpha P} \cos \theta + 0 \sin \theta + \displaystyle \frac{1}{\alpha}$

$= \displaystyle \frac{1}{\alpha} \left( 1 + \frac{\alpha-P}{P} \cos \theta \right)$

$= \displaystyle \frac{1 + \epsilon \cos \theta}{\alpha}$ ,

where $\epsilon = \displaystyle \frac{\alpha - P}{P}$ .

Finally, since $r = 1/u$ , we see that the planet’s orbit satisfies

$r = \displaystyle \frac{\alpha}{1 + \epsilon \cos \theta}$ ,

so that, as shown earlier in this series, the orbit is an ellipse with eccentricity $\epsilon$ .

Confirming Einstein’s Theory of General Relativity With Calculus, Part 5a: Confirming Orbits under Newtonian Mechanics with Calculus

$u''(\theta) + u(\theta) = \displaystyle \frac{1}{\alpha}$ ,

$u(0) = \displaystyle \frac{1}{P} \qquad \hbox{and} \qquad u'(0) = 0$ ;

the second equation arises since $u$ has a local extremum at $\theta = 0$ .

In the next few posts, we’ll discuss the solution of this initial-value problem. Today’s post would be appropriate for calculus students, which is confirming that

$u(\theta) = \displaystyle \frac{1 + \epsilon \cos \theta}{\alpha}$

solves this initial-value problem, where $\epsilon = \displaystyle \frac{\alpha-P}{P}$ . Since $r$ is the reciprocal of $u$ , we infer that

$r = \displaystyle \frac{\alpha}{1 + \epsilon \cos \theta}$ .

As we’ve already seen in this series, this means that the orbit of the planet is a conic section — either a circle, ellipse, parabola, or hyperbola. Since the orbit of a planet is stable and $\epsilon = 0$ is extremely unlikely, this means that the planet orbits the Sun in an ellipse, with the Sun at one focus of the ellipse.

So, for a calculus student to verify that planets move in ellipses, one must check that

$u(\theta) = \displaystyle \frac{1 + \epsilon \cos \theta}{\alpha}$

is a solution of the initial-value problem

$u''(\theta) + u(\theta) = \displaystyle \frac{1}{\alpha}$ ,

$u(0) = \displaystyle \frac{1}{P}$ ,

$u'(0) = 0$ .

The second line is easy to check:

$u(0) = \displaystyle \frac{1 + \epsilon \cos 0}{\alpha}$

$= \displaystyle \frac{1 + \epsilon}{\alpha}$

$= \displaystyle \frac{1 + \displaystyle \frac{\alpha-P}{P}}{\alpha}$

$= \displaystyle \frac{1}{\alpha} \frac{P + \alpha - P}{P}$

$= \displaystyle \frac{1}{\alpha} \frac{\alpha}{P}$

$= \displaystyle \frac{1}{P}$ .

The third line is also easy to check:

$u'(\theta) = \displaystyle \frac{-\epsilon \sin \theta}{\alpha}$

$u'(0) = \displaystyle \frac{-\epsilon \sin 0}{\alpha} = 0$ .

To check the first line, we first find $u''(\theta)$ :

$u''(\theta) = \displaystyle \frac{-\epsilon \cos \theta}{\alpha}$ ,

so that

$u''(\theta) + u(\theta) = \displaystyle \frac{-\epsilon \cos \theta}{\alpha} + \frac{1 + \epsilon \cos \theta}{\alpha} = \frac{1}{\alpha}$ ,

thus confirming that $u(\theta) = \displaystyle \frac{1 + \epsilon \cos \theta}{\alpha}$ solves the initial-value problem.

While the above calculations are well within the grasp of a good Calculus I student, I’ll be the first to admit that this solution is less than satisfying. We just mysteriously proposed a solution, seemingly out of thin air, and confirmed that it worked. In the next post, I’ll proposed a way that calculus students can be led to guess this solution. Then, we talk about finding the solution of this nonhomogeneous initial-value problem using standard techniques from differential equations.

Confirming Einstein’s Theory of General Relativity With Calculus, Part 4c: Newton’s Second Law and Newton’s Law of Gravitation

In this post, following from the previous two posts, we will show that if the motion of a planet around the Sun is expressed in polar coordinates $(r,\theta)$ , with the Sun at the origin, then under Newtonian mechanics (i.e., without general relativity) the motion of the planet follows the differential equation

$u''(\theta) + u(\theta) = \displaystyle \frac{1}{\alpha}$ ,

where $u = 1/r$ and $\alpha$ is a certain constant. Deriving this governing differential equation will require some principles from physics. If you’d rather skip the physics and get to the mathematics, we’ll get to solving this differential equations in the next post.

From Newton’s second law, the gravitational force on the planet as it orbits the Sun satisfies

${\bf F} = m{\bf a}$ ,

where the force ${\bf F}$ and the acceleration ${\bf a}$ are vectors. When written in polar coordinates, this becomes

${\bf F} = m \displaystyle \left[ \frac{d^2r}{dt^2} - r \left( \frac{d\theta}{dt} \right)^2 \right] {\bf u}_r + m \left(r \frac{d^2 \theta}{d t^2} + 2 \frac{dr}{dt} \frac{d\theta}{dt} \right) {\bf u}_\theta$ ,

where ${\bf u}_r$ is a unit vector pointing away from the origin and ${\bf u}_\theta$ is a unit vector perpendicular to ${\bf u}_r$ that points in the direction of increasing $\theta$ .

Furthermore, from Newton’s Law of Gravitation, if the Sun is located at the origin, then the gravitational force on the planet is

${\bf F} = \displaystyle -\frac{GMm}{r^2} {\bf u}_r$ ,

where $M$ is the mass of the sun, $m$ is the mass of the planet, and $G$ is the gravitational constant of the universe (which is a constant, no matter what Q from Star Trek: The Next Generation says).

Since these are the same force, the ${\bf u}_r$ components must be the same. (Also, the ${\bf u}_\theta$ component must be zero, but we won’t need to use that fact.) Therefore,

$m \displaystyle \left[ \frac{d^2r}{dt^2} - r \left( \frac{d\theta}{dt} \right)^2 \right] = \displaystyle -\frac{GMm}{r^2}$ ,

$\displaystyle \frac{d^2r}{dt^2} - r \left( \frac{d\theta}{dt} \right)^2 = \displaystyle -\frac{GM}{r^2}$ .

In a previous post, we showed that

$\displaystyle \frac{d\theta}{dt} = \frac{\ell}{mr^2}$ ,

where $\ell$ is a constant, and

$\displaystyle \frac{d^2r}{dt^2} = - \frac{\ell^2}{m^2 r^2} \frac{d^2}{d\theta^2} \left( \frac{1}{r} \right)$ .

Substituting, we find

$\displaystyle - \frac{\ell^2}{m^2 r^2} \frac{d^2}{d\theta^2} \left( \frac{1}{r} \right) - r \left( \frac{\ell}{mr^2} \right)^2 = \displaystyle -\frac{GM}{r^2}$

$\displaystyle \frac{\ell^2}{m^2 r^2} \frac{d^2}{d\theta^2} \left( \frac{1}{r} \right) + r \left( \frac{\ell^2}{m^2 r^4} \right) = \displaystyle \frac{GM}{r^2}$

$\displaystyle \frac{\ell^2}{m^2 r^2} \frac{d^2}{d\theta^2} \left( \frac{1}{r} \right) +\frac{\ell^2}{m^2 r^3} = \displaystyle \frac{GM}{r^2}$

$\displaystyle \frac{\ell^2}{m^2 r^2} \left[ \frac{d^2}{d\theta^2} \left( \frac{1}{r} \right) + \frac{1}{r} \right] = \displaystyle \frac{GM}{r^2}$

$\displaystyle \frac{d^2}{d\theta^2} \left( \frac{1}{r} \right) + \frac{1}{r} = \displaystyle \frac{GM}{r^2} \cdot \frac{m^2 r^2}{\ell^2}$

$\displaystyle \frac{d^2}{d\theta^2} \left( \frac{1}{r} \right) + \frac{1}{r} = \displaystyle \frac{GMm^2}{\ell^2}$ .

So, substituting $u = 1/r$ and $\alpha = \displaystyle \frac{\ell^2}{GMm^2}$ , we finally obtain the governing equation

$\displaystyle \frac{d^2 u}{d\theta^2} + u = \displaystyle \frac{1}{\alpha}$ .

This is the governing differential equation of planetary motion under Newtonian mechanics. For now, it’s not obvious why we chose $\displaystyle \frac{1}{\alpha}$ as the constant on the right-hand side instead of just $\alpha$ , but the reason for this choice will become apparent in future posts.

In the next few posts, we use differential equations (or, if you’d prefer, just calculus) to show that Newtonian mechanics predicts that planets orbit the Sun in ellipses.

Confirming Einstein’s Theory of General Relativity With Calculus, Part 4b: Acceleration in Polar Coordinates

In this part of the series, we will show that if the motion of a planet around the Sun is expressed in polar coordinates $(r,\theta)$ , with the Sun at the origin, then under Newtonian mechanics (i.e., without general relativity) the motion of the planet follows the differential equation

$u''(\theta) + u(\theta) = \displaystyle \frac{1}{\alpha}$ ,

Part of the derivation of this governing differential equation will involve Newton’s Second Law

${\bf F} = m {\bf a}$ ,

where $m$ is the mass of the planet and the force ${\bf F}$ and the acceleration $a$ are vectors. In usual rectangular coordinates, the acceleration vector would be expressed as

${\bf a} = x''(t) {\bf i} + y''(t) {\bf j}$ ,

where the components of the acceleration in the $x-$ and $y-$ directors are $x''(t)$ and $y''(t)$ , and the unit vectors ${\bf i}$ and ${\bf j}$ are perpendicular, pointing in the positive $x$ and positive $y$ directions.

Unfortunately, our problem involves polar coordinates, and rewriting the acceleration vector in polar coordinates, instead of rectangular coordinates, is going to take some work.

Suppose that the position of the planet is $(r,\theta)$ in polar coordinates, so that the position in rectangular coordinates is ${\bf r} = (r\cos \theta, r \sin \theta)$ . This may be rewritten as

${\bf r} = r \cos \theta {\bf i} + r \sin \theta {\bf j} = r ( \cos \theta {\bf i} + \sin \theta {\bf j}) = r {\bf u}_r$ ,

where

${\bf u}_r = \cos \theta {\bf i} + \sin \theta {\bf j}$

is a unit vector that points away from the origin. We see that this is a unit vector since

$\parallel {\bf u}_r \parallel = {\bf u}_r \cdot {\bf u}_r = \cos^2 \theta + \sin^2 \theta =1$ .

We also define

${\bf u}_\theta = -\sin \theta {\bf i} + \cos \theta {\bf j}$

to be a unit vector that is perpendicular to ${\bf u}_r$ ; it turns out that ${\bf u}_\theta$ points in the direction of increasing $\theta$ . To see that ${\bf u}_r$ and ${\bf u}_\theta$ are perpendicular, we observe

${\bf u}_r \cdot {\bf u}_\theta = -\sin \theta \cos \theta + \sin \theta \cos \theta = 0$ .

Computing the velocity and acceleration vectors in polar coordinates will have a twist that’s not experienced with rectangular coordinates since both ${\bf u}_r$ and ${\bf u}_\theta$ are functions of $\theta$ . Indeed, we have

$\displaystyle \frac{d{\bf u}_r}{d\theta} = \frac{d \cos \theta}{d\theta} {\bf i} + \frac{d\sin \theta}{d\theta} {\bf j} = -\sin \theta {\bf i} + \cos \theta {\bf j} = {\bf u}_\theta$ .

Furthermore,

$\displaystyle \frac{d{\bf u}_\theta}{d\theta} = -\frac{d \sin \theta}{d\theta} {\bf i} + \frac{d\cos \theta}{d\theta} {\bf j} = -\cos \theta {\bf i} - \sin \theta {\bf j} = -{\bf u}_r$ .

These two equations will be needed in the derivation below.

We are now in position to express the velocity and acceleration of the orbiting planet in polar coordinates. Clearly, the position of the planet is $r {\bf u}_r$ , or a distance $r$ from the origin in the direction of ${\bf u}_r$ . Therefore, by the Product Rule, the velocity of the planet is

${\bf v} = \displaystyle \frac{d}{dt} (r {\bf u}_r) = \displaystyle \frac{dr}{dt} {\bf u}_r + r \frac{d {\bf u}_r}{dt}$

We now apply the Chain Rule to the second term:

${\bf v} = \displaystyle \frac{dr}{dt} {\bf u}_r + r \frac{d {\bf u}_r}{d\theta} \frac{d\theta}{dt}$

$= \displaystyle \frac{dr}{dt} {\bf u}_r + r \frac{d\theta}{dt} {\bf u}_\theta$ .

Differentiating a second time with respect to time, and again using the Chain Rule, we find

${\bf a} = \displaystyle \frac{d {\bf v}}{dt} = \displaystyle \frac{d^2r}{dt^2} {\bf u}_r + \frac{dr}{dt} \frac{d{\bf u}_r}{dt} + \frac{dr}{dt} \frac{d\theta}{dt} {\bf u}_\theta + r \frac{d^2\theta}{dt^2} {\bf u}_\theta + r \frac{d\theta}{dt} \frac{d{\bf u}_\theta}{dt}$

$= \displaystyle \frac{d^2r}{dt^2} {\bf u}_r + \frac{dr}{dt} \frac{d{\bf u}_r}{d\theta} \frac{d\theta}{dt} + \frac{dr}{dt} \frac{d\theta}{dt} {\bf u}_\theta + r \frac{d^2\theta}{dt^2} {\bf u}_\theta + r \frac{d\theta}{dt} \frac{d{\bf u}_\theta}{d\theta} \frac{d\theta}{dt}$

$= \displaystyle \frac{d^2r}{dt^2} {\bf u}_r + \frac{dr}{dt} \frac{d\theta}{dt} {\bf u}_\theta + \frac{dr}{dt} \frac{d\theta}{dt} {\bf u}_\theta + r \frac{d^2\theta}{dt^2} {\bf u}_\theta - r \left(\frac{d\theta}{dt} \right)^2 {\bf u}_r$

$= \displaystyle \left[ \frac{d^2r}{dt^2} - r \left(\frac{d\theta}{dt} \right)^2 \right] {\bf u}_r + \left[ 2\frac{dr}{dt} \frac{d\theta}{dt} + r \frac{d^2\theta}{dt^2} \right] {\bf u}_\theta$ .

This will be needed in the next post, when we use both Newton’s Second Law and Newton’s Law of Gravitation, expressed in polar coordinates.