Topics

- Historical modelling on the estimation of the figure of the earth
- Use of calculus and series to derive formulas

Further reading

- The quest for the true figure of the Earth : ideas and expeditions in four centuries of geodesy by Michael Rand Hoare, 2004. ISBN:0754650200
- The problem of the earth's shape from Newton to Clairaut by John Leonard Greenberg, 1995. ISBN:0521385415
- Stanford Lecture: Donald Knuth - "The Joy of Asymptotics"

One of the great tools of applied mathematics has been the *art* of approximation. As you may have noticed in the past, complicated equations often have special cases that are easier to solve. A classical example would be a hyperbola like \[x^2 + 0.01 x -1 = 0\] can be solved by quadratic formula or completing the square to show \(x \approx -1.00501\) or \(x \approx 0.99501\). But this equation is also "close" to the equation \(x^2 - 1 = 0\), which we can solve immediately by factorization to get \(x \in \{-1,1\}\), values that differ from the exact solutions by fractions of a percent. This is a toy example -- we don't really need the simplify anything to find the solutions. But wouldn't it be nice if we had a systematic way to derive approximate equations whose solutions were stayed near to the solutions of the original equation? That's approximation theory. Here, we'll see some applications of the simpliest form of analytic approximation -- Taylor series -- to geometric models for Watt's curve and inference of the earth's figure.

The theory of approximation has roots in the ancient study of hyperbola's by Apollonius of Perga. Here's the core idea. If you want to calculate \(y\) but it's too hard, find something nearby like \(x\) that is easier to calculate and then use a method to estimate the value of \(y\) based on \(x\). The name comes from the lines used to bound the shape of a hyperbola. The most basic form of analytic approximation is a Taylor series.

Let's consider the quadratic equation given above. We can introduce a small parameter into our original equation by letting \(\epsilon = 0.01\) so that \[x^2 + \epsilon x - 1 = 0.\] Now, let's treat our unknown as a function the small parameter, \(x(\epsilon)\), and imagine that this function has a McLaurin series of the form \[x(\epsilon) = x_0 + \epsilon x_1 + \epsilon^2 x_2 \ldots.\] If we substitute this series into our quadratic equation, we will get a polynomial in \(\epsilon\). Since this polynomial should equal \(0\) for all small values of \(\epsilon\), then the coefficient of each \(\epsilon\)-monomial should also vanish. \[\begin{gather} 0 = (x_0^2 - 1) + \epsilon x_0 (2 x_1 + 1) + \epsilon^2 ( 2 x_0 x_2 + x_1^2 + x_1 ) + \ldots \end{gather}\]Since \(\epsilon\) is small, the most important term is \(x_0^2 - 1\), which vanishes for \(x_0 = 1\) or \(x_0 = -1\), the two values we initially guessed. The next most important term is \(x_0 (2 x_1 + 1)\), which will vanish for \(x_1 = -1/2\). Thus, this approximation method suggests \[x(\epsilon) \approx \pm 1 - \frac{1}{2} \epsilon.\] This very simple approximation gives \(x\approx -1.00500\) and \(x\approx 0.99500\), 4 decimal places of accuracy when \(\epsilon = 1/100\).

To determine the parameter conditions when Watt's linkage mobility is approximately straight around its midpoint, we can apply a series approximation. Recall that the implicit formula for Watt's curve is \[0 = 4 y^{2} \left(x^{2} + y^{2} - r^{2}\right) + \left(x^{2} + y^{2}\right) \left(\frac{\ell^{2}}{4} - r^{2} - 1 + x^{2} + y^{2} \right)^{2}.\] We are interested in Watt's special case where \(\ell = 2 \sqrt{1 - r^2}\), when the piston motion is said to be approximately straight, and the implicit configuration curve simplifies to \[0 = 4 y^{2} \left(x^{2} + y^{2} - r^{2}\right) + \left(x^{2} + y^{2}\right) \left(x^{2} + y^{2} - 2 r^{2}\right)^{2}.\] We can check by substitution and inspection that \((x,y) = (0,0)\) lies on the curve. Now, using a substituting \[y \approx a_1 x + a_2 x^2 + a_3 x^3 + a_4 x^4 + a_5 x^5,\] and picking the \(a_i\)'s so that coefficients vanish one by one, we can show that link is approximately give by \[y(x) = a_1 x + a_5 x^5 + O(x^6), \quad a_1 = \pm \frac{r}{\sqrt{1-r^2}}, \quad a_5 = \frac{-(1 + a_1^2)^3}{8 a_1 r^2 (1 - r^2)}.\] This shows that Watt's linkage motion is approximately a straight-line motion to fifth-order when the linking bar length \(\ell\) is picked right.

