Lagrange Points and Polynomial Equations: Part 4

From Wikipedia, Lagrange points are points of equilibrium for small-mass objects under the gravitational influence of two massive orbiting bodies. There are five such points in the Sun-Earth system, called L_1, L_2, L_3, L_4, and L_5.

The stable equilibrium points L_4 and L_5 are easiest to explain: they are the corners of equilateral triangles in the plane of Earth’s orbit. The points L_1 and L_2 are also equilibrium points, but they are unstable. Nevertheless, they have practical applications for spaceflight.

As we’ve seen, the positions of L_1 and L_2 can be found by numerically solving the fifth-order polynomial equations

t^5 - (3-\mu) t^4 + (3-2\mu)t^3 - \mu t^2 + 2\mu t - \mu = 0

and

t^5 + (3-\mu) t^4 + (3-2\mu)t^3 - \mu t^2 - 2\mu t - \mu = 0,

respectively. In these equations, \mu = \displaystyle \frac{m_2}{m_1+m_2} where m_1 is the mass of the Sun and m_2 is the mass of Earth. Also, t is the distance from the Earth to L_1 or L_2 measured as a proportion of the distance from the Sun to Earth.

We’ve also seen that, for the Sun and Earth, mu \approx 3.00346 \times 10^{-6}, and numerically solving the above quintics yields t \approx 0.000997 for L_1 and t \approx 0.01004 for L_2. In other words, L_1 and L_2 are approximately the same distance from Earth but in opposite directions.

There’s a good reason why the positive real roots of these two similar quintics are almost equal. We know that t will be a lot closer to 0 than 1 because, for gravity to balance, the Lagrange points have to be a lot closer to Earth than the Sun. For this reason, the terms \mu t^2 and 2\mu t will be a lot smaller than \mu, and so those two terms can be safely ignored in a first-order approximation. Also, the terms t^5 and (3-\mu)t^4 will be a lot smaller than (3-2\mu)t^3, and so those two terms can also be safely ignored in a first-order approximation. Furthermore, since \mu is also close to 0, the coefficient (3-2\mu) can be safely replaced by just 3.

Consequently, the solution of both quintic equations should be close to the solution of the cubic equation

3t^3  - \mu = 0,

which is straightforward to solve:

3t^3 = \mu

t^3 = \displaystyle \frac{\mu}{3}

t = \displaystyle \sqrt[3]{ \frac{\mu}{3} }.

If \mu = 3.00346 \times 10^{-6}, we obtain t \approx 0.010004, which is indeed reasonably close to the actual solutions for L_1 and L_2. Indeed, this may be used as the first approximation in Newton’s method to quickly numerically evaluate the actual solutions of the two quintic polynomials.

Lagrange Points and Polynomial Equations: Part 3

From Wikipedia, Lagrange points are points of equilibrium for small-mass objects under the gravitational influence of two massive orbiting bodies. There are five such points in the Sun-Earth system, called L_1, L_2, L_3, L_4, and L_5.

The stable equilibrium points L_4 and L_5 are easiest to explain: they are the corners of equilateral triangles in the plane of Earth’s orbit. The points L_1 and L_2 are also equilibrium points, but they are unstable. Nevertheless, they have practical applications for spaceflight.

The position of L_2 can be found by numerically solving the fifth-order polynomial equation

(m_1+m_2)t^5+(3m_1+2m_2)t^4+(3m_1+m_2)t^3

-(m_2+3m_3)t^2-(2m_2+3m_3)t-(m_2+m_3)=0.

In this equation, m_1 is the mass of the Sun, m_2 is the mass of Earth, m_3 is the mass of the spacecraft, and t is the distance from the Earth to L_2 measured as a proportion of the distance from the Sun to Earth. In other words, if the distance from the Sun to Earth is 1 unit, then the distance from the Earth to L_2 is t units. The above equation is derived using principles from physics which are not elaborated upon here.

We notice that the coefficients of t^5, t^4, and t^3 are all positive, while the coefficients of t^2, t, and the constant term are all negative. Therefore, since there is only one change in sign, this equation has only one positive real root by Descartes’ Rule of Signs.

