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)

is

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

= 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^2 \theta\cos 2 \theta}{\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.

Leave a comment

This site uses Akismet to reduce spam. Learn how your comment data is processed.