The hypothesis that earth is a ball with no edges traces back at least as far as the pre-Socratic Greek philosophers(c. 500 BC) but was first confirmed by the Magellan-Elcano circumnavigation completed in 1522. But with the increased scientific sophistication of the age of enlightenment, it became a matter of great debate if the earth's figure was actually spherical, and if not, was it more like a lemon sinched around the equator or an onion squached at the poles? How would you go about determining the figure of the earth? If the earth's cross-section is not a circle, how can we tell? It was easy to answer this question for the planet Jupiter -- astronomers can directly observe the figure of Jupiter through a telescope. Giovanni Cassini did this and determined that Jupiter was wider than tall. But we can not do the same for the planet we live on. This puzzle forces us to take a closer examination of what we mean by "figure" and "shape".

Our descriptive language of the differences between the shapes of a lemon and onion suggests a way forward. An important local characteristic of shape is curvature -- how unlike a straight line a curve or surface. A circle has the same curvature all the way around as you go, but other shapes like a lemon become more or less curved and for some shapes even reverse the direction of curvature. If we can measure curvature at different locations, then we might be able to stitch togethor those local curvatures determine the overall shape.

Curvature can be thought of as the rate of direction change per unit of distance travelled. Suppose we start at some point \(p\) on a planar curve. Draw a tangent vector to the curve at the point we are starting. Now, we walk along the curve, keeping track of the distance \(s\) we have travelled. As long as the curve we are walking along is smooth enough that we can measure changes in the tangent vector to the curve, and we have a universal reference frame like the night sky, we can also recording the angle \(\lambda(s)\) between the tangent vector to the curve at our current location and the tangent vector at our start location. The local curvature, then, can be thought of as \(d\lambda/ds\). On a circle, this implies that the curvature is constant everywhere and proportional the reciprocal of the radius. As we might expect, the bigger the circles get, the smaller their curvatures. On a surface, we can get the same results by tracking changes in the direction of an oriented unit normal vector.

Now, let's get down to specifics. In order to apply these curvature ideas, particularly back in the 1600's when this was first done, we have to have a relatively simple parameterizable model that can accomodate both lemon and onion shapes. The only model widely known to mathematicians of the time was that of ellipsoids of revolution -- bodies defined by rotating an ellipse around one of it's axes. All of spheres, lemon, and onion shapes could be constructed by rotating various ellipses. Of course, the earth is not a perfect ellipsoid of revolution -- it has mountains and plateaus and valleys and icecaps -- but that's a natural start.

We can represent an ellipsoid of revolution in terms of an ellipse \[x^2 + (1+\epsilon)y^2 = r^2\] where \(r\) is the radius at the equator and \(\epsilon\) represents the distortion of the ellipse from a perfect circle. When the ellipse is rotated around the vertical axis, it creates a lemon (prolate ellipsoid) when \(\epsilon < 0\) and an onion (oblate ellipsoid) when \(\epsilon > 0\).

The normal vector perpendicular to a point on this ellipse can be found by treating this ellipse as a single contour line of the function \(F(x,y) = x^2 + (1+\epsilon) y^2\). The gradient can be determined by differentiation, \(\nabla F = [ 2 x, 2 (1+ \epsilon ) y ]\), and at each point, the gradient vector will be perpendicular to the ellipse curve forming a contour of constant value. Let \(\lambda\) be the angle between the normal vector and the equatorial plane, *e.g.* \(\lambda\) is the celestial latitude of a position \((x,y)\). Then from elementary trigonometry, \[\tan \lambda = \frac{(1 + \epsilon)y}{x}.\] We can use this equation to represent both \(x\) and \(y\) parametrically as functions of \(\lambda\). Solving for \(y\) and substituting into our equation for the ellipse, \[\begin{gather*}
x^2 + \frac{1}{1 + \epsilon} x^2 \tan^2 \lambda = r^2,
\\
x(\lambda) = r \cos(\lambda) \sqrt{\frac{1+\epsilon}{1+\epsilon \cos^2 \lambda}}.
\end{gather*}\] Similarly, solving our tangent equation for \(x\) and substituting leads to \[y(\lambda) = r \sin(\lambda) \sqrt{\frac{1 }{(1+\epsilon)(1 + \epsilon \cos^{2}(\lambda ))}}.\]