Since m_3 is orders of magnitude smaller than both m_1 and m_2, this may safely approximated by

(m_1+m_2)t^5+(3m_1+2m_2)t^4+(3m_1+m_2)t^3 - m_2 t^2- 2m_2 t-m_2=0.

This new equation can be rewritten as

t^5 + \displaystyle \frac{3m_1 + 2m_2}{m_1 + m_2} t^4 + \frac{3m_1+m_2}{m_1+m_2} t^3 - \frac{m_2}{m_1+m_2} t^2 - \frac{2m_2}{m_1+m_2} t- \frac{m_2}{m_1+m_2} = 0.

If we define

\mu = \displaystyle \frac{m_2}{m_1+m_2},

we see that

\displaystyle \frac{3m_1 + 2m_2}{m_1 + m_2} = \frac{3m_1 + 3m_2}{m_1 + m_2} - \frac{m_2}{m_1 + m_2} = 3 - \mu

and

\displaystyle \frac{3m_1 + m_2}{m_1 + m_2} = \frac{3m_1 + 3m_2}{m_1 + m_2} - \frac{2m_2}{m_1 + m_2} = 3 - 2\mu,

so that the equation may be written as

t^5 + (3-\mu) t^4 + (3-2\mu) - \mu t^2 - 2\mu t - \mu = 0,

matching the equation found at Wikipedia.

For the Sun and Earth, m_1 \approx 1.9885 \times 10^{30} ~ \hbox{kg} and m_2 \approx 5.9724 \times 10^{24} ~ \hbox{kg}, so that

mu = \displaystyle \frac{5.9724 \times 10^{24}}{1.9885 \times 10^{30} + 5.9724 \times 10^{24}} \approx 3.00346 \times 10^{-6}.

This yields a quintic equation that is hopeless to solve using standard techniques from Precalculus, but the root can be found graphically by seeing where the function crosses the x-axis (or, in this case, the t-axis):

As it turns out, the root is t \approx 0.01004, so that L_2 is located 1.004\% of the distance from the Earth to the Sun in the direction away from the Sun.

Lagrange Points and Polynomial Equations: Part 2

From Wikipedia, Lagrange points are points of equilibrium for small-mass objects under the gravitational influence of two massive orbiting bodies. There are five such points in the Sun-Earth system, called L_1, L_2, L_3, L_4, and L_5.

The stable equilibrium points L_4 and L_5 are easiest to explain: they are the corners of equilateral triangles in the plane of Earth’s orbit. The points L_1 and L_2 are also equilibrium points, but they are unstable. Nevertheless, they have practical applications for spaceflight.

We begin with L_1, whose position can be found by numerically solving the fifth-order polynomial equation

(m_1+m_3)x^5+(3m_1+2m_3)x^4+(3m_1+m_3)x^3

-(3m_2+m_3)x^2-(3m_2+m_3)x-(m_2+m_3)=0.

In this equation, m_1 is the mass of the Sun, m_2 is the mass of Earth, m_3 is the mass of the spacecraft, and x is the distance from the Earth to L_1 measured as a proportion of the distance from the Sun to L_1. In other words, if the distance from the Sun to L_1 is 1 unit, then the distance from the Earth to L_1 is x units. The above equation is derived using principles from physics which are not elaborated upon here.

We notice that the coefficients of x^5, x^4, and x^3 are all positive, while the coefficients of x^2, x, and the constant term are all negative. Therefore, since there is only one change in sign, this equation has only one positive real root by Descartes’ Rule of Signs.

Since m_3 is orders of magnitude smaller than both m_1 and m_2, this may safely approximated by

m_1 x^5 + 3m_1 x^4 + 3m_1 x^3 - 3m_2 x^2 - 3m_2x - m_2=0.

Unfortunately, the unit x is not as natural for Earth-bound observers as t, the proportion of the distance of L_1 to Earth as a proportion of the distance from the Sun to Earth. Since L_1 is between the Sun and Earth, the distance from the Sun to Earth is x+1 units, so that t = x/(x+1). We then solve for x in terms of t (just like finding an inverse function):

t = \displaystyle \frac{x}{x+1}

