The historical roots of the calculus of variations trace back to Isaac Newton's Principia Mathematica problem concerning the shape of a solid of revolution that experiences minimal resistance to rapid motion through a "rare medium" consisting of elastic particles. Newton's stated solution to this problem barely hinted at the complexity and sophistication of his underlying analysis, and remained enigmatic to most readers of the Principia, at least until the late nineteenth century when rigorous foundations for the variational calculus began to be laid. This article exploits modern computer algebra to explore the meaning and origin of Newton's analysis and solution. This is a Mathematica notebook version of a paper that originally appeared in The Mathematica Journal, Volume 7 (Winter 1997), pages 64-71.
The study of motion of bodies in vacuo dates back at least to Galileo . The effect of air resistance was first analyzed by Newton in Book 2 of the Principia Mathematica [Newton 1687]. But his original analysis was so intricate and so sparsely revealed that much of what Newton said (or hinted at) remained only partly understood for well over two hundred years. Now the availability of a computer algebra system using traditional notation for input and output especially encourages an attempt to retrace his steps.
Figure 1 illustrates a solid of revolution consisting of a symmetrical "nose cone" attached to a cylindrical body, moving (downward) with velocity v through air. Alternatively, we can consider the solid to be stationary with the air streaming upward. Newton considered the air to consist of small elastic particles (such as molecules), each of mass m, uniformly distibuted with N particles per unit volume. The density of the air is then r = Nm and the flux of air particles across a unit area normal to the flow is Nv (particles per second per unit area). When one of these particles strikes the nose cone at angle f to the surface normal, in rebounding elastically it experiences a change 2mv cos(f) in its momentum. The resulting change in the vertical component of momentum is 2mv cos2(f). The rate at which such particles strike an area element dA of the nose cone surface S is Nv cos(f) dA, so the vertical force (rate of change of momentum, by Newton's own second law of motion) of air resistance experienced by this area element is the product dF = 2 r v2 cos3(f) dA. The total force of air resistance experienced by the solid is therefore given by the surface integral
evaluated over the surface S of the nose cone. We therefore see that (under Newton's assumptions) the air resistance is proportional to the density r of the air and to the square of the velocity v, with drag coefficient R given by
If the nose cone of the solid is a flat circle of radius r, not rounded at all, then
the area of the circle. Hence we may regard the drag coefficient R as the body's "effective cross-sectional area for resistance."
In the first three examples below, we will consider a solid with a cylindrical body with radius r = 1, having a nose cone of height h = 1 that is generated by revolving the curve y = y(x), 0 < x < 1 around the y-axis. Using cylindrical shells and the fact that
where ds = (dx2 + dy2)1/2 is the curve's arclength element, it is an elementary calculus exercise to show that the "reduced" drag coefficient Rp = R/p is given by the integral
In the computations that follow, we will for convenience delete the subscript and use R to denote this reduced drag coefficient.
Figure 2 shows a pointed cone defined by
The reduced drag coefficient R is given by
Figure 3 shows a rounded hemispherical nose cone defined by
The drag coefficient is
Thus the pointed conical nose and the rounded hemispherical nose provide the same air resistance, namely, half that of a cylinder with a flat circular "nose." This observation (for the hemispherical nose) was Newton's starting point in his Principia discussion of the solid of revolution of least resistance. (See Proposition 34 in Section VII of Book 2 in [Newton 1687]; this was Proposition 35 in the first edition of the Principia.)
Consider now a parabolic nose cone defined by
Then R is
Thus the parabolic nose cone offers less air resistance (about 20% less) than either the actual cone or the hemisphere.
Newton then raised the question of the shape of the optimal nose cone that provides the least possible air resistance. This was the first historical instance of a "calculus of variations" problem (predating the brachistochrone problem of the 1690s). He had the prescience to consider the possibility of a flat-tipped nose cone.
Figure 4 shows a "conical frustum" nose cone with a flat tip of radius a and with a denoting the angle of inclination of its side from the vertical. If the conical frustum has height h and top radius b, then
and cos f = sin a. Adding the reduced drag coefficient a2 of the flat tip and the integral corresponding to the curved frustum, we find that the reduced drag coefficient is given by
At this point Newton asked for the optimal conical frustum of height h and top radius b. Differentiating with respect to a,
we therefore see that the air resistance is minimal when tan(2a ) = 2b/h:
Thus a flat-tipped nose cone can offer less air resistance than either a pointed one or a smoothly rounded one!
Note that if the height h -> 0 while the base radius b remains fixed, then the relation
Thus an "infinitesimally thin" conical frustum nose cone has exterior base angle a = p/4. This result will play an important role in the following sections.
Ignoring the flat tip, the reduced drag coefficient of the curved surface itself is
But , where is the slant height of the conical frustum, so
This formula for the reduced drag coefficient R of the curved surface of a conical frustum of slant height s and base radii a and b is crucial to the investigation that follows.
After calculating the drag coefficient of the conical frustum, Newton reasoned that the optimal nose cone of least resistance would have a flat circular tip joined to the cylindrical body by a curvilinear band. What should be the radius a of the flat tip and what should be the shape y = y(x) of the arc generating this optimal band by revolution around the y-axis? Nowadays we would regard this as a calculus of variations problem (complicated by a variable endpoint condition) and proceed to set up the appropriate Euler-Lagrange equation.
Not having this Euler-Lagrange equation at his disposal, Newton proceeded directly, envisioning the curvilinear band as comprised of very narrow conical frusta (see his notes and preliminary analyses in [Newton 1974, 456--465]. Consider an adjacent pair of such conical frusta, obtained by revolution around the y-axis of the segment from P0 = (x-h, y+z-k) to P1 = (x, y) and the segment from P1 to P2 = (x+h, y+z+k) (Figure 5). Intuitively, Newton wanted to impose the condition that -- in the limit as h, k -> 0 with (x, y) fixed -- these two frusta offer miminal combined air resistance. Applying the conical frustum drag resistance formula, we find that their separate reduced drag coefficients are
Hence the total resistance coefficent of the two frusta is
We want to minimize R12 as a function of z with h, k, s, and x held constant, so we calculate the derivative with respect to z. To simplify the expression, we first divide by the factor 4h3.
Let's collect coefficients of powers of z in the numerator.
Now, reasoning (as Newton does) that z is very very small if both h and k are very small, let's "blot out" the higher powers of z,
set the result equal to 0, and solve for z.
We get our differential equation by noting from Figure 5 that when h is very small,
(thinking of the slope of the segment P0P2) and
(thinking of the second-degree term in the "Newton-Taylor series" which gives the departure z from the linear approximation).
So finally Newton's differential equation finally takes the "normal" form
How did Newton state this differential equation, given that Euler had not yet introduced the functional notation we are using here? In the Scholium to Proposition 34 in Book II of the Principia [Newton 1687, 333], we see a drawing equivalent to Figure 6, where BG designates the radius x0 of the nose cone's flat tip, and the curve GN (which we parametrize by y = y(x)) generates its curved surface. Newton asserts (without proof) that the shape of the optimal curved surface is determined by the fact that, given the fixed length BR, the point N where the tangent line is parallel to the segment GR is determined by the condition
If we write this condition in the form
and then substitute
the result is the first-order differential equation
(where C = BG/4 is constant). Then differentiation with respect to x yields precisely the second-order differential equation deq derived above.
Thus Newton's geometrically stated condition amounts to a first integral of deq. He asserts further that the base angle indicated in Figure 6 is , consistent with our result in the previous section.
Let us digress from Newton's line of investigation to explore a modern numerical solution of his differential equation for the surface of least resistance. (Though Newton did not have a computer algebra system at his disposal, he was no stranger to numerical integration, as one might infer from the nomenclature of "Newton-Cotes" methods.)
Suppose the cylindrical body has radius and height r = h = 1, and denote by x0 the (unknown) radius of the flat tip of the optimal nose cone. Then our solution y(x) should satisfy the endpoint conditions
Moreover, since the inclination of the curved surface at the tip is , we want to determine the initial point x0 so that the solution on [ x0,1] satisfying these endpoint conditions also satisfies the additional initial condition
y' (x0) = 1.
(This latter condition was the first historical instance of a "transversality condition" in a calculus of variations problem.) Thus our task is to find the solution of the initial value problem
such that y(1) = 1.
Let's first try x0 = 0.5 as an initial guess, solve our system using NDSolve, and then evaluate y(1).
From the differential equation with y' (x) > 0 we see that y'' (x) > 0, so the curve y = y(x) is convex upward, with its graph appearing as pictured in Figure 6. Evidently, we need to start "sooner" in order to end higher at x = 1. So let's try x0 = 0.3:
Now we've started too soon. The method of bisection would call for us to split the difference between x0 = 0.3 and x0 = 0.5, and try x0 = 0.4 next. Instead, let's first define the pertinent function f such that f (x0) = y(1), so we can apply FindRoot to find the desired value of x0.
First we check the two values calculated previously,
to verify that there is, indeed, a solution of the equation f (x0) = 1 between x0 = 0.3 and x0 = 0.5. We can use FindRoot to find this solution.
Thus it appears that the desired value -- the radius of the circular flat spot at the tip of the optimal nose cone -- is approximately
(Recall that we are investigating a nose-cone of fixed radius and height both equal to 1.)
Finally, we solve Newton's initial value problem with the value of x0 just found, and then calculate the corresponding minimal resistance coefficient.
Thus the minimum possible value of the reduced drag coefficient is about 75% that of a rounded hemispherical nose.
To retrace in modern notation the steps we believe Newton may have taken, let us write the differential equation,
Noting that the dependent variable y is missing from the equation, we begin with the usual reduction-of-order substitution u = y'.
Obviously we can separate the variables.
We'll let Mathematica do the integration by partial fractions.
Recalling our substitution, we get
Note that the dependent variable y is missing from this first-order equation which gives x explicitly in terms of y' (x) . We can solve in a standard fashion for x and y in terms of the parameter v = y'.
The usual chain rule device gives
Finally we want to impose the boundary conditions
The second of these says that x = x0 when v = 1, and then the first says that y = 0 when v = 1. To satisfy the left endpoint condition, we need to choose v1 so that x(v1) = y(v1) = 1. First, we find the values of C2 and C3.
With these substitutions, our two equations
define the desired solution curve parametrically in terms of the slope v = y' once x0 has been determined. We proceed to solve numerically, initially estimating x0 = 0.5 (something between 0 and 1) and v = 2 (something larger than 1) for the value of the parameter where x = y = 1.
So we pick out the initial x-value (tip radius) and the final v-value (slope)
and plot the axial profile shape of the optimal nose cone.
Finally, let's plot a 3-dimensional picture of Newton's optimal nose cone of least resistance, flat tip and all.
On what basis can we plausibly hope to have retraced -- using a computer algebra system -- the steps that Newton himself took? As Whiteside remarks (see page 466 of ), "The immediate reaction of Newton's contemporaries to this scholium [on the solid of least resistance] on its publication in the 1687 Principia was one of near-total incomprehension." Sometime in the mid-1690s, Newton composed an intended addendum to clarify the matter (though it never appeared in any subsequent edition of the Principia). Figure 9 is an elaboration of Figure 6 above and corresponds closely to Newton's figure as reproduced by Whiteside [Newton 1974, 479] (though with the x-axis now vertical and the y-axis points to the left; the positive z-axis is the negative y-axis). Given the radius x0 = BG of the nose-cone's flat circular tip and a fixed segment BE that serves to define the ratio v = BE/BG, the question (as explained in our discussion of Figure 6) is how to determine the point N on the nose-cone curve where the slope of the curve is y' = v.
Newton provides the following steps for the geometric construction of the point N. We translate Whiteside's Latin transcription of the hand-written manuscript, and use Mathematica to evaluate the results as we proceed through the construction.
1. Noting that
pick the point S on the x-axis such that
Then x = BS will be the x-coordinate of the desired point N. Note that this construction agrees with the expression for x(v) given in our evaluation of eq1 above. The y-coordinate of N is determined by the remaining steps.
2. Pick the point F on the x-axis such that
3. Pick the point A on the positive z-axis (the negative y-axis) such that BA = BG = x0. Then
4. Erect at A the perpendicular AH of length AH = BG + FS, that is,
and then complete the rectangle HAEI shown in Figure 9.
5. Construct the rectangular hyperbola xz = c that is tangent to the segment AG. An elementary computation shows that the hyperbola will be tangent to the hypotenuse AG of the triangle ABG provided that
Then denote by K and L the indicated points of intersection of this hyperbola with the vertical sides of the rectangle HAEI.
6. Finally, erect the segment NS perpendicular to the x-axis with length y = NS given by
The curvilinear area in the numerator on the right-hand side is given by
so the y-coordinate of N is
When we compare this with our previous result
se see that Newton's geometric construction (which he sets forth without any indication of its origins) duplicates the parametrization that we have (and surely he) derived by solution of his differential equation.
In this article we have used formidable computational machinery to explicate just a few pages of Newton's published work and preliminary notes on the solid of least resistance. Paraphrasing Whiteside, we can hardly fail to note how little justice the published Principia Mathematica does to the complexity and sophistication of the preliminary analyses on which that founding document of modern exact science was based, especially when we consider the geometric machinery he himself employed [Newton 1974, 457]. As William Whewell wrote in 1837 (Whewell 1837, 167],
"The ponderous instrument of synthesis, so effective in his [Newton's] hands, has never since been grasped by [any]one who could use it for such purposes; and we gaze at it with admiring curiosity, as on some gigantic implement of war, which stands idle among the memorials of ancient days, and makes us wonder what manner of man he was who could wield as a weapon what we can hardly lift as a burden."
Galileo Galilei. 1638. Dialogues Concerning the Two New Sciences. Translated by Henry Crew and Alfonso de Salvio, Great Books of the Western World vol. 28. Encyclopedia Britannica, Inc., 1952.
Newton, Isaac. 1687. Philosophiae Naturalis Principia
Mathematica (Mathematical Principles of Natural Philosophy).
Translated by Andrew Motte, revised by Florian Cajori. University of
California Press, 1966.
Newton, Isaac. 1974. The Mathematical Papers of Isaac Newton Volume VI. Edited by D. T. Whiteside. Cambridge University Press.
Whewell, William. 1837. History of the Inductive Sciences, from the Earliest to the Present Time. London.