Now, having position and latitude, we would like to relate latitude as a function of arclength to get curvature. However, we can already see it is going to be very hard to solve for latitude in terms of arclength. Instead, let's flip things around -- reciprocal curvature describes the change in arc length as a function of change in latitude. If we parameterize the arclength \(s\) as a function of the latitude \(\lambda\), standard calculus references tell us \[\frac{ds}{d\lambda} = \sqrt{ \left( \frac{dx}{d\lambda} \right)^2 + \left( \frac{dy}{d\lambda} \right)^2 }.\] If we calculate the derivatives, substitute, and simplify, we find the reciprocal curvature \[\frac{ds}{d\lambda} = \sqrt{\frac{(1 + \epsilon) r^2}{\left(1 + \epsilon \cos^2(\lambda)\right)^3}}.\]

This formula is closely related to one obtained by Maupertius in 1736 for variation in the earth's radius as a function of the latitude and elliptic eccentricity.

Maupertius's equation for the reciprocal curvature is a laborious nonlinear equation, so working with it directly at the time was not feasible. But since the earth is approximately spherical, \(\epsilon \approx 0\). The first two terms in the McLaurin series expansion in \(\epsilon\) give \[\frac{ds}{d\lambda} \approx r + \frac{1}{2} \epsilon r - \frac{3}{2} \epsilon r \cos^{2}{\left (\lambda \right )},\] or if we transform to more convenient and conventional form, \[m \sin^2(\lambda) + c \approx \frac{ds}{d\lambda}\] where \(c = r (1 - \epsilon)\) and \(m = 3/2 \epsilon r\). This remarkably little equation relating the latitude and distance on an ellipsoidal planet was applied by Maupertuis, Boscovich, Laplace, and Legendre (with minor variation), when estimating the figure of the earth. The left-hand side has units of distance per angle, the reciprocal of curvature. The ellipse parameters \(r\) and \(\epsilon\) on the right-hand side are the what we seek to estimate, while the free variables \(\lambda\) and \(ds/d\lambda\) can be measured directly by observation.

As we'll see when we start using linear least squares, things aren't quite so simple, but we'll make them work.

The three examples we've given above are solved by what's called "regular approximations" -- we applied some familiar ideas from calculus in a common-sense way, and got out useful approximations without anything weird happening. However, don't get over confident fishing for solutions this way -- the water gets deep quickly and there are monster's leaking below the surface. Consider the deceptive problem of approximating a positive solution to the cubic equation \[0.01 x^3 - x = 1.\] Since this is cubic, we expect three solutions, one of which should be positive. But the obvious approximation leads to \(x = -1\), which has only one solution instead of 3 and *is not even positive*! Sometimes small things are less important, but sometimes small things have big important consequences. You cann't always tell ahead of time which is which. We'll see more examples as we go forward.

The story of the debate over the figure of the earth begins with one of the most practical applied-mathematics problems of the day -- ship navigation, and the longitude problem in-particular. As nations grew and trade expanded during the Renaissance, shipwrecks became an important source of financial and personal risk. The risks became larger after 1492, as explorers discovered new trade routes consisting of long, expensive, and uncertain ocean crossings. Accurate navigational charts were being drawn with the help of trigonometry, but could not help ships unable to determine their own positions in open ocean. Latitude could be approximated based on observations of the stars, but longitude could not because the earth's daily rotation blurred the most accessible celestial references.

In 1656, Christiaan Huygens, a remarkable Dutch natural philosopher whom we will encounter again later, patented the first pendulum clock. One of the uses Huygens proposed for his clock was the determination of longitude. If you know what time it was in New York when it was noon in Paris, you would know the time-difference between the two, and by multiplication, the longitude of New York relative to Paris. Huygen's own clocks were not precise enough to practically solve the longitude problem -- John Harrison's chronometer designs would eventually do the job a century later. But the idea already intrigued many, and in 1670, when the recently created AcadÃ©mie des Sciences in Paris decided to send it's second scientific expedition to Cayenne in French Guyana, clock-testing was part of the plan. Two years later, Jean Richer reported back that the length of a seconds-pendulum was not the same in Cayenne and Paris! In Cayenne, the seconds-pendulum was slightly shorter than its 24.84 centimeters in Paris. Parisian clocks would lose about 2 minutes and 28 seconds each day when used in Cayenne.