t(x+1) = x

tx + t = x

t = x - tx

t= x(1-t)

\displaystyle \frac{t}{1-t} = x.

Substituting into the above equation, we find an equation for t:

\displaystyle \frac{m_1t^5}{(1-t)^5}  + \frac{3m_1t^4}{(1-t)^4} + \frac{3m_1t^3}{(1-t)^3} - \frac{3m_2t^2}{(1-t)^2} -  \frac{3m_2t}{1-t} - m_2=0

m_1t^5  + 3m_1t^4(1-t) + 3m_1t^3(1-t)^2 - 3m_2t^2(1-t)^3 -  3m_2t(1-t)^4 - m_2(1-t)^5=0

Expanding, we find

m_1 t^5 + 3m_1 (t^4 - t^5) + 3m_1 (t^3-2t^4+t^5) - 3m_2 (t^2-3t^3+3t^4-t^5)

-3m_2(t - 4t^2 + 6t^3 - 4t^4 + t^5) - m_2(1 - 5t + 10t^2 - 10 t^3 + 5t^4 + t^5) = 0

Collecting like terms, we find

(m_1 - 3m_1 + 3m_1 + 3m_2 - 3m_2 + m_2)t^5 + (3m_1-6m_1-9m_2+12m_2-5m_2)t^4

+ (3m_1+9m_2-18m_2+10m_2)t^3 + (-3m_2+12m_2-10m_2) t^2

+ (-3m_2+5m_2)t - m_2 = 0,

or

(m_1+m_2) t^5 - (3m_1 +2m_2) t^4 + (3m_1 + m_2) t^3 - m_2 t^2 + 2m_2 t- m_2 = 0.

Again, this equation has only one positive real root since the original quintic in x only had one positive real root. This new equation can be rewritten as

t^5 - \displaystyle \frac{3m_1 + 2m_2}{m_1 + m_2} t^4 + \frac{3m_1+m_2}{m_1+m_2} t^3 - \frac{m_2}{m_1+m_2} t^2 + \frac{2m_2}{m_1+m_2} t- \frac{m_2}{m_1+m_2} = 0.

If we define

\mu = \displaystyle \frac{m_2}{m_1+m_2},

we see that

\displaystyle \frac{3m_1 + 2m_2}{m_1 + m_2} = \frac{3m_1 + 3m_2}{m_1 + m_2} - \frac{m_2}{m_1 + m_2} = 3 - \mu

and

\displaystyle \frac{3m_1 + m_2}{m_1 + m_2} = \frac{3m_1 + 3m_2}{m_1 + m_2} - \frac{2m_2}{m_1 + m_2} = 3 - 2\mu,

so that the equation may be written as

t^5 + (\mu-3) t^4 + (3-2\mu) - \mu t^2 + 2\mu t - \mu = 0,

matching the equation found at Wikipedia.

For the Sun and Earth, m_1 \approx 1.9885 \times 10^{30} ~ \hbox{kg} and m_2 \approx 5.9724 \times 10^{24} ~ \hbox{kg}, so that

\mu = \displaystyle \frac{5.9724 \times 10^{24}}{1.9885 \times 10^{30} + 5.9724 \times 10^{24}} \approx 3.00346 \times 10^{-6}.

This yields a quintic equation that is hopeless to solve using standard techniques from Precalculus, but the root can be found graphically by seeing where the function crosses the x-axis (or, in this case, the t-axis):

As it turns out, the root is t \approx 0.00997, so that L_1 is located 0.997\% of the distance from the Earth to the Sun in the direction of the Sun.

Lagrange Points and Polynomial Equations: Part 1

I recently read a terrific article in the American Mathematical Monthly about Lagrange points, which are (from Wikipedia) “points of equilibrium for small-mass objects under the gravitational influence of two massive orbiting bodies.” There are five such points in the Sun-Earth system, called L_1, L_2, L_3, L_4, and L_5.

To describe these Lagrange points, I can do no better than the estimable Isaac Asimov. I quote from his essay “Colonizing the Heavens” from his book The Beginning and the End, which was published in 1977. I read the book over and over again as a boy in the mid-1980s. (Asimov’s essay originally concerned the Earth-Moon system; in the quote below, I changed the words to apply to the Sun-Earth system.)

