Philosophiæ Naturalis
Principia Mathematica
Physics · Mathematics · Proofs · Computations
Sir Isaac Newton · London, 1687 · Third Edition 1726
I. Introduction & Historical Context
The Principia Mathematica — published 5 July 1687 — is the single most influential work in the history of science. In three books Newton established the mathematical framework for classical mechanics, derived the universal law of gravitation, and explained the motions of planets, comets, moons, the tides, and the precession of the equinoxes from first principles.
Newton developed his mathematics in the period 1665–1666 while Cambridge was closed due to plague, calling it his annus mirabilis. He referred to his calculus as the Method of Fluxions, and used a purely geometric exposition in the Principia, believing (correctly) that geometric proofs were more rigorous to his contemporaries than his new analytic methods.
| Book | Latin Title | Topic |
|---|---|---|
| I | De Motu Corporum I | Motion in non-resisting media; orbital theory; inverse-square force law |
| II | De Motu Corporum II | Motion in resisting media; fluid mechanics; wave propagation |
| III | De Mundi Systemate | System of the World — planets, moons, comets, tides, precession |
II. The Eight Definitions
Newton opens with eight fundamental definitions. These establish the ontology of his mechanics.
The quantity of matter (mass, \(m\)) is the measure of it arising conjointly from its density and bulk (volume \(V\)).
$$m = \rho V$$Mass is an intrinsic, invariant property of a body — distinct from weight, which varies with location.
The quantity of motion is the product of mass and velocity:
$$\mathbf{p} = m\mathbf{v}$$This is the modern linear momentum. Newton's great insight was that force changes momentum, not velocity directly.
The innate force of matter (\(vis\ insita\)) is a power of resisting by which every body, as much as in it lies, endeavours to persevere in its present state — rest or uniform rectilinear motion.
$$\text{If } \mathbf{F} = \mathbf{0},\quad \frac{d\mathbf{v}}{dt} = \mathbf{0}$$An impressed force is an action exerted upon a body to change its state. It exists only in the action; the body persists in its new state by inertia alone.
A centripetal force is that by which bodies are drawn or impelled towards some point as to a centre. Examples: gravity (towards Earth's centre), magnetism, the string of a sling.
$$\mathbf{F}_c = -\frac{mv^2}{r}\hat{r} = -m\omega^2 r\,\hat{r}$$Newton distinguishes three measures of centripetal force:
- Absolute quantity — analogous to the strength of a magnet; the source strength (e.g., \(GM\)).
- Accelerative quantity — acceleration it produces: \(g = GM/r^2\).
- Motive quantity — the actual force on the body: \(F = mg = GMm/r^2\).
III. The Three Laws of Motion
"Every body perseveres in its state of rest, or of uniform motion in a right line, unless it is compelled to change that state by forces impressed thereon."
$$\mathbf{F} = \mathbf{0} \implies \frac{d\mathbf{v}}{dt} = \mathbf{0} \implies \mathbf{v} = \text{const}$$"The alteration of motion is ever proportional to the motive force impressed; and is made in the direction of the right line in which that force is impressed."
$$\mathbf{F} = \frac{d\mathbf{p}}{dt} = \frac{d(m\mathbf{v})}{dt}$$For constant mass \(m\):
$$\boxed{\mathbf{F} = m\mathbf{a}}$$Note: Newton stated the law in terms of impulse and change of momentum, not instantaneous force — the modern form \(\mathbf{F}=m\mathbf{a}\) emerges from taking the continuous limit.
"To every action there is always opposed an equal reaction; or the mutual actions of two bodies upon each other are always equal, and directed to contrary parts."
$$\mathbf{F}_{12} = -\mathbf{F}_{21}$$This implies conservation of total momentum in an isolated system: if \(\mathbf{F}_{12}+\mathbf{F}_{21}=0\), then \(\frac{d(\mathbf{p}_1+\mathbf{p}_2)}{dt}=0\).
Corollaries to the Laws
A body acted upon by two forces simultaneously will describe the diagonal of a parallelogram in the same time it would have described the sides by those forces separately.
$$\mathbf{F}_\text{net} = \mathbf{F}_1 + \mathbf{F}_2$$The common centre of gravity of two or more bodies does not alter its state of motion or rest by the actions of the bodies among themselves — proto-conservation of momentum.
$$\mathbf{p}_\text{total} = \sum_i m_i \mathbf{v}_i = \text{const when }\sum \mathbf{F}_\text{ext}=0$$IV. The Method of Fluxions (Calculus)
Newton called time-derivatives fluxions and quantities that flow fluents. His notation used dots (still used in mechanics today): \(\dot{x} = dx/dt\), \(\ddot{x} = d^2x/dt^2\).
Fluxions and Modern Notation
| Newton's Term | Notation | Modern Meaning |
|---|---|---|
| Fluent | \(x\) | Function of time |
| Fluxion | \(\dot{x}\) | \(dx/dt\) — rate of change |
| Second fluxion | \(\ddot{x}\) | \(d^2x/dt^2\) — acceleration |
| Moment | \(o\cdot\dot{x}\) | \(dx\) — infinitesimal increment |
| Genita | products, powers | Functions formed algebraically |
Newton's Rule for Differentiation
Newton derived what we now call the power rule. Let \(y = x^n\). Replacing \(x\) by \(x + o\cdot\dot{x}\):
$$y + o\cdot\dot{y} = (x + o\dot{x})^n = x^n + n x^{n-1}(o\dot{x}) + \tfrac{n(n-1)}{2}x^{n-2}(o\dot{x})^2 + \cdots$$Subtracting \(y = x^n\), dividing by \(o\), and discarding higher-order terms in \(o\):
$$\boxed{\dot{y} = n x^{n-1}\dot{x}} \quad\Longleftrightarrow\quad \frac{d(x^n)}{dt} = n x^{n-1}\frac{dx}{dt}$$Let \(z = xy\). Then \(z + o\dot{z} = (x+o\dot{x})(y+o\dot{y}) = xy + o(x\dot{y}+y\dot{x}) + o^2\dot{x}\dot{y}\). Subtracting \(z=xy\), dividing by \(o\), neglecting \(o^2\):
$$\frac{d(xy)}{dt} = x\frac{dy}{dt} + y\frac{dx}{dt}$$Fundamental Theorem of Calculus (Newton's Version)
Newton stated: if \(y\) is a fluent and \(\dot{y}=f(x)\dot{x}\), then the area under the curve \(f\) is \(y\). The area function \(A(x) = \int_a^x f(t)\,dt\) satisfies \(\frac{dA}{dx}=f(x)\).
$$\frac{d}{dx}\int_a^x f(t)\,dt = f(x)$$V. Method of First and Last Ratios (Limits)
In Lemmas I–XI, Newton establishes the theory of limits geometrically — predating Cauchy's epsilon-delta formalization by 150 years.
"Quantities, and the ratios of quantities, which in any finite time converge continually to equality, and before the end of that time approach nearer to each other than by any given difference, are ultimately equal."
This is Newton's definition of a limit: \(\lim_{t\to T} f(t) = L\) iff for every \(\varepsilon>0\) there exists \(\delta\) such that \(|f(t)-L|<\varepsilon\) when \(|t-T|<\delta\).
The area under a curve \(AB\) is the limit of the sum of inscribed (or circumscribed) rectangles as the number of rectangles \(n\to\infty\) and their width \(\to 0\).
$$A = \lim_{n\to\infty}\sum_{k=1}^n f(x_k)\Delta x = \int_a^b f(x)\,dx$$The ultimate ratio of an arc, its chord, and the tangent is one of equality. Let arc \(AB\), chord \(AB\), and tangent \(AT\) approach zero together. As \(B\to A\):
$$\lim_{B\to A}\frac{\text{arc }AB}{\text{chord }AB} = 1,\quad \lim_{B\to A}\frac{\text{chord }AB}{\text{tangent }AT} = 1$$This justifies replacing arc-length with tangent-length for infinitesimal motions, central to the orbital mechanics proofs.
The Subleme Ratio — Instantaneous Velocity
From Lemma VII, the instantaneous velocity at \(t\) is:
$$v(t) = \lim_{\Delta t\to 0}\frac{\Delta x}{\Delta t} = \dot{x}$$Newton calls this the "first ratio of nascent quantities" or the "last ratio of evanescent quantities."
VI. Infinite Series & Binomial Theorem
Newton's 1665 discovery of the generalized binomial series was a key tool. For any real exponent \(n\) and \(|x|<1\):
$$\boxed{(1+x)^n = \sum_{k=0}^\infty \binom{n}{k}x^k = 1 + nx + \frac{n(n-1)}{2!}x^2 + \frac{n(n-1)(n-2)}{3!}x^3 + \cdots}$$Newton used this to expand functions needed in orbital calculations. For example, for the gravitational potential with \(r\gg d\):
$$\frac{1}{r+d} = \frac{1}{r}\left(1+\frac{d}{r}\right)^{-1} = \frac{1}{r}\left(1 - \frac{d}{r} + \frac{d^2}{r^2} - \cdots\right)$$Newton's Series for sin and cos
Newton derived these from the binomial series and his integration of \((1-x^2)^{-1/2}\):
$$\sin x = x - \frac{x^3}{3!} + \frac{x^5}{5!} - \cdots,\qquad \cos x = 1 - \frac{x^2}{2!} + \frac{x^4}{4!} - \cdots$$VII. Vector Kinematics & Polar Coordinates
For orbital mechanics, polar coordinates \((r,\theta)\) are natural. Position in polar form:
$$\mathbf{r} = r\hat{r}$$Since \(\hat{r}\) rotates with the particle, differentiating requires chain rule over \(\theta\). Define unit vectors:
$$\hat{r} = \cos\theta\,\hat{i}+\sin\theta\,\hat{j},\qquad \hat{\theta} = -\sin\theta\,\hat{i}+\cos\theta\,\hat{j}$$Key identities: \(\dot{\hat{r}} = \dot{\theta}\hat{\theta}\) and \(\dot{\hat{\theta}} = -\dot{\theta}\hat{r}\).
Velocity and Acceleration in Polar Coordinates
The \(\hat{r}\)-component gives the radial equation; the \(\hat{\theta}\)-component gives the transverse equation. For a central force \(\mathbf{F}=F(r)\hat{r}\):
$$\text{Transverse: }\quad r\ddot{\theta}+2\dot{r}\dot{\theta} = \frac{1}{r}\frac{d}{dt}(r^2\dot{\theta}) = 0 \implies r^2\dot{\theta} = h = \text{const}$$This constant \(h\) is the specific angular momentum — and its constancy is Kepler's Second Law.
VIII. Prop. I — Kepler's Second Law (Equal Areas)
"The areas which revolving bodies describe by radii drawn to an immovable centre of force lie in the same immovable planes, and are proportional to the times of their description."
This is a purely kinematic consequence of any central force — it does not require gravity specifically.
Setup: Divide time into equal intervals \(\Delta t\). During the first interval, with no force, the body moves from \(A\) to \(B\). The triangle \(SAB\) (S = centre) has area \(\tfrac{1}{2}|SA||AB|\sin\angle SAB\).
Impulse step: At \(B\), a central impulse \(\mathbf{J}\) acts (towards \(S\)). Without it, the body would reach \(c\) (so \(Bc = AB\)). With it, the body reaches \(C = c + \delta c\), where \(\delta c\) is parallel to \(BS\) (the force is central).
Area comparison: Triangles \(SBc\) and \(SBC\) share base \(SB\) and their vertices \(c,C\) differ only by \(\delta c \parallel SB\), so they have equal area. Thus \(\text{Area}(SAB) = \text{Area}(SBC)\).
Limit: Passing to the continuous limit (\(\Delta t\to 0\)), the swept area per unit time is constant:
$$\frac{dA}{dt} = \frac{1}{2}r^2\dot{\theta} = \frac{h}{2} = \text{const}$$Modern Proof via Angular Momentum
The areal velocity is:
$$\frac{dA}{dt} = \frac{1}{2}\left|\mathbf{r}\times\mathbf{v}\right| = \frac{L}{2m}$$where \(L = m|\mathbf{r}\times\mathbf{v}|\) is angular momentum. Differentiating:
$$\frac{dL}{dt} = \mathbf{r}\times\mathbf{F}$$For a central force \(\mathbf{F} = f(r)\hat{r}\), the torque \(\mathbf{r}\times\mathbf{F} = r\hat{r}\times f(r)\hat{r} = \mathbf{0}\). Hence \(L = \text{const}\) and \(\frac{dA}{dt} = \text{const}\).
IX. Prop. VI — Centripetal Force Formula
For a body moving on a curve with a central force directed to pole \(S\), the magnitude of the centripetal acceleration is:
$$\boxed{a_c = \frac{v^2 \cdot QR}{SP^2 \cdot QT^2} \longrightarrow \frac{h^2}{r^2 \cdot p}}$$where \(p\) is the perpendicular distance from \(S\) to the tangent, \(h = r^2\dot\theta\) is specific angular momentum.
In modern terms, for circular-ish motion with curvature radius \(\rho\):
$$F = \frac{mv^2}{\rho}$$For a general central force orbit in polar form \(r=r(\theta)\), setting \(u=1/r\):
$$F(r) = \frac{mh^2 u^2}{r^2}\left(u + \frac{d^2u}{d\theta^2}\right) = \frac{mh^2}{r^2}\left(\frac{1}{r} + \frac{d^2(1/r)}{d\theta^2}\right)$$This is the Binet equation, central to Prop. XI.
X. Prop. XI — Ellipse Implies Inverse Square (The Crown Jewel)
"If a body revolves in an ellipse, it is required to find the law of the centripetal force tending to the focus of the ellipse."
This is arguably Newton's greatest mathematical achievement — deriving that any elliptical orbit around a focus requires an inverse-square force.
Setup: Geometry of the Ellipse
An ellipse with semi-major axis \(a\), semi-minor axis \(b\), eccentricity \(e\), and focus at origin \(S\):
$$r = \frac{\ell}{1+e\cos\theta}, \qquad \ell = \frac{b^2}{a} = a(1-e^2)$$where \(\ell\) is the semi-latus rectum.
Let \(u = 1/r\). From the ellipse equation \(r = \ell/(1+e\cos\theta)\):
$$u = \frac{1}{r} = \frac{1+e\cos\theta}{\ell}$$Differentiate twice with respect to \(\theta\):
$$\frac{du}{d\theta} = -\frac{e\sin\theta}{\ell},\qquad \frac{d^2u}{d\theta^2} = -\frac{e\cos\theta}{\ell}$$Now substitute into the Binet equation \(F = mh^2 u^2(u + d^2u/d\theta^2)\):
$$u + \frac{d^2u}{d\theta^2} = \frac{1+e\cos\theta}{\ell} - \frac{e\cos\theta}{\ell} = \frac{1}{\ell}$$Therefore:
$$F = mh^2 u^2 \cdot \frac{1}{\ell} = \frac{mh^2}{\ell}\cdot\frac{1}{r^2}$$Since \(mh^2/\ell = GM\cdot m\) (shown by matching to circular orbit), we obtain:
$$\boxed{F = \frac{GMm}{r^2}}$$The force law is precisely inverse square in distance — and the constant of proportionality involves only the masses and Newton's constant.
Converse (Prop. XIII): Inverse Square implies Conics
Newton also proved the converse: an inverse-square centripetal force produces conic section orbits. The shape depends on the total energy \(E = \frac{1}{2}mv^2 - \frac{GMm}{r}\):
| Energy \(E\) | Orbit Type | Eccentricity |
|---|---|---|
| \(E < 0\) | Ellipse (includes circle) | \(0 \le e < 1\) |
| \(E = 0\) | Parabola | \(e = 1\) |
| \(E > 0\) | Hyperbola | \(e > 1\) |
XI. Conic Sections — Complete Treatment
Polar Equations (Focus at Origin)
$$r(\theta) = \frac{\ell}{1 + e\cos\theta}, \qquad \ell = \frac{h^2}{GM}$$Orbital Parameters
| Parameter | Symbol | Formula |
|---|---|---|
| Semi-major axis | \(a\) | \(a = \ell/(1-e^2)\) (ellipse) |
| Semi-minor axis | \(b\) | \(b = a\sqrt{1-e^2}\) |
| Semi-latus rectum | \(\ell\) | \(\ell = b^2/a = h^2/(GM)\) |
| Periapsis distance | \(r_p\) | \(r_p = \ell/(1+e) = a(1-e)\) |
| Apoapsis distance | \(r_a\) | \(r_a = \ell/(1-e) = a(1+e)\) |
| Orbital period | \(T\) | \(T = 2\pi a^{3/2}/\sqrt{GM}\) |
| Orbital energy | \(E\) | \(E = -GMm/(2a)\) |
Vis-Viva Equation
The total orbital energy is conserved:
$$E = \frac{1}{2}mv^2 - \frac{GMm}{r} = -\frac{GMm}{2a}$$ $$\boxed{v^2 = GM\left(\frac{2}{r} - \frac{1}{a}\right)}$$This relates orbital speed to position at any point on the orbit, without solving differential equations.
Kepler's Equation (Time–Angle Relation)
The true anomaly \(\theta\) relates to the eccentric anomaly \(E_a\) and mean anomaly \(M = n(t-t_0)\) (where \(n=2\pi/T\)) by:
$$M = E_a - e\sin E_a \qquad\text{(Kepler's Equation)}$$ $$\tan\frac{\theta}{2} = \sqrt{\frac{1+e}{1-e}}\tan\frac{E_a}{2}$$Kepler's Equation must be solved numerically or by series inversion.
XII. The Shell Theorem (Prop. LXXI & LXXII)
Prop. LXXI: A spherical shell of uniform density attracts an external point mass as if all its mass were concentrated at its centre.
Prop. LXXII: A spherical shell exerts no gravitational force on an internal point mass (complete cancellation).
Consider a thin shell of radius \(R\), surface density \(\sigma\), total mass \(M = 4\pi R^2\sigma\). For a point \(P\) at distance \(d\) from the centre, integrate over all mass elements \(dm = \sigma\,dA\).
By symmetry, only the component along \(PO\) survives. Parametrize by angle \(\psi\) from the axis:
$$F = \int_0^\pi \frac{G\,dm}{s^2}\cos\alpha$$where \(s\) is the distance from \(P\) to the ring, and \(\alpha\) is the angle the force makes with axis \(PO\).
Using the law of cosines: \(s^2 = R^2+d^2-2Rd\cos\psi\), and geometry:
$$dm = 2\pi R^2\sigma\sin\psi\,d\psi,\qquad \cos\alpha = \frac{d-R\cos\psi}{s}$$ $$F = \frac{G\sigma \cdot 2\pi R^2}{2}\int_{-1}^{1}\frac{(d-Ru)}{(R^2+d^2-2Rdu)^{3/2}}\,du$$External case \(d>R\): The integral evaluates to \(\frac{2}{d^2}\), giving \(F = \frac{GM}{d^2}\). ✓
Internal case \(d
Consequently, for a uniform sphere of radius \(R_E\) and mass \(M_E\):
$$g(r) = \begin{cases} \dfrac{GM_E}{r^2} & r \ge R_E \\[6pt] \dfrac{GM_E\, r}{R_E^3} & r < R_E \end{cases}$$XIII. Law of Universal Gravitation
Every particle of matter in the universe attracts every other particle with a force proportional to the product of their masses and inversely proportional to the square of the distance between them.
\(G = 6.674\times10^{-11}\ \text{N\,m}^2/\text{kg}^2\) (Newton's constant, determined experimentally).
Newton's Derivation Strategy
- Show all moons and planets sweep equal areas (Kepler II holds → central force).
- Show orbital periods obey \(T^2\propto a^3\) (Kepler III).
- From \(T^2\propto a^3\) for circular orbits: \(F\propto 1/r^2\).
- From Prop. XI: any elliptical orbit requires \(F\propto 1/r^2\).
- Apply Shell Theorem: spherical bodies act as point masses.
- The Moon test: the same \(F\propto 1/r^2\) law that holds here extends to the Moon.
Gravitational Potential
$$\Phi(r) = -\frac{GM}{r}$$ $$\mathbf{F} = -m\nabla\Phi = -\frac{GMm}{r^2}\hat{r}$$The potential energy of a two-body system:
$$U(r) = m\Phi(r) = -\frac{GMm}{r}$$The negative sign means it takes energy to separate two masses — they are in a bound state when \(E=T+U<0\).
XIV. Kepler's Three Laws — Derived from Newton
Kepler's First Law: Orbits are Conics
The radial equation of motion for a central \(1/r^2\) force is:
$$\ddot{r} - r\dot{\theta}^2 = -\frac{GM}{r^2}$$Substitute \(u=1/r\), use \(\dot{\theta}=h/r^2=hu^2\) and \(\dot{r}=-h\,du/d\theta\):
$$\ddot{r} = -h^2u^2\frac{d^2u}{d\theta^2}$$The radial equation becomes:
$$-h^2u^2\frac{d^2u}{d\theta^2} - h^2u^3 = -GMu^2$$ $$\frac{d^2u}{d\theta^2} + u = \frac{GM}{h^2}$$This is a linear ODE. The general solution is:
$$u = \frac{GM}{h^2}(1 + e\cos(\theta-\theta_0))$$ $$\boxed{r = \frac{h^2/GM}{1+e\cos\theta}}$$This is the polar equation of a conic section with eccentricity \(e\) determined by initial conditions. ✓
Kepler's Second Law: Equal Areas
Already proved above as Prop. I — it follows from angular momentum conservation for any central force.
Kepler's Third Law: \(T^2 \propto a^3\)
The area of an ellipse is \(A = \pi ab\). Since \(dA/dt = h/2\):
$$T = \frac{A}{h/2} = \frac{2\pi ab}{h}$$Using \(b = a\sqrt{1-e^2}\) and \(h = \sqrt{GM\ell} = \sqrt{GMa(1-e^2)}\):
$$T = \frac{2\pi a \cdot a\sqrt{1-e^2}}{\sqrt{GMa(1-e^2)}} = \frac{2\pi a^{3/2}}{\sqrt{GM}}$$ $$\boxed{T^2 = \frac{4\pi^2}{GM}\,a^3}$$The constant of proportionality depends only on \(GM\) of the central body — same for all satellites.
| Planet | \(a\) (AU) | \(T\) (yr) | \(T^2/a^3\) |
|---|---|---|---|
| Mercury | 0.387 | 0.241 | 1.002 |
| Venus | 0.723 | 0.615 | 1.001 |
| Earth | 1.000 | 1.000 | 1.000 |
| Mars | 1.524 | 1.881 | 0.999 |
| Jupiter | 5.203 | 11.86 | 1.000 |
| Saturn | 9.537 | 29.46 | 1.000 |
XV. Orbital Mechanics
Orbital Velocity
From the vis-viva equation:
| Orbit | Speed at radius \(r\) |
|---|---|
| Circular (\(a=r\)) | \(v_c = \sqrt{GM/r}\) |
| Elliptical | \(v = \sqrt{GM(2/r-1/a)}\) |
| At periapsis | \(v_p = \sqrt{GM(1+e)/(a(1-e))}\) |
| At apoapsis | \(v_a = \sqrt{GM(1-e)/(a(1+e))}\) |
Hohmann Transfer Orbit
To transfer from circular orbit \(r_1\) to \(r_2 > r_1\), use an ellipse with \(a=(r_1+r_2)/2\). Two impulses:
$$\Delta v_1 = \sqrt{\frac{GM}{r_1}}\left(\sqrt{\frac{2r_2}{r_1+r_2}}-1\right), \quad \Delta v_2 = \sqrt{\frac{GM}{r_2}}\left(1-\sqrt{\frac{2r_1}{r_1+r_2}}\right)$$Conservation Laws Summary
| Conserved Quantity | Expression | Physical Origin |
|---|---|---|
| Angular momentum | \(L = mr^2\dot{\theta}\) | Central force (no torque) |
| Energy | \(E = \frac{1}{2}mv^2 - \frac{GMm}{r}\) | Conservative force |
| Laplace-Runge-Lenz | \(\mathbf{A} = \mathbf{v}\times\mathbf{L} - GMm\hat{r}\) | Inverse-square only |
XVI. Escape Velocity
For a body to escape to infinity (\(r\to\infty\), \(v\to 0\)), total energy must be \(E\ge 0\):
$$\frac{1}{2}mv_e^2 - \frac{GMm}{R} = 0$$ $$\boxed{v_e = \sqrt{\frac{2GM}{R}}}$$Note \(v_e = \sqrt{2}\,v_c\) where \(v_c = \sqrt{GM/R}\) is the circular orbital speed. For Earth:
$$v_e = \sqrt{\frac{2\times6.674\times10^{-11}\times5.972\times10^{24}}{6.371\times10^6}} \approx 11.19\ \text{km/s}$$XVII. Book II — Motion in Resisting Media
Newton considered drag forces proportional to powers of velocity. The two simplest cases:
Linear Drag (Stokes' Law — viscous regime)
$$m\dot{v} = mg - bv \qquad (b = 6\pi\eta R \text{ for a sphere})$$Let \(v_T = mg/b\) (terminal velocity). Then \(\dot{v} = g(1-v/v_T)\):
$$\int_0^v \frac{dv'}{1-v'/v_T} = gt \implies v(t) = v_T(1-e^{-gt/v_T})$$Position:
$$y(t) = v_T t + \frac{v_T^2}{g}(e^{-gt/v_T}-1)$$Quadratic Drag (Newton — turbulent regime)
$$m\dot{v} = mg - cv^2 \qquad (c = \tfrac{1}{2}C_d\rho A)$$ $$v(t) = v_T\tanh\left(\frac{gt}{v_T}\right), \qquad v_T = \sqrt{\frac{mg}{c}}$$Projectile with Drag (2D)
$$m\ddot{x} = -c|\mathbf{v}|v_x,\qquad m\ddot{y} = -mg - c|\mathbf{v}|v_y$$These coupled ODEs require numerical integration.
XVIII. The Pendulum (Book II & III)
Newton used pendulum experiments extensively to measure \(g\), test drag laws, and verify the equivalence of inertial and gravitational mass.
Equation of Motion
For a simple pendulum of length \(L\), angle \(\theta\) from vertical:
$$\ddot{\theta} + \frac{g}{L}\sin\theta = 0$$Small-Angle Approximation
For \(\theta\ll 1\), \(\sin\theta\approx\theta\):
$$\ddot{\theta} + \omega_0^2\theta = 0, \qquad \omega_0 = \sqrt{g/L}$$ $$\theta(t) = \theta_0\cos(\omega_0 t + \phi), \qquad T = 2\pi\sqrt{L/g}$$Period Correction for Finite Amplitude
The exact period requires an elliptic integral. Newton used a series expansion:
$$T = 2\pi\sqrt{\frac{L}{g}}\left(1 + \frac{1}{16}\theta_0^2 + \frac{11}{3072}\theta_0^4 + \cdots\right)$$Multiply \(\ddot\theta + \omega_0^2\sin\theta = 0\) by \(\dot\theta\) and integrate:
$$\frac{1}{2}\dot\theta^2 - \omega_0^2\cos\theta = \text{const} = -\omega_0^2\cos\theta_0$$ $$\dot\theta = \omega_0\sqrt{2(\cos\theta-\cos\theta_0)}$$ $$T = 4\int_0^{\theta_0}\frac{d\theta}{\omega_0\sqrt{2(\cos\theta-\cos\theta_0)}} = \frac{4}{\omega_0}\int_0^{\theta_0}\frac{d\theta}{\sqrt{4\sin^2(\theta_0/2)-4\sin^2(\theta/2)}}$$Setting \(\sin(\theta/2) = k\sin\phi\), \(k=\sin(\theta_0/2)\):
$$T = \frac{2\pi}{\omega_0}\cdot\frac{2}{\pi}K(k) = 2\pi\sqrt{L/g}\cdot\frac{2}{\pi}K(\sin(\theta_0/2))$$where \(K\) is the complete elliptic integral of the first kind.
XIX. Book III — Tides & Precession of the Equinoxes
Tidal Forces
The tidal force arises from the differential gradient of gravity across a body of size \(d\ll r\):
$$\Delta F = F(r-d) - F(r) \approx \frac{dF}{dr}\cdot d = -\frac{2GMm\,d}{r^3}$$Newton gave the first quantitative theory of ocean tides, explaining the semi-diurnal pattern and the relative contributions of the Moon and Sun:
$$\frac{F_\text{Moon}}{F_\text{Sun}} \approx \frac{M_\text{Moon}/r_\text{Moon}^3}{M_\text{Sun}/r_\text{Sun}^3} \approx 2.2$$Precession of the Equinoxes
Newton showed the Earth's oblate shape (equatorial bulge) causes the Moon and Sun to exert a torque, producing lunisolar precession with period ~26,000 years:
$$\frac{d\psi}{dt} = -\frac{3}{2}\frac{n^2}{Ω}\cos\varepsilon\cdot\frac{C-A}{C}$$where \(n\) = Moon's mean motion, \(\varepsilon\) = obliquity, \(C,A\) = moments of inertia.
XX. The Moon Test
Newton's celebrated confirmation that terrestrial and celestial gravity are the same force:
Moon's acceleration from orbital data:
$$a_M = \frac{4\pi^2 r_M}{T_M^2} = \frac{4\pi^2(60.27\times6.371\times10^6)}{(27.32\times86400)^2} \approx 2.72\times10^{-3}\ \text{m/s}^2$$Prediction from inverse square:
$$g_E/60.27^2 = 9.81/3632.5 \approx 2.70\times10^{-3}\ \text{m/s}^2$$Agreement to within 0.7% — conclusive confirmation! Newton writes: "the force by which the Moon is retained in its orbit is that very force which we commonly call gravity."
XXI. GNU Octave Examples
1 — Kinematics: Free Fall & Parabolic Trajectory
%% Projectile motion — Newton's Second Law, F = mg downward %% Demonstrates kinematic equations derived from integration of a = -g ĵ clear; clc; %% ─── Constants ─────────────────────────────────────────────── g = 9.81; % gravitational acceleration (m/s^2) v0 = 30; % initial speed (m/s) angles = [15 30 45 60 75]; % launch angles (degrees) colors = {'b', 'g', 'r', 'm', 'c'}; figure(1); hold on; grid on; title('Projectile Motion — Newton Principia Style', 'FontSize', 14); xlabel('Range x (m)'); ylabel('Height y (m)'); for k = 1:length(angles) alpha = angles(k) * pi / 180; vx = v0 * cos(alpha); vy = v0 * sin(alpha); %% Analytic solution: integrate a_x=0, a_y=-g twice T_flight = 2 * vy / g; % time of flight t = linspace(0, T_flight, 500); x = vx .* t; % x(t) = vx·t y = vy .* t - 0.5 .* g .* t.^2; % y(t) = vy·t − ½g·t² R = v0^2 * sin(2*alpha) / g; % range formula H = vy^2 / (2*g); % max height plot(x, y, colors{k}, 'LineWidth', 1.8, ... 'DisplayName', sprintf('θ=%d° R=%.1fm H=%.1fm', angles(k), R, H)); end legend('Location', 'northwest'); axis([0 100 0 50]); %% ─── Numerical Euler integration (verify analytic) ────────── alpha = 45 * pi / 180; dt = 0.001; state = [0; 0; v0*cos(alpha); v0*sin(alpha)]; % [x y vx vy] trajectory = state'; while state(2) >= 0 ax = 0; ay = -g; state(1) += state(3) * dt; state(2) += state(4) * dt; state(3) += ax * dt; state(4) += ay * dt; trajectory = [trajectory; state']; end plot(trajectory(:,1), trajectory(:,2), 'k--', ... 'LineWidth', 2, 'DisplayName', 'Euler 45° (numerical)'); printf('Range at 45°: analytic=%.3f, numerical=%.3f m\n', ... v0^2/g, trajectory(end,1));
2 — Orbit Integration (RK4)
%% Numerical integration of Newton's equations of motion for two-body problem %% F = -GM/r² r̂ → x'' = -GM·x/r³, y'' = -GM·y/r³ %% Uses RK4 to verify Kepler's laws numerically clear; clc; GM = 1.0; % gravitational parameter (normalised) e = 0.5; % eccentricity a = 1.0; % semi-major axis l = a*(1-e^2); % semi-latus rectum %% Initial conditions at periapsis (θ=0) r0 = a*(1-e); v0 = sqrt(GM*(1+e)/(a*(1-e))); state0 = [r0; 0; 0; v0]; % [x; y; vx; vy] %% Derivatives function (Newton's 2nd Law) function ds = deriv(s, GM) r3 = (s(1)^2 + s(2)^2)^(3/2); ds = [s(3); s(4); -GM*s(1)/r3; -GM*s(2)/r3]; end %% RK4 integrator T = 2*pi*a^(3/2)/sqrt(GM); % Kepler period dt = T/10000; N = round(T/dt); traj = zeros(N,4); s = state0; for i = 1:N traj(i,:) = s'; k1 = deriv(s, GM); k2 = deriv(s+dt/2*k1, GM); k3 = deriv(s+dt/2*k2, GM); k4 = deriv(s+dt*k3, GM); s = s + (dt/6)*(k1 + 2*k2 + 2*k3 + k4); end %% Plot orbit figure(2); plot(traj(:,1), traj(:,2), 'b', 'LineWidth', 1.5); hold on; plot(0,0,'r*', 'MarkerSize', 12); % focus (Sun) axis equal; grid on; title(sprintf('Keplerian Orbit e=%.2f, T={%.4f (analytic), dt=%g', e, T, dt)); xlabel('x (AU)'); ylabel('y (AU)'); %% Verify Kepler's Third Law: T² ∝ a³ a_vals = linspace(0.5, 5, 20); T_analytic = 2*pi .* a_vals.^(3/2) ./ sqrt(GM); figure(3); loglog(a_vals, T_analytic.^2, 'b-o', a_vals, a_vals.^3, 'r--', 'LineWidth', 2); legend({'T² (computed)', 'a³ (reference)'}, 'Location', 'northwest'); title("Kepler's Third Law: T² = 4π²a³/(GM)"); xlabel('Semi-major axis a'); ylabel('T² and a³'); grid on; %% Verify Energy Conservation r = sqrt(traj(:,1).^2 + traj(:,2).^2); v2 = traj(:,3).^2 + traj(:,4).^2; E = 0.5.*v2 - GM./r; % specific energy printf('Energy drift: max|ΔE/E₀| = %.2e\n', max(abs(E-E(1)))/abs(E(1))); %% Verify Angular Momentum Conservation (Kepler II) L = traj(:,1).*traj(:,4) - traj(:,2).*traj(:,3); % r × v printf('Angular momentum drift: max|ΔL/L₀| = %.2e\n', max(abs(L-L(1)))/abs(L(1)));
3 — Pendulum ODE (Large Amplitude)
%% Nonlinear pendulum: θ'' + (g/L)sin(θ) = 0 %% Compare small-angle approximation to exact (RK4) solution %% Plot phase portrait (θ vs θ̇) clear; clc; pkg load odepkg; % or use lsode / manual RK4 g = 9.81; L = 1.0; omega0 = sqrt(g/L); % small-angle angular frequency T0 = 2*pi/omega0; % small-angle period %% ODE: state = [θ; θ̇] function ds = pendulum(t, s, g, L) ds = [s(2); -(g/L)*sin(s(1))]; end figure(4); hold on; grid on; title('Pendulum Phase Portrait (θ vs θ̇)'); xlabel('θ (rad)'); ylabel('dθ/dt (rad/s)'); amplitudes = [5 20 60 120 170] * pi/180; % initial angles colors2 = {'b','g','r','m','c'}; for k = 1:length(amplitudes) th0 = amplitudes(k); s0 = [th0; 0]; tspan = [0, 4*T0]; %% Use lsode for ODE integration f = @(s,t) pendulum(t, s, g, L); t_out = linspace(0, 4*T0, 4000); sol = lsode(f, s0, t_out); plot(sol(:,1), sol(:,2), colors2{k}, 'LineWidth', 1.4, ... 'DisplayName', sprintf('θ₀ = %d°', round(amplitudes(k)*180/pi))); end legend(); axis equal; %% ─── Period vs Amplitude ─────────────────────────────────── th_range = linspace(1, 175, 50) * pi/180; T_exact = zeros(1,length(th_range)); for k = 1:length(th_range) % Numerical integration: T = 4∫₀^θ₀ dθ/sqrt(2(cos θ - cos θ₀)) th0 = th_range(k); th = linspace(0, th0-1e-6, 10000); integrand = 1 ./ sqrt(2*(cos(th) - cos(th0))); T_exact(k) = 4/sqrt(g/L) * trapz(th, integrand); end T_small = T0 * ones(1, length(th_range)); T_corr = T0 * (1 + (1/16)*th_range.^2); % Newton's series correction figure(5); plot(th_range*180/pi, T_exact/T0, 'b', 'LineWidth', 2); hold on; plot(th_range*180/pi, T_small/T0, 'r--', 'LineWidth', 1.5); plot(th_range*180/pi, T_corr/T0, 'g-.', 'LineWidth', 1.5); legend({'Exact (elliptic integral)', 'Small angle T₀', ... "Newton's series T₀(1+θ₀²/16)"}); title('Pendulum Period vs Amplitude'); xlabel('Amplitude θ₀ (degrees)'); ylabel('T/T₀'); grid on;
4 — Projectile with Quadratic Drag
%% Newton Book II — projectile with air resistance F_drag = -c|v|v %% Compare vacuum vs linear drag vs quadratic drag trajectories clear; clc; g = 9.81; m = 0.1; % mass (kg) c = 0.005; % quadratic drag coeff (kg/m) b = 0.01; % linear drag coeff (kg/s) v0 = 50; % initial speed (m/s) alpha = 45*pi/180; %% RK4 for drag ODE: state = [x; y; vx; vy] function ds = drag_ode(state, g, m, c, b, mode) vx = state(3); vy = state(4); speed = sqrt(vx^2+vy^2); switch mode case 0, Fx = 0; Fy = -m*g; % vacuum case 1, Fx = -b*vx; Fy = -m*g - b*vy; % linear case 2, Fx = -c*speed*vx; Fy = -m*g - c*speed*vy;% quadratic end ds = [vx; vy; Fx/m; Fy/m]; end function traj = integrate(s0, g, m, c, b, mode) dt = 0.002; s = s0; traj = s'; while s(2) >= 0 k1 = drag_ode(s, g,m,c,b,mode); k2 = drag_ode(s+dt/2*k1, g,m,c,b,mode); k3 = drag_ode(s+dt/2*k2, g,m,c,b,mode); k4 = drag_ode(s+dt*k3, g,m,c,b,mode); s = s + (dt/6)*(k1+2*k2+2*k3+k4); traj = [traj; s']; end end s0 = [0; 0; v0*cos(alpha); v0*sin(alpha)]; t0 = integrate(s0, g, m, c, b, 0); % vacuum t1 = integrate(s0, g, m, c, b, 1); % linear drag t2 = integrate(s0, g, m, c, b, 2); % quadratic drag figure(6); plot(t0(:,1),t0(:,2),'b', t1(:,1),t1(:,2),'g', t2(:,1),t2(:,2),'r', 'LineWidth',2); legend({'Vacuum (Principia Prop.I-II)', 'Linear drag (Stokes)', ... 'Quadratic drag (Newton Bk II)'}); title('Effect of Air Resistance on Projectile Range'); xlabel('x (m)'); ylabel('y (m)'); grid on; printf('Range vacuum: %.2f m\n', t0(end,1)); printf('Range lin drag: %.2f m\n', t1(end,1)); printf('Range quad drag: %.2f m\n', t2(end,1));
5 — N-Body Gravitational Simulation
%% N-body gravitational simulation using RK4 %% Newton's Law: F_ij = G·mᵢ·mⱼ/|rᵢ-rⱼ|² directed r̂_ij %% Demonstrates Sun-Earth-Moon or figure-8 three-body orbit clear; clc; G = 1.0; % normalised gravitational constant %% ─── Set up Sun-Earth-Jupiter system ───────────────────── m = [1.0; 3e-6; 1e-3]; % masses: Sun, Earth, Jupiter N = length(m); %% Initial positions (AU) and velocities (AU/yr, normalised) r = [0 0; 1.0 0; 5.2 0]; v = [0 0; 0 6.28; 0 2.76]; %% RK4 N-body derivatives function [dr,dv] = nbody_deriv(r, v, m, G) N = length(m); dr = v; dv = zeros(N,2); for i = 1:N for j = 1:N if i ~= j diff = r(j,:) - r(i,:); dist3 = (norm(diff))^3; dv(i,:) += G*m(j)*diff/dist3; end end end end dt = 0.002; Tfinal = 12; % ~12 Earth years steps = round(Tfinal/dt); hist = zeros(steps, N*2); for step = 1:steps hist(step,:) = [r(1,:) r(2,:) r(3,:)]; [dr1,dv1] = nbody_deriv(r, v, m, G); [dr2,dv2] = nbody_deriv(r+dt/2*dr1, v+dt/2*dv1, m, G); [dr3,dv3] = nbody_deriv(r+dt/2*dr2, v+dt/2*dv2, m, G); [dr4,dv4] = nbody_deriv(r+dt*dr3, v+dt*dv3, m, G); r += (dt/6)*(dr1+2*dr2+2*dr3+dr4); v += (dt/6)*(dv1+2*dv2+2*dv3+dv4); end figure(7); plot(hist(:,1),hist(:,2),'y.', hist(:,3),hist(:,4),'b', ... hist(:,5),hist(:,6),'r', 'MarkerSize',3); legend({'Sun','Earth','Jupiter'}); axis equal; grid on; title('Sun-Earth-Jupiter: Newton N-Body (12 yr)'); xlabel('x (AU)'); ylabel('y (AU)');
6 — Solving Kepler's Equation & Visualising the Orbit
%% Kepler's Equation: M = E - e·sin(E) %% Solve for eccentric anomaly E given mean anomaly M %% Then convert to true anomaly θ and plot orbit position vs time clear; clc; function Ea = kepler_solve(M, e, tol) % Newton-Raphson iteration: f(E)=E-e·sin(E)-M, f'(E)=1-e·cos(E) Ea = M; % initial guess for iter = 1:100 f = Ea - e*sin(Ea) - M; fp = 1 - e*cos(Ea); dE = f ./ fp; Ea = Ea - dE; if max(abs(dE)) < tol, break; end end end e = 0.6; % eccentricity (comet-like orbit) a = 1.0; % semi-major axis (AU) GM = 4*pi^2; % in AU³/yr² units T = 2*pi*a^(3/2)/sqrt(GM); M = linspace(0, 2*pi, 1000); Ea = kepler_solve(M, e, 1e-12); % Convert eccentric → true anomaly theta = 2*atan2(sqrt(1+e)*sin(Ea/2), sqrt(1-e)*cos(Ea/2)); r = a*(1-e^2) ./ (1+e*cos(theta)); x = r.*cos(theta); y = r.*sin(theta); figure(8); subplot(1,2,1); plot(x, y, 'b-', 'LineWidth', 2); hold on; plot(0,0,'r*','MarkerSize',14); axis equal; grid on; title(sprintf('Orbit e=%.1f',e)); xlabel('x (AU)'); ylabel('y (AU)'); % Mark equal time intervals (Kepler's 2nd law visually) N_pts = 12; idx = round(linspace(1,1000,N_pts+1)); for k = 1:N_pts i1 = idx(k); i2 = idx(k+1); fill([0 x(i1:i2)' 0],[0 y(i1:i2)' 0],'b','FaceAlpha',0.15,'EdgeColor','none'); end subplot(1,2,2); t_norm = M/2/pi; plot(t_norm, theta*180/pi, 'b', t_norm, t_norm*360, 'r--', 'LineWidth', 2); legend({'True anomaly θ', 'Uniform (circular)'}); title('True anomaly vs time — Kepler II'); xlabel('t/T'); ylabel('θ (degrees)'); grid on; printf('Kepler converged. Max residual: %.2e\n', max(abs(Ea - e*sin(Ea) - M)));
XXII. Interactive Simulations
Simulation I — Orbital Conic Sections
Vary the initial velocity to see how the orbit changes from circle to ellipse to parabola to hyperbola, exactly as Newton's Prop. XI–XIII predict.
Simulation II — Projectile Motion
Newton's Lemma I applied: parabolic trajectory as superposition of uniform horizontal and uniformly accelerated vertical motion (Principia, Book I, Cor. II to Laws).
Simulation III — Pendulum & Phase Space
Compare small-angle (simple harmonic) vs large-angle motion. The phase space portrait shows the separatrix at \(E = 2mgL\) between oscillation and rotation.
Simulation IV — Kepler's Second Law (Equal Areas)
Visualises Proposition I: shaded sectors swept in equal time intervals are equal in area, regardless of where the planet is in its orbit.
Simulation V — Shell Theorem Visualisation
Plots \(g(r)\) inside and outside a uniform sphere, illustrating Props. LXXI–LXXII: \(g\propto r\) inside, \(g\propto 1/r^2\) outside.