This one piece of modest data was the pebble that started an avalanche of debate. Giovanni Cassini, the new influential Director of the Paris observatory, interpreted Richer's observation in the context of Huygen's pendulum theory (see later). There, the slower period implied that gravity was weaker in Cayenne, probably because earth was smaller around the equator than it was around the poles. Alternatively, in 1683, Robert Hooke had hypothesized that the earth was shaped like an onion rather than a sphere. Isaac Newton ran with this idea in the first edition of *Principia* (1687). He argued that the slowing of the pendulum was due to the spin of the earth counteracting gravitational forces, and that the slowing would coincide with a bulging of the earth around the equator. Cassini's and Newton's contradictory positions, proposed by influential scientists on opposite sides of the English channel, made the figure of the earth a hot topic.

In one of his famous annecdotes, Richard Feynman and an abacus salesman were challenged to compute something like \(x = \sqrt[3]{1729.03}\). Feynman observed that this was the same as solving \(x^3 - 12^3 - \epsilon = 0\) when \(\epsilon = 1.03\), and derived a linear approximation to \(x(\epsilon)\approx x_0 + \epsilon x_1\) quickly in his head.

Using our methods above, find \(x_0\) and \(x_1\), estimate \(x\), and determine how accurate Feynman's estimated answer was.Asymptotic approximations can be used to show spherical trigonometry agrees with planar trigonometry when the triangles are much smaller than the radius of the sphere they lie on. Consider a spherical triangle with angles \(A\),\(B\),\(\Gamma\) opposite sides (measured in radians) \(a\),\(b\), and \(c\). Let \(\alpha\), \(\beta\), and \(\gamma\) be the lengths of the sides of a triangle on a sphere of radius \(r\), such that the corresponding angle measures of those sides are \(a = \alpha/r\), \(b = \beta/r\), and \(c=\gamma/r\). Derive approximations for each of the following in the limit of large a sphere (\(r \rightarrow \infty\), or \(\epsilon := 1/r \rightarrow 0\)), and name the associated identity from planar trigonometry.

- \(\sin a \sin B = \sin b \sin A\)
- \(\cos \Gamma = \sin A \sin B \cos c - \cos A \cos B.\)
- \(\cos c = \cos a \cos b\) when the angle \(\Gamma\) opposite \(c\) is \(90^{\circ}\).

- Previously, we found the following implicit solution for the path traced by the handle of Archimedes's trammel was \[r^2 = x^2 + \left(\frac{1}{\frac{s}{r}+1}\right)^2 y^2.\]
- Describe the asymptotic nature of the solution when \(s=1\) and \(r \rightarrow \infty\).
- Describe the asymptotic nature of the solution when \(r=1\) and \(s \rightarrow 0\).
- Describe the asymptotic nature of the solution when \(r=0\).
- Describe the asymptotic nature of the solution when \(r=-s\).
- Describe the asymptotic nature of the solution when \(r=-s/2\).

- The theory of asymptotic approximation has roots in the ancient study of hyperbola's by Apollonius of Perga. Consider the hyperbolas of the form \(x^2 - m^2 y^2 + ax + by + c = 0\), which we would like to approximate when \(x\) and \(y\) are large.
- Using the change-of-coordinates \(u=1/x\) and \(v=1/y\), transform this hyperbola to a quartic polynomial.
- Assume \(u\) and \(v\) are small, discard the smallest terms and find all tangent lines that approximate part of the solution and \((u,v) =(0,0)\).
- Which parameters do the asymptotic approximations depend on, and which do not matter at all? (This is how asymptotic approximations can make problems simplier!)
- Transform your tangent-line approximations back to \(x\) and \(y\) coordinates.
- Use the same style of analysis to obtain asymptotic approximations to \(x^3 - 64 y^2 + 4 = 0\) when \(x\) and \(y\) are large.

- Our introductory example of asymptotic approximation can be arrived at by other means.
- Solve \(x^2 + \epsilon x - 1=0\) for \(\epsilon(x)\).
- Recalling the calculus rules that, \[\frac{dx}{d\epsilon} = \left(\frac{d\epsilon}{dx}\right)^{-1}, \quad \frac{d^2x}{d\epsilon^2} = - \left(\frac{dx}{d\epsilon}\right)^3 \frac{d^2\epsilon}{dx^2},\] use Taylor series to find the first 3 terms of the approximation of \(x(\epsilon)\) when \(x(0) = 1\).

Assuming latitude, measured in degrees, falls between \(0\) and \(90\), plot the maximum error between our exact formula for reciprocal curvature and it's first-order McLaurin series approximation for \(\epsilon\) varying from \(0.01\) to \(1\).

Find the \(\epsilon^2\) term in the McLaurin series for \(ds/d\lambda\) when \(\epsilon \approx 0\).