Imagine the Sun at zenith, exactly overhead. Trace a line due eastward from the Sun down to the horizon. Two-thirds of the way along that line, one-third of the way up from the horizon, is one of those places. Trace another line westward away from the Sun down to the horizon. Two-thirds of the way along that line, one-third of the way up from the horizon, is another of those places.

Put an object in either place and it will form an equilateral triangle with the Sun and Earth…

What is so special about those places? Back in 1772, the astronomer Joseph Louis Lagrange showed that in those places any object remained stationary with respect to the Sun. As the Earth moved about the Sun, any object in either of those places would also move about the Sun in such a way as to keep perfect step with the Earth. The competing gravities of the Sun and Earth would keep it where it was. If anything happened to push it out of place it would promptly move back, wobbling back and forth a bit (“librating”) as it did so. The two places are called “Lagrangian points” or “libration points.”

Lagrange discovered five such places altogether, but three of them are of no importance since they don’t represent stable conditions. An object in those three places, once pushed out of place, would continue to drift outward and would never return.

The last paragraph of the above quote represents a rare failure of imagination by Asimov, who wrote prolifically about the future of spaceflight. Points L_4 and L_5 are indeed stable equilibria, and untold science fiction stories have placed spacecraft or colonies at these locations. (The rest of Asimov’s essay speculates about using these points in the Earth-Moon system for space colonization.) However, while the points L_1 and L_2 are unstable equilibria, they do have practical applications for spacecraft that can perform minor course corrections to stay in position. (The point L_3 is especially unstable to outside gravitational influences and thus seems unsuitable for spacecraft.) Again from Wikipedia,

Sun–Earth L1 is suited for making observations of the Sun–Earth system. Objects here are never shadowed by Earth or the Moon and, if observing Earth, always view the sunlit hemisphere… Solar and heliospheric missions currently located around L1 include the Solar and Heliospheric Observatory, Wind, Aditya-L1 Mission and the Advanced Composition Explorer. Planned missions include the Interstellar Mapping and Acceleration Probe(IMAP) and the NEO Surveyor.

Sun–Earth L2 is a good spot for space-based observatories. Because an object around L2 will maintain the same relative position with respect to the Sun and Earth, shielding and calibration are much simpler… The James Webb Space Telescope was positioned in a halo orbit about L2 on January 24, 2022.

Earth–Moon L1 allows comparatively easy access to Lunar and Earth orbits with minimal change in velocity and this has as an advantage to position a habitable space station intended to help transport cargo and personnel to the Moon and back. The SMART-1 Mission passed through the L1 Lagrangian Point on 11 November 2004 and passed into the area dominated by the Moon’s gravitational influence.

Earth–Moon L2 has been used for a communications satellite covering the Moon’s far side, for example, Queqiao, launched in 2018, and would be “an ideal location” for a propellant depot as part of the proposed depot-based space transportation architecture.

While the locations L_4 and L_5 are easy to describe, the precise locations of L_1 and L_2 are found by numerically solving a fifth-order polynomial equation. This was news to me when I read that article from the American Mathematical Monthly. While I had read years ago that finding the positions of the other three Lagrange points wasn’t simple, I did not realize that it was no more complicated that numerically finding the roots of a polynomial.

The above article from the American Mathematical Monthly concludes…

[t]he mathematical tools that Lagrange uses to arrive at a solution to this three-body problem lie entirely within the scope of modern courses in algebra, trigonometry, and first-semester calculus. But surely no ordinary person could have pursued the many extraordinarily complicated threads in his work to their ends, let alone woven them together into a magnificent solution to the problem as he has done. Lagrange noted in the introduction to his essay, “This research is really no more than for pure curiosity …” If only he could have watched on Christmas Day as the James Webb Space Telescope began its journey to the Lagrange point L_2!

In this short series, we discuss the polynomial equations for finding L_1 and L_2.

Confirming Einstein’s Theory of General Relativity With Calculus, Part 6i: Rationale for Method of Undetermined Coefficients VI

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.

In the last few posts, I’ve used a standard technique from differential equations: to solve the nth order homogeneous differential equation with constant coefficients

a_n y^{(n)} + \dots + a_3 y''' + a_2 y'' + a_1 y' + a_0 y = 0,

we first solve the characteristic equation

a_n r^n + \dots + a_3 r^3 + a_2 r^2 + a_1 r + a_0 = 0

using techniques from Precalculus. The form of the roots r determines the solutions of the differential equation.

While this is a standard technique from differential equations, the perspective I’m taking in this series is scaffolding the techniques used to predict the precession in a planet’s orbit using only techniques from Calculus and Precalculus. So let me discuss why the above technique works, assuming that the characteristic equation does not have repeated roots. (The repeated roots case is a little more complicated but is not needed for the present series of posts.)

We begin by guessing that the above differential equation has a solution of the form y = e^{rt}. Differentiating, we find y' = re^{rt}, y'' = r^2 e^{rt}, etc. Therefore, the differential equation becomes

a_n r^n e^{rt} + \dots + a_3 r^3 e^{rt} + a_2 r^2 e^{rt} + a_1 r e^{rt} + a_0 e^{rt} = 0

e^{rt} \left(a_n r^n  + \dots + a_3 r^3 + a_2 r^2 + a_1 r  + a_0 \right) = 0

a_n r^n  + \dots + a_3 r^3 + a_2 r^2 + a_1 r  + a_0 = 0

The last step does not “lose” any possible solutions for r since e^{rt} can never be equal to 0. Therefore, solving the differential equation reduces to finding the roots of this polynomial, which can be done using standard techniques from Precalculus.

For example, one of the differential equations that we’ve encountered is y''+y=0. The characteristic equation is r^2+1=0, which has roots r=\pm i. Therefore, two solutions to the differential equation are e^{it} and e^{-it}, so that the general solution is

y = c_1 e^{it} + c_2 e^{-it}.

To write this in a more conventional way, we use Euler’s formula e^{ix} = \cos x + i \sin x, so that

y = c_1 (\cos t + i \sin t) + c_2 (\cos (-t) + i \sin (-t))

= c_1 \cos t + i c_1 \sin t + c_2 \cos t - i c_2 \sin t

= (c_1 + c_2) \cos t + (ic_1 - ic_2) \sin t

= C_1 \cos t + C_2 \sin t.

Likewise, in the previous post, we encountered the fourth-order differential equation y^{(4)}+5y''+4y = 0. To find the roots of the characteristic equation, we factor:

r^4 + 5r^2 + 4r = 0

(r^2+1)(r^2+4) = 0

r^2 +1 = 0 \qquad \hbox{or} \qquad \hbox{or} r^2 + 4 = 0

r = \pm i \qquad \hbox{or} \qquad r = \pm 2i.

Therefore, four solutions of this differential equation are e^{it}, e^{-it}, e^{2it}, and e^{-2it}, so that the general solution is

y = c_1 e^{it} + c_2 e^{-it} + c_3 e^{2it} + c_4 e{-2it}.

Using Euler’s formula as before, this can be rewritten as

y = C_1 \cos t + C_2 \sin t + C_3 \cos 2t + C_4 \sin 2t.

Confirming Einstein’s Theory of General Relativity With Calculus, Part 6g: Rationale for Method of Undetermined Coefficients IV

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.

We have shown that the motion of a planet around the Sun, expressed in polar coordinates (r,\theta) with the Sun at the origin, under general relativity follows the initial-value problem

u''(\theta) + u(\theta) = \displaystyle \frac{1}{\alpha} + \frac{\delta}{\alpha^2} + \frac{\delta \epsilon^2}{2\alpha^2} + \frac{2\delta \epsilon \cos \theta}{\alpha^2} + \frac{\delta \epsilon^2 \cos 2\theta}{2\alpha^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).

In this post, we will use the guesses

u(\theta) = f(\theta) \cos \theta \qquad \hbox{or} u(\theta) = f(\theta) \sin \theta

that arose from the technique/trick of reduction of order, where f(\theta) is some unknown function, to find the general solution of the differential equation

u^{(4)} + 2u'' + u = 0.

To do this, we will need to use the Product Rule for higher-order derivatives that was derived in the previous post:

(fg)'' = f'' g + 2 f' g' + f g''

and

(fg)^{(4)} = f^{(4)} g + 4 f''' g' + 6 f'' g'' + 4f' g''' + f g^{(4)}.

In these formulas, Pascal’s triangle makes a somewhat surprising appearance; indeed, this pattern can be proven with mathematical induction.

We begin with u(\theta) = f(\theta) \cos \theta. If g(\theta) = \cos \theta, then

g'(\theta) = - \sin \theta,

g''(\theta) = -\cos \theta,

g'''(\theta) = \sin \theta,

g^{(4)}(\theta) = \cos \theta.

Substituting into the fourth-order differential equation, we find the differential equation becomes

(f \cos \theta)^{(4)} + 2 (f \cos \theta)'' + f \cos \theta = 0

f^{(4)} \cos \theta - 4 f''' \sin \theta - 6 f'' \cos \theta + 4 f' \sin \theta + f \cos \theta + 2 f'' \cos \theta - 4 f' \sin \theta - 2 f \cos \theta + f \cos \theta = 0

f^{(4)} \cos \theta - 4 f''' \sin \theta - 6 f'' \cos \theta  + 2 f'' \cos \theta  = 0

f^{(4)} \cos \theta - 4 f''' \sin \theta - 4 f'' \cos \theta = 0

The important observation is that the terms containing f and f' cancelled each other. This new differential equation doesn’t look like much of an improvement over the original fourth-order differential equation, but we can make a key observation: if f'' = 0, then differentiating twice more trivially yields f''' = 0 and f^{(4)} = 0. Said another way: if f'' = 0, then u(\theta) = f(\theta) \cos \theta will be a solution of the original differential equation.

Integrating twice, we can find f:

f''(\theta) = 0

f'(\theta) = c_1

f(\theta) = c_1 \theta + c_2.

Therefore, a solution of the original differential equation will be

u(\theta) = c_1 \theta \cos \theta + c_2 \cos \theta.

We now repeat the logic for u(\theta) = f(\theta) \sin \theta:

(f \sin \theta)^{(4)} + 2 (f \sin \theta)'' + f \sin \theta = 0

f^{(4)} \sin \theta + 4 f''' \cos \theta - 6 f'' \sin \theta - 4 f' \cos\theta + f \sin \theta + 2 f'' \sin \theta + 4 f' \cos \theta - 2 f \sin \theta + f \sin \theta = 0

f^{(4)} \sin\theta + 4 f''' \cos \theta - 6 f'' \sin \theta + 2 f'' \sin \theta = 0

f^{(4)} \sin\theta - 4 f''' \cos\theta - 4 f'' \sin\theta = 0.

Once again, a solution of this new differential equation will be f(\theta) = c_3 \theta + c_4, so that f'' = f''' = f^{(4)} = 0. Therefore, another solution of the original differential equation will be

u(\theta) = c_3 \theta \sin \theta + c_4 \sin \theta.

Adding these provides the general solution of the differential equation:

u(\theta) = c_1 \theta \cos \theta + c_2 \cos \theta + c_3 \theta \sin \theta + c_4 \sin \theta.

Except for the order of the constants, this matches the solution that was presented earlier by using techniques taught in a proper course in differential equations.

Confirming Einstein’s Theory of General Relativity With Calculus, Part 6f: Rationale for Method of Undetermined Coefficients III

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.

We have shown that the motion of a planet around the Sun, expressed in polar coordinates (r,\theta) with the Sun at the origin, under general relativity follows the initial-value problem

u''(\theta) + u(\theta) = \displaystyle \frac{1}{\alpha} + \frac{\delta}{\alpha^2} + \frac{\delta \epsilon^2}{2\alpha^2} + \frac{2\delta \epsilon \cos \theta}{\alpha^2} + \frac{\delta \epsilon^2 \cos 2\theta}{2\alpha^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).

In the previous post, I used a standard technique from differential equations to find the general solution of

u^{(4)} + 2u'' + u = 0.

to be

u(theta) = c_1 \cos \theta + c_2 \sin \theta + c_3 \theta \cos \theta + c_4 \theta \sin \theta.

However, as much as possible in this series, I want to take the perspective of a talented calculus student who has not yet taken differential equations — so that the conclusion above is far from obvious. How could this be reasonable coaxed out of such a student?

To begin, we observe that the characteristic equation is

r^4 + 2r^2 + 1 = 0,

or

(r^2 + 1)^2 = 0.

Clearly this has the same roots as the simpler equation r^2 + 1 = 0, which corresponds to the second-order differential equation u'' + u = 0. We’ve already seen that u_1(\theta) = \cos \theta and u_2(\theta) = \sin \theta are solutions of this differential equation; perhaps they might also be solutions of the more complicated differential equation also? The answer, of course, is yes:

u_1^{(4)} + 2 u_1'' + u_1 = \cos \theta - 2 \cos \theta + \cos \theta = 0

and

u_2^{(4)} + 2u_2'' + u_2 = \sin \theta - 2 \sin \theta + \sin \theta = 0.

The far trickier part is finding the two additional solutions. To find these, we use a standard trick/technique called reduction of order. In this technique, we guess that any additional solutions much have the form of either

u(\theta) = f(\theta) \cos \theta \qquad \hbox{or} \qquad  u(\theta) = f(\theta) \sin \theta,

where f(\theta) is some unknown function that we’re multiplying by the solutions we already have. We then substitute this into the differential equation u^{(4)} + 2u'' + u = 0 to form a new differential equation for the unknown f, which we can (hopefully) solve.

Doing this will require multiple applications of the Product Rule for differentiation. We already know that

(fg)' = f' g + f g'.

We now differentiate again, using the Product Rule, to find (fg)'':

(fg)'' = ( [fg]')' = (f'g)' + (fg')'

= f''g + f' g' + f' g' + f g''

= f'' g + 2 f' g' + f g''.

We now differential twice more to find (fg)^{(4)}:

(fg)''' = ( [fg]'')' = (f''g)' + 2(f'g')' +  (fg'')'

= f'''g + f'' g' + 2f'' g' + 2f' g'' + f' g'' + f g'''

= f''' g + 3 f'' g' + 3 f' g'' + f g'''.

A good student may be able to guess the pattern for the next derivative:

(fg)^{(4)} = ( [fg]''')' = (f'''g)' + 3(f''g')' +3(f'g'')' + (fg''')'

= f^{(4)}g + f''' g' + 3f''' g' + 3f'' g'' + 3f'' g'' + 3f'g''' + f' g''' + f g^{(4)}

= f^{(4)} g + 4 f''' g' + 6 f'' g'' + 4f' g''' + f g^{(4)}.

In this way, Pascal’s triangle makes a somewhat surprising appearance; indeed, this pattern can be proven with mathematical induction.

In the next post, we’ll apply this to the solution of the fourth-order differential equation.

Confirming Einstein’s Theory of General Relativity With Calculus, Part 6e: Rationale for Method of Undetermined Coefficients II

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.

We have shown that the motion of a planet around the Sun, expressed in polar coordinates (r,\theta) with the Sun at the origin, under general relativity follows the initial-value problem

u''(\theta) + u(\theta) = \displaystyle \frac{1}{\alpha} + \frac{\delta}{\alpha^2} + \frac{\delta \epsilon^2}{2\alpha^2} + \frac{2\delta \epsilon \cos \theta}{\alpha^2} + \frac{\delta \epsilon^2 \cos 2\theta}{2\alpha^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).

In the previous post, we derived the method of undetermined coefficients for the simplified differential equation

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

In this post, we consider the simplified differential equation if the right-hand side has only the fourth term,

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

Let v(\theta) =  \displaystyle \frac{2\delta \epsilon }{\alpha^2}\cos \theta. Then v satisfies the new differential equation v'' + v = 0. Since u'' + u = v, we may substitute:

(u''+u)'' + (u'' + u) = 0

u^{(4)} + u'' + u'' + u = 0

u^{(4)} + 2u'' + u = 0.

The characteristic equation of this homogeneous differential equation is r^4 + 2r^2 + 1 = 0, or (r^2+1)^2 = 0. Therefore, r = i and r = -i are both double roots of this quartic equation. Therefore, the general solution for u is

u(\theta) = c_1 \cos \theta + c_2 \sin \theta + c_3 \theta \cos \theta + c_4 \theta \sin \theta.

Substituting into the original differential equation will allow for the computation of c_3 and c_4:

u''(\theta) + u(\theta) = -c_1 \cos \theta - c_2 \sin \theta - 2c_3 \sin \theta - c_3 \theta \cos \theta + 2c_4 \cos \theta - c_4 \theta \sin \theta

+   c_1 \cos \theta + c_2 \sin \theta + c_3 \theta \cos \theta + c_4 \theta \sin \theta

\displaystyle \frac{2\delta \epsilon }{\alpha^2}\cos \theta = - 2c_3 \sin \theta+ 2c_4 \cos \theta

Matching coefficients, we see that c_3 = 0 and c_4 = \displaystyle \frac{\delta \epsilon }{\alpha^2}. Therefore,

u(\theta) = c_1 \cos \theta + c_2 \sin \theta + \displaystyle \frac{\delta \epsilon }{\alpha^2} \theta \sin \theta

is the general solution of the simplified differential equation. Setting c_1 = c_2 = 0, we find that

u(\theta) =  \displaystyle \frac{\delta \epsilon }{\alpha^2} \theta \sin \theta

is one particular solution of this simplified differential equation. Not surprisingly, this matches the result is the method of undetermined coefficients had been blindly followed.

As we’ll see in a future post, the presence of this \theta \sin \theta term is what predicts the precession of a planet’s orbit under general relativity.

Confirming Einstein’s Theory of General Relativity With Calculus, Part 6d: Rationale for Method of Undetermined Coefficeints I

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.

We have shown that the motion of a planet around the Sun, expressed in polar coordinates (r,\theta) with the Sun at the origin, under general relativity 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,

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 constant coefficients. To use this technique, we first expand the right-hand side of the differential equation and then apply a power-reduction trigonometric identity:

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

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

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

This is now in the form for using the method of undetermined coefficients. However, in this series, I’d like to take some time to explain why this technique actually works. To begin, we look at a simplified differential equation using only the first three terms on the right-hand side:

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

Let v(\theta) =\displaystyle \frac{1}{\alpha} + \frac{\delta}{\alpha^2} + \frac{\delta \epsilon^2}{2\alpha^2}  . Since v is a constant, this function satisfies the simple differential equation v' = 0. Since u''+u=v, we can substitute:

(u'' + u)' = 0

u''' + u' = 0

(We could have more easily said, “Take the derivative of both sides,” but we’ll be using a more complicated form of this technique in future posts.) The characteristic equation of this differential equation is r^3 + r = 0. Factoring, we obtain r(r^2 + 1) = 0, so that the three roots are r = 0 and r = \pm i. Therefore, the general solution of this differential equation is

u(\theta) = c_1 \cos \theta + c_2 \sin \theta + c_3.

Notice that this matches the outcome of blindly using the method of undetermined coefficients without conceptually understanding why this technique works.

The constants c_1 and c_2 are determined by the initial conditions. To find c_3, we observe

u''(\theta) +u(\theta) =  -c_1 \cos \theta - c_2 \sin \theta +c_1 \cos \theta + c_2 \sin \theta + c_3

\displaystyle \frac{1}{\alpha} + \frac{\delta}{\alpha^2} + \frac{\delta \epsilon^2}{2\alpha^2}  = c_3.

Therefore, the general solution of this simplified differential equation is

u(\theta) = c_1 \cos \theta + c_2 \sin \theta + \displaystyle \frac{1}{\alpha} + \frac{\delta}{\alpha^2} + \frac{\delta \epsilon^2}{2\alpha^2}.

Furthermore, setting c_1 = c_2 = 0, we see that

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

is a particular solution to the differential equation

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

In the next couple of posts, we find the particular solutions associated with the other terms on the right-hand side.