11. Vectors and Space Curves
Chapter 10 described curves in the plane with two coordinate functions. This chapter moves to three-dimensional space and brings in the tool that makes 3D calculus manageable: the vector, a quantity with both size and direction. It first builds the geometry of space (vectors, dot and cross products, lines, planes, and quadric surfaces). It then does calculus with vector-valued functions to measure how curves bend and twist, and applies this to motion, ending with a vector derivation of Kepler's first law. These ideas are the working language of physics, computer graphics, robotics, and the multivariable calculus of Chapters 12 to 14.
<a id="111-vectors-and-three-dimensional-space"></a>
11.1 Vectors and Three-Dimensional Space
This subchapter sets up 3D coordinates and the algebra of vectors, then uses the dot and cross products to describe lines, planes, and the curved surfaces given by second-degree equations.
Coordinates in Three Dimensions
Theory
A point on a line needs one number, and a point in the plane needs two. A point in space needs three: a drone's position is not fixed until you know how far east, how far north, and how high it is. We add a third axis, the $z$-axis, perpendicular to both the $x$- and $y$-axes at the origin $O$. A point is then an ordered triple $P = (x, y, z)$, and the set of all triples is $\mathbb{R}^3$.
A good model is the corner of a room. The corner is the origin, the floor edges are the positive $x$- and $y$-axes, and the vertical edge is the positive $z$-axis. The floor is the $xy$-plane ($z = 0$), and the two walls are the $xz$-plane ($y = 0$) and the $yz$-plane ($x = 0$). The three coordinate planes cut space into eight octants. The one with all coordinates positive is the first octant.
The axes have to be arranged in a fixed way, because there are two mirror-image choices. Mathematics, physics, and most graphics libraries use a right-handed system: curl the fingers of your right hand from the positive $x$-axis toward the positive $y$-axis, and your thumb points along the positive $z$-axis. (Some engines, such as Unity and older DirectX code, use left-handed axes. Mixing the two flips the sign of every cross product in Section 11.1.)
z
|
| P = (x, y, z)
| *
| /:
| / : z = height above the xy-plane
| / :
| / :
| / :
|/ :
O------:----------- y
/ \ :
/ \ :
/ \ :
/ \ :
/ * (x, y, 0), the "shadow" of P on the floor
x
Right-handed: curling x toward y, the thumb points along +z.
Distance. The distance between $P_1 = (x_1, y_1, z_1)$ and $P_2 = (x_2, y_2, z_2)$ comes from applying the Pythagorean theorem twice. On the floor, the horizontal distance $d_h$ satisfies $d_h^2 = (x_2 - x_1)^2 + (y_2 - y_1)^2$. The full distance is the hypotenuse of a vertical right triangle with legs $d_h$ and $|z_2 - z_1|$:
$$|P_1P_2| = \sqrt{(x_2 - x_1)^2 + (y_2 - y_1)^2 + (z_2 - z_1)^2}.$$
The midpoint of $P_1P_2$ is $\left(\frac{x_1 + x_2}{2}, \frac{y_1 + y_2}{2}, \frac{z_1 + z_2}{2}\right)$.
Spheres. A sphere is the set of points at a fixed distance $r$ from a center $(h, k, l)$. Squaring the distance formula gives
$$(x - h)^2 + (y - k)^2 + (z - l)^2 = r^2.$$
When a sphere's equation is given in expanded form, such as $x^2 + y^2 + z^2 - 4x + 6y - 3 = 0$, complete the square in each variable to read off the center and radius. If the right side comes out negative, there are no points at all; if it is zero, the "sphere" is a single point.
One equation in 3D is a surface. The biggest adjustment when moving from the plane to space is that the same equation describes a higher-dimensional object:
| Equation | In $\mathbb{R}^2$ | In $\mathbb{R}^3$ |
|---|---|---|
| $x = 3$ | Vertical line | Plane parallel to the $yz$-plane |
| $y = x$ | Line through the origin | Vertical plane containing the $z$-axis |
| $x^2 + y^2 = 1$ | Circle | Infinite circular cylinder around the $z$-axis |
| $x^2 + y^2 + z^2 = 1$ | (not applicable) | Unit sphere |
| $x = 3$ and $y = 2$ | Single point | Vertical line through $(3, 2, 0)$ |
The rule of thumb: each independent equation removes one degree of freedom. In 3D, one equation leaves a 2D surface, two equations leave a 1D curve, and three equations usually leave isolated points. A variable that does not appear is free: $x^2 + y^2 = 1$ says nothing about $z$, so every height is allowed.
Inequalities describe solid regions. For example, $x^2 + y^2 + z^2 \le 4$ is a solid ball, and $1 \le z \le 3$ is a slab between two horizontal planes.
Examples
Distance between $P = (1, -2, 3)$ and $Q = (4, 2, 15)$:
$$|PQ| = \sqrt{3^2 + 4^2 + 12^2} = \sqrt{9 + 16 + 144} = \sqrt{169} = 13.$$
Center and radius of $x^2 + y^2 + z^2 - 4x + 6y - 3 = 0$. Completing the square:
$$(x^2 - 4x + 4) + (y^2 + 6y + 9) + z^2 = 3 + 4 + 9 \quad\Longrightarrow\quad (x - 2)^2 + (y + 3)^2 + z^2 = 16,$$
so the center is $(2, -3, 0)$ and the radius is $4$.
import numpy as np
def distance(P, Q):
return np.sqrt(np.sum((np.asarray(Q, float) - np.asarray(P, float))**2))
P, Q = (1, -2, 3), (4, 2, 15)
print(distance(P, Q)) # 13.0
print((np.array(P) + np.array(Q)) / 2) # [2.5 0. 9. ] midpoint
# Sphere from the expanded form x^2 + y^2 + z^2 + Dx + Ey + Fz + G = 0
def sphere_from_general(D, E, F, G):
center = np.array([-D / 2, -E / 2, -F / 2]) # completing the square in each variable
r2 = center @ center - G # (D/2)^2 + (E/2)^2 + (F/2)^2 - G
if r2 < 0:
return center, None # no real points
return center, np.sqrt(r2)
print(sphere_from_general(-4, 6, 0, -3)) # (array([ 2., -3., 0.]), np.float64(4.0))
print(sphere_from_general(2, 0, 0, 5)) # (array([-1., 0., 0.]), None) -> (x+1)^2 + y^2 + z^2 = -4
# Right-handedness check: the cross product of the x and y unit vectors should be +z
i, j = np.array([1, 0, 0]), np.array([0, 1, 0])
print(np.cross(i, j)) # [0 0 1]
Key Takeaways
- Points in space are triples $(x, y, z)$ in a right-handed system: $x$ curled toward $y$ points the thumb along $z$.
- Distance: $\sqrt{\Delta x^2 + \Delta y^2 + \Delta z^2}$, which is the Pythagorean theorem applied twice.
- A sphere is $(x - h)^2 + (y - k)^2 + (z - l)^2 = r^2$. Complete the square to find the center and radius.
- In $\mathbb{R}^3$, one equation gives a surface and two equations give a curve. A missing variable is free to take any value.
🧪 Practice
- Find the distance from $(3, -1, 2)$ to (a) the origin, (b) the $xy$-plane, (c) the $z$-axis, and (d) the $yz$-plane.
- Find the center and radius of $x^2 + y^2 + z^2 + 8x - 2y + 6z + 1 = 0$. Does the sphere intersect the $xz$-plane? If so, describe the intersection.
- Describe in words and sketch the regions: (a) $y = -2$; (b) $x^2 + z^2 = 9$; (c) $1 \le x^2 + y^2 + z^2 \le 4$, $z \ge 0$; (d) the set of points equidistant from $(0, 0, 0)$ and $(2, 0, 0)$.
- Interview-style: A collision checker compares the distance between the centers of thousands of pairs of spheres against the sum of their radii every frame. How would you speed up the comparison without changing its result? Hint: the square root is monotonic, so compare $\Delta x^2 + \Delta y^2 + \Delta z^2$ with $(r_1 + r_2)^2$ instead.
Vector Operations
Theory
Many quantities are not described by a single number. "The wind blows at 40 km/h" is incomplete until you add "from the northwest". A displacement, a velocity, a force, and an electric field all have a magnitude (how much) and a direction (which way). A vector packages both.
Geometrically, a vector is an arrow. Two arrows with the same length and direction represent the same vector, wherever they are drawn. A vector records a move, not a location: "3 steps east and 4 steps north" is the same instruction from any starting point. The arrow from $A$ to $B$ is written $\overrightarrow{AB}$.
Algebraically, place the tail at the origin. The head lands at a point $(a_1, a_2, a_3)$, and we write the vector by its components:
$$\mathbf{a} = \langle a_1, a_2, a_3 \rangle.$$
The angle brackets distinguish the vector (a move) from the point $(a_1, a_2, a_3)$ (a location). The vector from the origin to a point $P$ is the position vector of $P$. For $A = (x_1, y_1, z_1)$ and $B = (x_2, y_2, z_2)$,
$$\overrightarrow{AB} = \langle x_2 - x_1,\ y_2 - y_1,\ z_2 - z_1 \rangle \qquad (\text{head minus tail}).$$
Operations. All are done component by component.
| Operation | Formula | Geometric meaning |
|---|---|---|
| Addition | $\mathbf{a} + \mathbf{b} = \langle a_1 + b_1, a_2 + b_2, a_3 + b_3 \rangle$ | Place the tail of $\mathbf{b}$ at the head of $\mathbf{a}$ (triangle law); the diagonal of the parallelogram |
| Scalar multiple | $c\mathbf{a} = \langle ca_1, ca_2, ca_3 \rangle$ | Stretch by $\lvert c \rvert$; reverse direction if $c < 0$ |
| Subtraction | $\mathbf{a} - \mathbf{b} = \mathbf{a} + (-1)\mathbf{b}$ | The arrow from the head of $\mathbf{b}$ to the head of $\mathbf{a}$ |
| Magnitude | $\lvert \mathbf{a} \rvert = \sqrt{a_1^2 + a_2^2 + a_3^2}$ | Length of the arrow |
| Unit vector | $\mathbf{u} = \dfrac{\mathbf{a}}{\lvert \mathbf{a} \rvert}$ for $\mathbf{a} \ne \mathbf{0}$ | Direction of $\mathbf{a}$ with length $1$ |
Triangle law: a + b Subtraction: a - b
* head of a + b a - b
.'| *----------------->*
.' | \ /
a + b .' | b \ b / a
.' | \ /
.' | \ /
O---------->* \ /
a \ /
\ /
Parallelogram: drawing b first, then a, O
reaches the same point, so a + b = b + a. a - b runs from the head of b to the head of a
Two vectors are parallel if one is a scalar multiple of the other. The zero vector $\mathbf{0} = \langle 0, 0, 0 \rangle$ has length $0$ and no direction, and it counts as parallel to every vector.
Standard basis vectors. The unit vectors along the axes are
$$\mathbf{i} = \langle 1, 0, 0 \rangle, \quad \mathbf{j} = \langle 0, 1, 0 \rangle, \quad \mathbf{k} = \langle 0, 0, 1 \rangle,$$
and every vector splits uniquely as $\langle a_1, a_2, a_3 \rangle = a_1\mathbf{i} + a_2\mathbf{j} + a_3\mathbf{k}$. Both notations are common, and they mean the same thing.
Algebraic rules. For vectors $\mathbf{a}, \mathbf{b}, \mathbf{c}$ and scalars $c, d$:
- $\mathbf{a} + \mathbf{b} = \mathbf{b} + \mathbf{a}$ and $(\mathbf{a} + \mathbf{b}) + \mathbf{c} = \mathbf{a} + (\mathbf{b} + \mathbf{c})$
- $\mathbf{a} + \mathbf{0} = \mathbf{a}$ and $\mathbf{a} + (-\mathbf{a}) = \mathbf{0}$
- $c(\mathbf{a} + \mathbf{b}) = c\mathbf{a} + c\mathbf{b}$, $(c + d)\mathbf{a} = c\mathbf{a} + d\mathbf{a}$, and $(cd)\mathbf{a} = c(d\mathbf{a})$
- $\lvert c\mathbf{a} \rvert = \lvert c \rvert,\lvert \mathbf{a} \rvert$ and the triangle inequality $\lvert \mathbf{a} + \mathbf{b} \rvert \le \lvert \mathbf{a} \rvert + \lvert \mathbf{b} \rvert$
These are exactly the rules of ordinary arithmetic, so vector equations can be manipulated like algebraic ones, as long as you never divide by a vector.
Physics: forces add as vectors. When several forces act on an object, their combined effect is the single resultant force, their vector sum. An object is in equilibrium when the resultant is $\mathbf{0}$.
Examples
An airplane flies with airspeed $500$ km/h on a heading of N30°E (30° east of north). A wind blows due east at $60$ km/h. With $x$ pointing east and $y$ north, the air velocity is $500\langle \sin 30°, \cos 30° \rangle = \langle 250, 433.01 \rangle$, and the ground velocity is
$$\mathbf{v}_{\text{ground}} = \langle 250, 433.01 \rangle + \langle 60, 0 \rangle = \langle 310, 433.01 \rangle,$$
so the ground speed is $\sqrt{310^2 + 433.01^2} \approx 532.5$ km/h on a bearing of about N35.6°E.
import numpy as np
# Airplane: airspeed 500 km/h on heading N30E; wind 60 km/h toward the east
heading = np.radians(30) # measured from north toward east
v_air = 500 * np.array([np.sin(heading), np.cos(heading), 0]) # x = east, y = north, z = up
v_wind = np.array([60, 0, 0])
v_ground = v_air + v_wind # velocities add as vectors
speed = np.linalg.norm(v_ground)
bearing = np.degrees(np.arctan2(v_ground[0], v_ground[1])) # atan2(east, north) = angle from north
print(np.round(v_ground, 2), round(speed, 2), round(bearing, 2))
# [310. 433.01 0. ] 532.54 35.6
# Magnitude and unit vector
v = np.array([3.0, 4.0, 12.0])
print(np.linalg.norm(v), v / np.linalg.norm(v))
# 13.0 [0.23076923 0.30769231 0.92307692]
# Resultant of three forces (newtons), one per row
F = np.array([[10, 0, 0], [0, -5, 5], [-4, 2, -1]])
print(F.sum(axis=0), np.linalg.norm(F.sum(axis=0)))
# [ 6 -3 4] 7.810249675906654 -> a fourth force of <-6, 3, -4> would give equilibrium
Key Takeaways
- A vector has magnitude and direction but no fixed position. $\overrightarrow{AB}$ is "head minus tail".
- Addition, subtraction, and scalar multiplication work componentwise and obey the usual algebraic rules.
- $\lvert \mathbf{a} \rvert = \sqrt{a_1^2 + a_2^2 + a_3^2}$, and $\mathbf{a}/\lvert \mathbf{a} \rvert$ is the unit vector in the direction of $\mathbf{a}$.
- Every vector is $a_1\mathbf{i} + a_2\mathbf{j} + a_3\mathbf{k}$. Forces and velocities combine by vector addition.
🧪 Practice
- For $\mathbf{a} = \langle 2, -1, 3 \rangle$ and $\mathbf{b} = \langle -4, 0, 1 \rangle$, compute $\mathbf{a} + \mathbf{b}$, $3\mathbf{a} - 2\mathbf{b}$, $\lvert \mathbf{a} - \mathbf{b} \rvert$, and a vector of length $5$ in the direction of $\mathbf{b}$.
- Show that $A = (1, 2, 3)$, $B = (3, 5, 4)$, and $C = (7, 11, 6)$ are collinear by comparing $\overrightarrow{AB}$ and $\overrightarrow{AC}$.
- A 100 N weight hangs from two ropes attached to the ceiling at angles of 30° and 45° from the horizontal, on opposite sides. Write each rope's tension as a vector and use equilibrium to find both tensions.
- Interview-style: In a physics engine, positions and directions are both stored as 3-component arrays. Why do graphics libraries often store them as 4-component arrays $(x, y, z, 1)$ for points and $(x, y, z, 0)$ for directions? Hint: a translation should move a point but not a direction. Look at what a $4 \times 4$ matrix whose last column holds a translation does to each kind.
Dot Product and Projections
Theory
Adding vectors combines them. Often, though, we want to compare two vectors: how much of a force acts along the direction of motion, or how closely two directions agree. The dot product answers this with a single number.
Definition. For $\mathbf{a} = \langle a_1, a_2, a_3 \rangle$ and $\mathbf{b} = \langle b_1, b_2, b_3 \rangle$,
$$\mathbf{a} \cdot \mathbf{b} = a_1b_1 + a_2b_2 + a_3b_3.$$
The result is a scalar, so it is also called the scalar product. Its algebraic properties are the familiar ones:
| Property | Statement |
|---|---|
| Length | $\mathbf{a} \cdot \mathbf{a} = \lvert \mathbf{a} \rvert^2$ |
| Commutative | $\mathbf{a} \cdot \mathbf{b} = \mathbf{b} \cdot \mathbf{a}$ |
| Distributive | $\mathbf{a} \cdot (\mathbf{b} + \mathbf{c}) = \mathbf{a} \cdot \mathbf{b} + \mathbf{a} \cdot \mathbf{c}$ |
| Scalars | $(c\mathbf{a}) \cdot \mathbf{b} = c(\mathbf{a} \cdot \mathbf{b}) = \mathbf{a} \cdot (c\mathbf{b})$ |
| Zero | $\mathbf{0} \cdot \mathbf{a} = 0$ |
Geometric meaning. If $\theta$ is the angle between $\mathbf{a}$ and $\mathbf{b}$ (with $0 \le \theta \le \pi$), then
$$\mathbf{a} \cdot \mathbf{b} = \lvert \mathbf{a} \rvert,\lvert \mathbf{b} \rvert \cos\theta.$$
Why: the vectors $\mathbf{a}$, $\mathbf{b}$, and $\mathbf{a} - \mathbf{b}$ form a triangle. The law of cosines gives $\lvert \mathbf{a} - \mathbf{b} \rvert^2 = \lvert \mathbf{a} \rvert^2 + \lvert \mathbf{b} \rvert^2 - 2\lvert \mathbf{a} \rvert\lvert \mathbf{b} \rvert\cos\theta$. Expanding the left side with the properties above gives $\lvert \mathbf{a} \rvert^2 - 2,\mathbf{a} \cdot \mathbf{b} + \lvert \mathbf{b} \rvert^2$. Comparing the two expressions proves the formula.
The sign of the dot product tells you the kind of angle:
| $\mathbf{a} \cdot \mathbf{b}$ | Angle | Interpretation |
|---|---|---|
| $> 0$ | Acute | The vectors point "roughly the same way" |
| $= 0$ | Right angle | Orthogonal (perpendicular) |
| $< 0$ | Obtuse | The vectors point "roughly opposite ways" |
The orthogonality test $\mathbf{a} \cdot \mathbf{b} = 0$ is the most used fact in this chapter. Since $\lvert \cos\theta \rvert \le 1$, the formula also gives the Cauchy-Schwarz inequality $\lvert \mathbf{a} \cdot \mathbf{b} \rvert \le \lvert \mathbf{a} \rvert,\lvert \mathbf{b} \rvert$.
Direction angles. The angles $\alpha, \beta, \gamma$ that a nonzero vector $\mathbf{a}$ makes with the positive axes satisfy $\cos\alpha = \frac{a_1}{\lvert \mathbf{a} \rvert}$, $\cos\beta = \frac{a_2}{\lvert \mathbf{a} \rvert}$, $\cos\gamma = \frac{a_3}{\lvert \mathbf{a} \rvert}$. These direction cosines are the components of the unit vector $\frac{\mathbf{a}}{\lvert \mathbf{a} \rvert}$, so $\cos^2\alpha + \cos^2\beta + \cos^2\gamma = 1$.
Projections. Shine a light perpendicular to $\mathbf{a}$ and look at the shadow that $\mathbf{b}$ casts on the line through $\mathbf{a}$. The signed length of that shadow is the scalar projection (or component) of $\mathbf{b}$ onto $\mathbf{a}$, and the shadow itself is the vector projection:
$$\operatorname{comp}{\mathbf{a}}\mathbf{b} = \frac{\mathbf{a} \cdot \mathbf{b}}{\lvert \mathbf{a} \rvert} = \lvert \mathbf{b} \rvert\cos\theta, \qquad \operatorname{proj}{\mathbf{a}}\mathbf{b} = \frac{\mathbf{a} \cdot \mathbf{b}}{\lvert \mathbf{a} \rvert^2},\mathbf{a} = \frac{\mathbf{a} \cdot \mathbf{b}}{\mathbf{a} \cdot \mathbf{a}},\mathbf{a}.$$
* head of b
/|
/ |
b / | b - proj_a(b) (orthogonal to a)
/ |
/ |
/ th |
O======*----------------> a
|------|
proj_a(b): the "shadow" of b on the line of a
The leftover piece $\mathbf{b} - \operatorname{proj}_{\mathbf{a}}\mathbf{b}$ is orthogonal to $\mathbf{a}$. So every vector splits as (part along $\mathbf{a}$) + (part perpendicular to $\mathbf{a}$). This orthogonal decomposition is used constantly: splitting gravity into components along and across a ramp, splitting acceleration into tangential and normal parts (Section 11.4), and least-squares fitting in linear algebra.
Work. A constant force $\mathbf{F}$ moving an object along a displacement $\mathbf{D}$ does work $W = \mathbf{F} \cdot \mathbf{D} = \lvert \mathbf{F} \rvert,\lvert \mathbf{D} \rvert\cos\theta$. Only the component of the force along the motion does work. This generalizes the $W = Fd$ of Section 7.3, and Section 14.2 extends it to curved paths.
Examples
The angle between a main diagonal of a cube and one of its edges. Take the unit cube with $\mathbf{a} = \langle 1, 1, 1 \rangle$ and $\mathbf{b} = \langle 1, 0, 0 \rangle$. Then $\cos\theta = \frac{1}{\sqrt{3}\cdot 1}$, so $\theta = \arccos\frac{1}{\sqrt{3}} \approx 54.74°$.
Projection of $\mathbf{b} = \langle 3, -3, 6 \rangle$ onto $\mathbf{a} = \langle 2, 1, -2 \rangle$. Here $\mathbf{a} \cdot \mathbf{b} = 6 - 3 - 12 = -9$ and $\lvert \mathbf{a} \rvert = 3$, so
$$\operatorname{comp}{\mathbf{a}}\mathbf{b} = -3, \qquad \operatorname{proj}{\mathbf{a}}\mathbf{b} = \frac{-9}{9}\langle 2, 1, -2 \rangle = \langle -2, -1, 2 \rangle.$$
The negative component means $\mathbf{b}$ leans away from $\mathbf{a}$. The perpendicular part is $\langle 3, -3, 6 \rangle - \langle -2, -1, 2 \rangle = \langle 5, -2, 4 \rangle$, and indeed $\langle 5, -2, 4 \rangle \cdot \langle 2, 1, -2 \rangle = 10 - 2 - 8 = 0$.
import numpy as np
def angle_deg(a, b):
c = a @ b / (np.linalg.norm(a) * np.linalg.norm(b)) # cos(theta)
return np.degrees(np.arccos(np.clip(c, -1.0, 1.0))) # clip guards against rounding just past +-1
def proj(a, b):
"""Vector projection of b onto a."""
return (a @ b) / (a @ a) * a # a @ a avoids a square root
diag, edge = np.array([1, 1, 1]), np.array([1, 0, 0])
print(angle_deg(diag, edge)) # 54.735610317245346
a, b = np.array([2, 1, -2]), np.array([3, -3, 6])
p = proj(a, b)
perp = b - p
print(a @ b / np.linalg.norm(a), p, perp, perp @ a)
# -3.0 [-2. -1. 2.] [ 5. -2. 4.] 0.0 component, projection, perpendicular part, check
# Work: pulling a sled 40 m with a 70 N force at 35 degrees above the horizontal
F = 70 * np.array([np.cos(np.radians(35)), np.sin(np.radians(35))])
D = np.array([40, 0])
print(round(F @ D, 1)) # 2293.6 joules (the vertical part does no work)
# Why clip: the angle between a vector and itself should be exactly 0
u = np.array([1.0, 1.0, 1.0])
c = u @ u / (np.linalg.norm(u) * np.linalg.norm(u))
print(repr(c), np.arccos(c)) # np.float64(1.0000000000000002) nan (plus a RuntimeWarning)
Key Takeaways
- $\mathbf{a} \cdot \mathbf{b} = a_1b_1 + a_2b_2 + a_3b_3 = \lvert \mathbf{a} \rvert,\lvert \mathbf{b} \rvert\cos\theta$, a scalar.
- $\mathbf{a} \cdot \mathbf{b} = 0$ exactly when the vectors are orthogonal. The sign tells you whether the angle is acute or obtuse.
- $\operatorname{comp}{\mathbf{a}}\mathbf{b} = \frac{\mathbf{a} \cdot \mathbf{b}}{\lvert \mathbf{a} \rvert}$ and $\operatorname{proj}{\mathbf{a}}\mathbf{b} = \frac{\mathbf{a} \cdot \mathbf{b}}{\mathbf{a} \cdot \mathbf{a}}\mathbf{a}$. The remainder $\mathbf{b} - \operatorname{proj}_{\mathbf{a}}\mathbf{b}$ is orthogonal to $\mathbf{a}$.
- Work done by a constant force is $W = \mathbf{F} \cdot \mathbf{D}$.
- In code, clip the cosine to $[-1, 1]$ before calling
arccos.
🧪 Practice
- Find the angle between $\langle 1, 2, -2 \rangle$ and $\langle 4, 0, -3 \rangle$, and find the direction angles of $\langle 1, 2, -2 \rangle$.
- Find $\operatorname{comp}{\mathbf{a}}\mathbf{b}$ and $\operatorname{proj}{\mathbf{a}}\mathbf{b}$ for $\mathbf{a} = \langle 1, 1, 0 \rangle$, $\mathbf{b} = \langle 2, -1, 4 \rangle$. Write $\mathbf{b}$ as the sum of a vector parallel to $\mathbf{a}$ and a vector orthogonal to $\mathbf{a}$.
- A box slides 8 m down a ramp inclined at 20°. Find the work done by gravity on a 15 kg box ($g = 9.8$ m/s²), first with the dot product and then by noting that only the vertical drop matters. Explain why the answers agree.
- Interview-style: A recommendation system ranks items by the cosine similarity $\frac{\mathbf{a} \cdot \mathbf{b}}{\lvert \mathbf{a} \rvert\lvert \mathbf{b} \rvert}$ of embedding vectors rather than by the raw dot product. When would the two rankings differ, and how would you make the cosine search as cheap as a dot-product search? Hint: the raw dot product rewards long vectors (often popular items). Normalize every vector to unit length once, when it is stored.
Cross Product and Applications
Theory
The dot product turns two vectors into a number. In 3D, there is a second product that turns two vectors into a new vector perpendicular to both. This is exactly what is needed to find the normal direction of a plane, the axis of a rotation, or the torque of a force.
Definition. For $\mathbf{a} = \langle a_1, a_2, a_3 \rangle$ and $\mathbf{b} = \langle b_1, b_2, b_3 \rangle$,
$$\mathbf{a} \times \mathbf{b} = \begin{vmatrix} \mathbf{i} & \mathbf{j} & \mathbf{k} \ a_1 & a_2 & a_3 \ b_1 & b_2 & b_3 \end{vmatrix} = \langle a_2b_3 - a_3b_2,\ a_3b_1 - a_1b_3,\ a_1b_2 - a_2b_1 \rangle.$$
The determinant is a memory aid, expanded along the first row (remember the minus sign on the $\mathbf{j}$ term). The cross product is defined only in three dimensions.
Geometric meaning.
- Direction: $\mathbf{a} \times \mathbf{b}$ is orthogonal to both $\mathbf{a}$ and $\mathbf{b}$ (check: $\mathbf{a} \cdot (\mathbf{a} \times \mathbf{b}) = 0$ by expansion). Of the two perpendicular directions, it points the way your right thumb points when your fingers curl from $\mathbf{a}$ toward $\mathbf{b}$.
- Length: $\lvert \mathbf{a} \times \mathbf{b} \rvert = \lvert \mathbf{a} \rvert,\lvert \mathbf{b} \rvert\sin\theta$, which is the area of the parallelogram spanned by $\mathbf{a}$ and $\mathbf{b}$. This follows from Lagrange's identity $\lvert \mathbf{a} \times \mathbf{b} \rvert^2 = \lvert \mathbf{a} \rvert^2\lvert \mathbf{b} \rvert^2 - (\mathbf{a} \cdot \mathbf{b})^2$, which can be checked by expanding both sides.
Consequently, $\mathbf{a} \times \mathbf{b} = \mathbf{0}$ exactly when $\mathbf{a}$ and $\mathbf{b}$ are parallel (the parallelogram is flat). The dot product detects perpendicularity, and the cross product detects parallelism.
Basis products: go around the cycle Parallelogram area
i a x b (perpendicular to a and b)
^ \ ^
/ \ forward: i x j = k | *--------------*
/ v j x k = i | / /
k <---- j k x i = j | / b /
| / area = /
backward gives a minus sign: |/ |a||b|sin th
j x i = -k, k x j = -i, i x k = -j O--------------*
a
Properties. Some familiar rules fail, so this table deserves care:
| Property | Holds? |
|---|---|
| $\mathbf{a} \times \mathbf{b} = -(\mathbf{b} \times \mathbf{a})$ | Yes (anticommutative: order matters) |
| $\mathbf{a} \times \mathbf{a} = \mathbf{0}$ | Yes |
| $\mathbf{a} \times (\mathbf{b} + \mathbf{c}) = \mathbf{a} \times \mathbf{b} + \mathbf{a} \times \mathbf{c}$ | Yes (distributive) |
| $(c\mathbf{a}) \times \mathbf{b} = c(\mathbf{a} \times \mathbf{b})$ | Yes |
| $(\mathbf{a} \times \mathbf{b}) \times \mathbf{c} = \mathbf{a} \times (\mathbf{b} \times \mathbf{c})$ | No (not associative) |
| $\mathbf{a} \times (\mathbf{b} \times \mathbf{c}) = (\mathbf{a} \cdot \mathbf{c})\mathbf{b} - (\mathbf{a} \cdot \mathbf{b})\mathbf{c}$ | Yes (the "BAC minus CAB" rule) |
Scalar triple product. The number $\mathbf{a} \cdot (\mathbf{b} \times \mathbf{c})$ equals the determinant of the $3 \times 3$ matrix with rows $\mathbf{a}, \mathbf{b}, \mathbf{c}$. Its absolute value is the volume of the parallelepiped (slanted box) spanned by the three vectors: the base area is $\lvert \mathbf{b} \times \mathbf{c} \rvert$, and the height is the component of $\mathbf{a}$ along the base's normal. In particular, three vectors are coplanar exactly when their triple product is $0$. Its sign tells you whether $(\mathbf{a}, \mathbf{b}, \mathbf{c})$ is right-handed (positive) or left-handed (negative).
Torque. When a force $\mathbf{F}$ acts at a point with position $\mathbf{r}$ relative to a pivot, the turning effect is the torque $\boldsymbol{\tau} = \mathbf{r} \times \mathbf{F}$. Its magnitude $\lvert \mathbf{r} \rvert,\lvert \mathbf{F} \rvert\sin\theta$ explains why you push a door far from the hinge and at right angles to it. Its direction is the axis of rotation. Angular momentum $\mathbf{r} \times m\mathbf{v}$, used in Section 11.4, has the same form.
Examples
Area of the triangle with vertices $P = (1, 0, 2)$, $Q = (3, 1, 0)$, $R = (0, 4, 1)$. Its two edge vectors are $\overrightarrow{PQ} = \langle 2, 1, -2 \rangle$ and $\overrightarrow{PR} = \langle -1, 4, -1 \rangle$, so
$$\overrightarrow{PQ} \times \overrightarrow{PR} = \langle (1)(-1) - (-2)(4),\ (-2)(-1) - (2)(-1),\ (2)(4) - (1)(-1) \rangle = \langle 7, 4, 9 \rangle.$$
The triangle is half the parallelogram, so its area is $\frac{1}{2}\sqrt{49 + 16 + 81} = \frac{\sqrt{146}}{2} \approx 6.04$. The vector $\langle 7, 4, 9 \rangle$ is also a normal to the plane through $P$, $Q$, $R$, which is used in the next topic.
import numpy as np
P, Q, R = np.array([1, 0, 2]), np.array([3, 1, 0]), np.array([0, 4, 1])
n = np.cross(Q - P, R - P) # normal to the triangle's plane
print(n, n @ (Q - P), n @ (R - P)) # [7 4 9] 0 0 orthogonal to both edges
print(0.5 * np.linalg.norm(n)) # 6.041522986797286 area of triangle PQR
# Scalar triple product = determinant = signed volume of the parallelepiped
a, b, c = np.array([1, 2, 0]), np.array([0, 1, 3]), np.array([2, 0, 1])
print(a @ np.cross(b, c), np.linalg.det(np.array([a, b, c]))) # 13 13.0
# Coplanarity test for four points: zero volume
def coplanar(A, B, C, D, tol=1e-9):
vol = (B - A) @ np.cross(C - A, D - A)
return abs(vol) <= tol # tolerance: floating-point inputs are never exact
print(coplanar(P, Q, R, P + 2 * (Q - P) - (R - P))) # True: D is built inside the plane of P, Q, R
print(coplanar(P, Q, R, np.array([0, 0, 0]))) # False
# Torque of a 50 N force on a 0.3 m wrench, applied at 60 degrees to the handle
r = np.array([0.3, 0, 0])
F = 50 * np.array([np.cos(np.radians(60)), np.sin(np.radians(60)), 0])
print(np.round(np.cross(r, F), 3)) # [ 0. 0. 12.99] N*m, about the z-axis
# Not associative: (i x i) x j != i x (i x j)
i, j = np.eye(3)[0], np.eye(3)[1]
print(np.cross(np.cross(i, i), j), np.cross(i, np.cross(i, j))) # [0. 0. 0.] [ 0. -1. 0.]
Key Takeaways
- $\mathbf{a} \times \mathbf{b}$ is a vector orthogonal to both factors, oriented by the right-hand rule, with length $\lvert \mathbf{a} \rvert\lvert \mathbf{b} \rvert\sin\theta$.
- $\lvert \mathbf{a} \times \mathbf{b} \rvert$ is the area of the parallelogram they span; half of it is the triangle's area.
- $\mathbf{a} \times \mathbf{b} = \mathbf{0}$ exactly when the vectors are parallel. The product is anticommutative and not associative.
- $\lvert \mathbf{a} \cdot (\mathbf{b} \times \mathbf{c}) \rvert$ is the parallelepiped's volume; zero means coplanar.
- Torque is $\boldsymbol{\tau} = \mathbf{r} \times \mathbf{F}$.
🧪 Practice
- Compute $\langle 1, -2, 3 \rangle \times \langle 4, 0, -1 \rangle$ and verify that the result is orthogonal to both factors. Then compute the product in the other order.
- Find a unit vector orthogonal to the plane through $(0, 0, 1)$, $(2, 1, 0)$, $(1, 3, 2)$, and find the area of the triangle with those vertices.
- Use the scalar triple product to decide whether $(1, 0, 1)$, $(2, 1, 3)$, $(0, 2, 1)$, and $(3, 3, 5)$ are coplanar. If not, find the volume of the tetrahedron they form (one sixth of the parallelepiped).
- Interview-style: In a triangle mesh, how would you decide whether a triangle's vertices $A, B, C$ are listed counterclockwise as seen from a camera at point $E$, and why does this matter for rendering? Hint: look at the sign of $\big((B - A) \times (C - A)\big) \cdot (E - A)$. Renderers use this sign for back-face culling, skipping triangles that face away from the camera.
Lines and Planes in Space
Theory
In the plane, a line is pinned down by a point and a slope. In space, a single slope cannot describe a direction, because directions have two degrees of freedom. Vectors fix this: a point and a direction vector determine a line, and a point and a normal vector determine a plane.
Lines. The line through the point $P_0$ (position vector $\mathbf{r}_0$) parallel to the vector $\mathbf{v} = \langle a, b, c \rangle$ consists of the points you reach by starting at $P_0$ and moving some multiple of $\mathbf{v}$:
$$\mathbf{r}(t) = \mathbf{r}_0 + t\mathbf{v}, \qquad t \in \mathbb{R}.$$
In components, this is the set of parametric equations $x = x_0 + at$, $y = y_0 + bt$, $z = z_0 + ct$. When $a, b, c$ are all nonzero, solving each one for $t$ gives the symmetric equations
$$\frac{x - x_0}{a} = \frac{y - y_0}{b} = \frac{z - z_0}{c}.$$
Restricting $t$ to $[0, 1]$ with $\mathbf{v} = \mathbf{r}_1 - \mathbf{r}_0$ gives the segment from $\mathbf{r}_0$ to $\mathbf{r}_1$: $\mathbf{r}(t) = (1 - t)\mathbf{r}_0 + t\mathbf{r}_1$. This is the linear interpolation ("lerp") used throughout graphics.
Two lines in space can relate in three ways, not two:
| Relationship | Test |
|---|---|
| Parallel (or identical) | Direction vectors are multiples: $\mathbf{v}_1 \times \mathbf{v}_2 = \mathbf{0}$ |
| Intersecting | Not parallel, and the equations $\mathbf{r}_1(t) = \mathbf{r}_2(s)$ have a solution |
| Skew | Not parallel and never meet. This is impossible in the plane: think of a road and a bridge crossing over it |
Planes. A plane is not determined by a single direction inside it, but it is determined by the one direction perpendicular to it. Let $\mathbf{n} = \langle a, b, c \rangle$ be a normal vector and $P_0 = (x_0, y_0, z_0)$ a point of the plane. A point $P$ lies in the plane exactly when the vector $\overrightarrow{P_0P}$ lies flat in it, that is, when it is orthogonal to $\mathbf{n}$:
$$\mathbf{n} \cdot (\mathbf{r} - \mathbf{r}_0) = 0 \quad\Longleftrightarrow\quad a(x - x_0) + b(y - y_0) + c(z - z_0) = 0 \quad\Longleftrightarrow\quad ax + by + cz = d,$$
with $d = ax_0 + by_0 + cz_0$. The coefficients of $x, y, z$ in a linear equation are a normal vector. This lets you read geometry straight from an equation.
n (normal)
^
|
+-------|------------------------+
/ | /
/ * P0 / every vector P0 -> P lying in
/ \ / the plane is orthogonal to n:
/ \ / n . (r - r0) = 0
/ * P /
+--------------------------------+
Standard constructions:
- Plane through three points $P, Q, R$: use $\mathbf{n} = \overrightarrow{PQ} \times \overrightarrow{PR}$ and the point $P$.
- Angle between two planes: the angle between their normals, $\cos\theta = \frac{\lvert \mathbf{n}_1 \cdot \mathbf{n}_2 \rvert}{\lvert \mathbf{n}_1 \rvert\lvert \mathbf{n}_2 \rvert}$. Planes are parallel when their normals are parallel.
- Line of intersection of two nonparallel planes: it lies in both, so it is orthogonal to both normals, and its direction is $\mathbf{v} = \mathbf{n}_1 \times \mathbf{n}_2$. Find one common point by setting a variable to $0$ and solving.
- Line meets plane: substitute $x(t), y(t), z(t)$ into $ax + by + cz = d$ and solve for $t$. If the coefficient of $t$ is $\mathbf{n} \cdot \mathbf{v} = 0$, the line is parallel to the plane (no solution, or the line lies in it).
Distances. Both formulas are projections onto a normal direction:
$$\text{point } (x_1, y_1, z_1) \text{ to plane } ax + by + cz = d: \quad D = \frac{\lvert ax_1 + by_1 + cz_1 - d \rvert}{\sqrt{a^2 + b^2 + c^2}},$$
$$\text{skew lines } \mathbf{r}_1 + t\mathbf{v}_1 \text{ and } \mathbf{r}_2 + s\mathbf{v}_2: \quad D = \frac{\lvert (\mathbf{r}_2 - \mathbf{r}_1) \cdot (\mathbf{v}_1 \times \mathbf{v}_2) \rvert}{\lvert \mathbf{v}_1 \times \mathbf{v}_2 \rvert}.$$
For skew lines, $\mathbf{v}_1 \times \mathbf{v}_2$ is perpendicular to both lines, and the shortest connecting segment points along it. Projecting any connecting vector onto that direction gives the gap.
Examples
The plane through $(1, 0, 0)$, $(0, 2, 0)$, $(0, 0, 3)$. The edge vectors $\langle -1, 2, 0 \rangle$ and $\langle -1, 0, 3 \rangle$ have cross product $\langle 6, 3, 2 \rangle$, so the plane is $6x + 3y + 2z = 6$, or $\frac{x}{1} + \frac{y}{2} + \frac{z}{3} = 1$ (the intercept form). The line $\mathbf{r}(t) = t\langle 1, 1, 1 \rangle$ meets it when $6t + 3t + 2t = 6$, at $t = \frac{6}{11}$. The distance from $(2, 2, 2)$ to the plane is $\frac{\lvert 12 + 6 + 4 - 6 \rvert}{7} = \frac{16}{7}$.
The planes $x + y + z = 1$ and $x - 2y + 3z = 1$ meet in a line with direction $\langle 1, 1, 1 \rangle \times \langle 1, -2, 3 \rangle = \langle 5, -2, -3 \rangle$. Setting $z = 0$ gives $x + y = 1$, $x - 2y = 1$, so $(1, 0, 0)$ is on both planes, and the line is $\mathbf{r}(t) = \langle 1 + 5t, -2t, -3t \rangle$. The angle between the planes is $\arccos\frac{2}{\sqrt{3}\sqrt{14}} \approx 72.0°$.
import numpy as np
def plane_through(P, Q, R):
n = np.cross(Q - P, R - P)
return n, n @ P # plane: n . x = d
def intersect_line_plane(P0, v, n, d, eps=1e-12):
denom = n @ v # zero when the line is parallel to the plane
if abs(denom) < eps:
return None
t = (d - n @ P0) / denom # solve n . (P0 + t v) = d for t
return P0 + t * v
def point_plane_distance(X, n, d):
return abs(n @ X - d) / np.linalg.norm(n)
def line_line_distance(P1, v1, P2, v2):
w = np.cross(v1, v2) # common perpendicular direction
if np.allclose(w, 0): # parallel lines: distance from P2 to line 1
return np.linalg.norm(np.cross(P2 - P1, v1)) / np.linalg.norm(v1)
return abs((P2 - P1) @ w) / np.linalg.norm(w)
n, d = plane_through(np.array([1., 0, 0]), np.array([0., 2, 0]), np.array([0., 0, 3]))
print(n, d) # [6. 3. 2.] 6.0 6x + 3y + 2z = 6
print(intersect_line_plane(np.array([0., 0, 0]), np.array([1., 1, 1]), n, d))
# [0.54545455 0.54545455 0.54545455] t = 6/11
print(point_plane_distance(np.array([2., 2, 2]), n, d)) # 2.2857142857142856 = 16/7
# Skew lines: L1 through (1,0,0) along <1,1,0>; L2 through (0,0,2) along <0,1,1>
print(line_line_distance(np.array([1., 0, 0]), np.array([1., 1, 0]),
np.array([0., 0, 2]), np.array([0., 1, 1])))
# 0.5773502691896258 = 1/sqrt(3)
Key Takeaways
- A line is $\mathbf{r}_0 + t\mathbf{v}$: a point plus multiples of a direction. Lines in space can be parallel, intersecting, or skew.
- A plane is $\mathbf{n} \cdot (\mathbf{r} - \mathbf{r}_0) = 0$, or $ax + by + cz = d$ with normal $\langle a, b, c \rangle$.
- Use cross products to build normals (plane through three points) and directions (intersection of two planes).
- Point-to-plane and skew-line distances are projections onto a normal direction.
🧪 Practice
- Find parametric and symmetric equations for the line through $(2, -1, 4)$ and $(5, 3, 2)$. Where does it cross the $xy$-plane?
- Find the equation of the plane that contains the point $(1, 2, -1)$ and the line $x = 3t$, $y = 1 + t$, $z = 2 - t$.
- Decide whether the lines $\langle 1, 2, 0 \rangle + t\langle 1, -1, 2 \rangle$ and $\langle 0, 0, 1 \rangle + s\langle 2, 1, 1 \rangle$ are parallel, intersecting, or skew. Then find either their point of intersection or the distance between them.
- Interview-style: Write the core of a ray tracer's ray-plane test: given a ray origin $O$, direction $\mathbf{d}$, and a plane $\mathbf{n} \cdot \mathbf{x} = c$, return the hit point or report a miss. What edge cases must you handle? Hint: $t = \frac{c - \mathbf{n} \cdot O}{\mathbf{n} \cdot \mathbf{d}}$. Handle $\mathbf{n} \cdot \mathbf{d} \approx 0$ (parallel ray) and $t < 0$ (plane behind the origin), and consider a small $t_{\min}$ to avoid a ray re-hitting the surface it just left.
Cylinders and Quadric Surfaces
Theory
Planes are the 3D version of lines: their equations are first degree. The next step up, second-degree equations in $x, y, z$, gives the quadric surfaces, the 3D relatives of the conic sections of Section 10.3. They appear as satellite dishes (paraboloids), cooling towers (hyperboloids), and the saddle-shaped surfaces of Section 12.6.
Strategy: traces. A surface is hard to draw directly, so slice it. The curve in which a surface meets a plane is a trace. Cutting with horizontal planes $z = k$ is like reading the contour lines on a topographic map, and cutting with $x = k$ and $y = k$ gives vertical cross-sections. For quadrics, every trace is a conic, so a few slices reveal the shape.
Cylinders. If one variable is missing from an equation, that variable is free (Section 11.1, first topic). The surface is a cylinder: take the curve in the plane of the two variables that appear and slide it along the axis of the missing one. The curve is like a cookie-cutter shape and the cylinder is the dough it stamps out.
- $x^2 + y^2 = 1$: a circular cylinder along the $z$-axis.
- $z = x^2$: a parabolic "trough" running along the $y$-axis.
- $y^2 - z^2 = 1$: a hyperbolic cylinder along the $x$-axis.
The six quadrics. After moving and rotating axes, every nondegenerate quadric that is not a cylinder has one of these standard forms ($a, b, c > 0$):
| Surface | Equation | Horizontal traces ($z = k$) | Vertical traces | Picture |
|---|---|---|---|---|
| Ellipsoid | $\frac{x^2}{a^2} + \frac{y^2}{b^2} + \frac{z^2}{c^2} = 1$ | Ellipses | Ellipses | A stretched sphere |
| Elliptic paraboloid | $\frac{z}{c} = \frac{x^2}{a^2} + \frac{y^2}{b^2}$ | Ellipses | Parabolas opening up | A bowl or dish |
| Hyperbolic paraboloid | $\frac{z}{c} = \frac{y^2}{b^2} - \frac{x^2}{a^2}$ | Hyperbolas (two lines at $k = 0$) | Parabolas, some up and some down | A saddle or potato chip |
| Elliptic cone | $\frac{z^2}{c^2} = \frac{x^2}{a^2} + \frac{y^2}{b^2}$ | Ellipses (a point at $k = 0$) | Hyperbolas (crossed lines through the axis) | Two cones tip to tip |
| Hyperboloid of one sheet | $\frac{x^2}{a^2} + \frac{y^2}{b^2} - \frac{z^2}{c^2} = 1$ | Ellipses, smallest at $z = 0$ | Hyperbolas | A cooling tower |
| Hyperboloid of two sheets | $-\frac{x^2}{a^2} - \frac{y^2}{b^2} + \frac{z^2}{c^2} = 1$ | Ellipses for $\lvert k \rvert > c$, empty for $\lvert k \rvert < c$ | Hyperbolas | Two separate bowls |
The axis of the surface is the variable that "behaves differently": the one appearing to the first power in a paraboloid, or the one with the odd sign in a hyperboloid or cone. Counting signs distinguishes the hyperboloids: in the form "$\ldots = 1$", one minus sign gives one sheet, and two minus signs give two sheets.
Hyperbolic paraboloid z = y^2 - x^2 Hyperboloid of one sheet x^2 + y^2 - z^2 = 1
up along y \ /
\ _____ / \ /
\/ \/ \_______/ ellipse at z = 1
/\ /\ <- saddle point at origin ( ) narrowest at z = 0
down along x / \
/ \
slice x = 0: z = y^2 (opens up) / \
slice y = 0: z = -x^2 (opens down) traces in vertical planes are hyperbolas
slice z = 0: y = +-x (two lines)
Identifying a quadric. For an equation such as $4x^2 + y^2 - z^2 - 8x + 4y + 4 = 0$, complete the square in every variable that appears both squared and to the first power, then divide to get $1$ on the right (or $0$ for a cone). Here $4(x - 1)^2 + (y + 2)^2 - z^2 = 4$, so $(x - 1)^2 + \frac{(y + 2)^2}{4} - \frac{z^2}{4} = 1$: a hyperboloid of one sheet, centered at $(1, -2, 0)$, with axis parallel to the $z$-axis.
A remarkable fact: the hyperboloid of one sheet and the hyperbolic paraboloid are doubly ruled. Through every point pass two straight lines lying entirely on the surface. Cooling towers are built from straight steel beams for this reason, and the saddle $z = xy$ is swept out by lines $x = k$, $z = ky$.
Examples
When a cross term like $xy$ appears, the axes of the surface are rotated, and sign counting no longer works directly. Writing the quadratic part as $\mathbf{x}^T M \mathbf{x}$ with a symmetric matrix $M$, the eigenvalues of $M$ are the coefficients in the rotated coordinates (a result from linear algebra). The code below uses this to classify central quadrics $\mathbf{x}^T M \mathbf{x} = k$.
import numpy as np
def classify_central(M, k):
"""Classify the surface x^T M x = k for a symmetric 3x3 matrix M (no linear terms)."""
lam = np.linalg.eigvalsh(M) # real eigenvalues = coefficients after rotating axes
if k == 0:
pos, neg = np.sum(lam > 1e-12), np.sum(lam < -1e-12)
if pos and neg:
return "elliptic cone" if pos + neg == 3 else "pair of planes"
return "single point or line (degenerate)"
s = lam / k # rewrite as sum (lam_i / k) u_i^2 = 1
pos, neg, zero = np.sum(s > 1e-12), np.sum(s < -1e-12), np.sum(np.abs(s) <= 1e-12)
if zero:
return "cylinder over a conic" # a missing variable after rotation
# number of positive coefficients in "... = 1" form decides the type
return {3: "ellipsoid", 2: "hyperboloid of one sheet", 1: "hyperboloid of two sheets", 0: "empty set"}[pos]
print(classify_central(np.diag([1/4, 1/9, 1]), 1)) # ellipsoid x^2/4 + y^2/9 + z^2 = 1
print(classify_central(np.diag([1, 1, -1]), 1)) # hyperboloid of one sheet x^2 + y^2 - z^2 = 1
print(classify_central(np.diag([1, 1, -1]), -1)) # hyperboloid of two sheets x^2 + y^2 - z^2 = -1
print(classify_central(np.diag([1, 1, -1]), 0)) # elliptic cone x^2 + y^2 = z^2
print(classify_central(np.diag([1, 4, 0]), 4)) # cylinder over a conic x^2 + 4y^2 = 4
# A cross term hides the shape: 2xy + z^2 = 1 has M = [[0,1,0],[1,0,0],[0,0,1]]
M = np.array([[0, 1, 0], [1, 0, 0], [0, 0, 1]])
print(np.linalg.eigvalsh(M), classify_central(M, 1)) # [-1. 1. 1.] hyperboloid of one sheet
# Traces of the saddle z = y^2 - x^2 in horizontal planes z = k
for k in [-1, 0, 1]:
kind = "hyperbola opening along x" if k < 0 else ("hyperbola opening along y" if k > 0 else "two lines y = +-x")
print(f"z = {k:2d}: y^2 - x^2 = {k:2d} -> {kind}")
# z = -1: y^2 - x^2 = -1 -> hyperbola opening along x
# z = 0: y^2 - x^2 = 0 -> two lines y = +-x
# z = 1: y^2 - x^2 = 1 -> hyperbola opening along y
Key Takeaways
- A missing variable means a cylinder: the curve in the other two variables, extended along the missing axis.
- Sketch surfaces by their traces in planes $x = k$, $y = k$, $z = k$. The traces of quadrics are conics.
- Know the six standard quadrics. The odd variable gives the axis, and the number of minus signs separates one-sheet from two-sheet hyperboloids.
- Complete the square to find the center. With cross terms, the eigenvalues of the symmetric coefficient matrix give the type.
🧪 Practice
- Identify and sketch: (a) $z = 4 - x^2$; (b) $x^2 + 4z^2 = 16$; (c) $y = x^2 + z^2$; (d) $x^2 - y^2 - z^2 = 1$.
- Complete the square to identify $x^2 + 2y^2 + z^2 - 2x + 8y + 5 = 0$ and $x^2 - y^2 + z^2 + 4y - 4 = 0$. Give the center and axis of each.
- Find the traces of $z^2 = x^2 + 4y^2$ in the planes $z = 2$, $x = 0$, and $y = 1$, and use them to sketch the surface.
- Interview-style: A lidar scanner returns points that should lie on a curved reflector, and you need to decide whether the reflector is a paraboloid or part of an ellipsoid. How would you fit a quadric to the point cloud and classify it? Hint: the general quadric is linear in its 10 coefficients, so fitting is a homogeneous least-squares problem (use the SVD). Then examine the eigenvalues of the $3 \times 3$ quadratic part: a near-zero eigenvalue signals a paraboloid or cylinder.
<a id="112-vector-valued-functions"></a>
11.2 Vector-Valued Functions
A vector-valued function outputs a vector instead of a number, and its arrow tip traces a curve in space. This subchapter extends limits, derivatives, and integrals to such functions. The rule that makes it all easy: work one component at a time.
Definition and Domain
Theory
The functions of earlier chapters take a number and return a number. To describe a moving object in space, we need a function that takes a time $t$ and returns a position: a point, or equivalently its position vector. A vector-valued function (or vector function) is a rule
$$\mathbf{r}(t) = \langle f(t), g(t), h(t) \rangle = f(t)\mathbf{i} + g(t)\mathbf{j} + h(t)\mathbf{k},$$
where the component functions $f, g, h$ are ordinary real functions. This is the parametric description of Section 10.1 with a third coordinate, written as a single vector so that vector algebra can be used on it.
An analogy: picture a flashlight beam anchored at the origin. As time passes, the beam swings and changes length, and the spot where it ends traces a path through the room. The function $\mathbf{r}$ is the beam; the curve is the path the spot traces.
Domain. Unless it is stated otherwise, the domain of $\mathbf{r}$ is the set of all $t$ for which every component is defined, that is, the intersection of the component domains. A single undefined component makes the whole vector undefined, just as a 3D point with an unknown height is not a point.
| Component | Typical restriction |
|---|---|
| $\sqrt{u(t)}$ | $u(t) \ge 0$ |
| $\ln u(t)$ | $u(t) > 0$ |
| $\frac{p(t)}{q(t)}$ | $q(t) \ne 0$ |
| $\arcsin u(t)$, $\arccos u(t)$ | $-1 \le u(t) \le 1$ |
| $\tan u(t)$ | $u(t) \ne \frac{\pi}{2} + k\pi$ |
In applications, the domain is often restricted further by the problem itself: a projectile's trajectory only makes sense from launch until it hits the ground.
Vector functions in two dimensions are written $\mathbf{r}(t) = \langle f(t), g(t) \rangle$. Everything in this subchapter works in any number of dimensions except the cross product, which needs three.
Examples
Find the domain of
$$\mathbf{r}(t) = \left\langle \sqrt{4 - t},\ \ln t,\ \frac{1}{t - 2} \right\rangle.$$
The components require $t \le 4$, $t > 0$, and $t \ne 2$, respectively. The domain is the intersection:
sqrt(4 - t): <=========================================] t <= 4
ln t: (=======================================> t > 0
1/(t - 2): <=====================)(=======================> t != 2
---+----+----+----+----+----+----+----+----+---
-1 0 1 2 3 4 5
domain: (=============)(========] (0, 2) U (2, 4]
import numpy as np
def r(t):
t = np.asarray(t, dtype=float)
with np.errstate(invalid="ignore", divide="ignore"): # silence warnings outside the domain
return np.stack([np.sqrt(4 - t), np.log(t), 1 / (t - 2)], axis=-1)
def in_domain(t):
"""Intersection of the component domains: t <= 4, t > 0, t != 2."""
t = np.asarray(t, dtype=float)
return (t <= 4) & (t > 0) & (t != 2)
ts = np.array([-1.0, 0.0, 0.5, 2.0, 3.0, 4.0, 5.0])
print(in_domain(ts)) # [False False True False True True False]
print(np.isfinite(r(ts)).all(axis=1)) # same pattern: nan or inf in some component outside the domain
print(r(np.array([1.0, 3.0])))
# [[ 1.73205081 0. -1. ] r(1) = <sqrt 3, 0, -1>
# [ 1. 1.09861229 1. ]] r(3) = <1, ln 3, 1>
# Sampled numerically, a vector function is an (n, 3) array: one row per time
t = np.linspace(0, 2 * np.pi, 5)
helix = np.column_stack([np.cos(t), np.sin(t), t / (2 * np.pi)])
print(np.round(helix, 3) + 0.0) # + 0.0 turns -0.0 into 0.0 for display
# [[ 1. 0. 0. ]
# [ 0. 1. 0.25]
# [-1. 0. 0.5 ]
# [ 0. -1. 0.75]
# [ 1. 0. 1. ]] one turn around the cylinder, rising by 1
Key Takeaways
- A vector function $\mathbf{r}(t) = \langle f(t), g(t), h(t) \rangle$ assigns a vector to each value of the parameter. Its tip traces a curve.
- Its domain is the intersection of the component domains: every component must be defined.
- A vector function is a parametric curve written as one vector, so vector algebra applies to it.
🧪 Practice
- Find the domain of $\mathbf{r}(t) = \left\langle \ln(9 - t^2),\ \frac{t}{\sqrt{t + 1}},\ \arcsin\frac{t}{2} \right\rangle$.
- For $\mathbf{r}(t) = \langle t, \cos t, \sin t \rangle$, compute $\mathbf{r}(0)$, $\mathbf{r}(\pi/2)$, $\mathbf{r}(\pi)$, and show that $\lvert \mathbf{r}(t) - t\mathbf{i} \rvert = 1$ for all $t$. What does that say about the curve?
- Write a vector function for the segment from $(1, 2, 3)$ to $(4, 0, -1)$, traversed in 2 seconds at constant velocity, and state its domain.
- Interview-style: You store a robot's trajectory as samples $\mathbf{r}(t_k)$ in an $(n, 3)$ array, but the sensor occasionally drops one coordinate and writes
NaN. How do you find the valid parts of the trajectory and estimate the missing values? Hint: a row is valid only if all three components are finite, which is the domain rule. Interpolate each component separately over time, because component functions are independent.
Limits and Continuity of Vector Functions
Theory
To build calculus, we first need to say what it means for $\mathbf{r}(t)$ to approach a vector $\mathbf{L}$. The natural geometric idea is that the arrow $\mathbf{r}(t)$ ends up arbitrarily close to the arrow $\mathbf{L}$, that is, the distance $\lvert \mathbf{r}(t) - \mathbf{L} \rvert$ goes to $0$:
$$\lim_{t \to a}\mathbf{r}(t) = \mathbf{L} \quad\text{means}\quad \lim_{t \to a}\lvert \mathbf{r}(t) - \mathbf{L} \rvert = 0.$$
Limits are computed componentwise. If $\mathbf{r}(t) = \langle f(t), g(t), h(t) \rangle$, then
$$\lim_{t \to a}\mathbf{r}(t) = \left\langle \lim_{t \to a}f(t),\ \lim_{t \to a}g(t),\ \lim_{t \to a}h(t) \right\rangle,$$
and the left side exists exactly when all three limits on the right exist. Why: write $\mathbf{L} = \langle L_1, L_2, L_3 \rangle$. One component of a vector is never longer than the whole vector, and the length is never more than the sum of the components:
$$\lvert f(t) - L_1 \rvert ;\le; \lvert \mathbf{r}(t) - \mathbf{L} \rvert ;\le; \lvert f(t) - L_1 \rvert + \lvert g(t) - L_2 \rvert + \lvert h(t) - L_3 \rvert.$$
The left inequality (and its versions for $g$ and $h$) shows that if the distance goes to $0$, every component converges. The right inequality shows the converse, by the Squeeze Theorem (Section 2.2).
So everything from Chapter 2 carries over: limit laws, L'Hôpital's rule, and the squeeze theorem are applied one component at a time. The limit laws also extend to vector operations: if $\mathbf{u}(t) \to \mathbf{A}$ and $\mathbf{v}(t) \to \mathbf{B}$, then $\mathbf{u} \cdot \mathbf{v} \to \mathbf{A} \cdot \mathbf{B}$ and $\mathbf{u} \times \mathbf{v} \to \mathbf{A} \times \mathbf{B}$, because both products are built from sums and products of components.
Continuity. $\mathbf{r}$ is continuous at $a$ if $\lim_{t \to a}\mathbf{r}(t) = \mathbf{r}(a)$. By the componentwise rule, $\mathbf{r}$ is continuous exactly when all of its components are continuous. Geometrically, a continuous vector function traces a curve without jumps: the moving point never teleports. One misbehaving component is enough to break continuity.
Examples
Compute $\displaystyle\lim_{t \to 0}\left\langle \frac{\sin t}{t},\ \frac{1 - \cos t}{t},\ e^t \right\rangle$. The first two components are $\frac{0}{0}$ forms whose limits are $1$ and $0$ (Section 2.2), and $e^t \to 1$ by continuity. So the limit is $\langle 1, 0, 1 \rangle$. Note that $\mathbf{r}(0)$ itself is undefined, so the function has a removable discontinuity at $0$. Defining $\mathbf{r}(0) = \langle 1, 0, 1 \rangle$ removes it.
On the other hand, $\mathbf{q}(t) = \left\langle t,\ t^2,\ \frac{\lvert t \rvert}{t} \right\rangle$ has no limit at $0$: the first two components behave well, but the third tends to $1$ from the right and $-1$ from the left. The point jumps from $(0, 0, -1)$ to $(0, 0, 1)$ as $t$ passes $0$.
import sympy as sp
import numpy as np
t = sp.symbols('t')
r = sp.Matrix([sp.sin(t) / t, (1 - sp.cos(t)) / t, sp.exp(t)])
print([sp.limit(c, t, 0) for c in r]) # [1, 0, 1] one ordinary limit per component
# Numerically: the distance |r(t) - L| shrinks to 0 (roughly in proportion to t here)
L = np.array([1.0, 0.0, 1.0])
for h in [1e-1, 1e-2, 1e-3]:
v = np.array([np.sin(h) / h, (1 - np.cos(h)) / h, np.exp(h)])
print(h, np.linalg.norm(v - L))
# 0.1 0.11644540981575136
# 0.01 0.011225226925597792
# 0.001 0.001118481367608234
# A jump in one component destroys the limit of the whole vector function
q = sp.Matrix([t, t**2, sp.Abs(t) / t])
print(sp.limit(q[2], t, 0, '+'), sp.limit(q[2], t, 0, '-')) # 1 -1 one-sided limits disagree
Key Takeaways
- $\mathbf{r}(t) \to \mathbf{L}$ means $\lvert \mathbf{r}(t) - \mathbf{L} \rvert \to 0$, which is equivalent to every component converging.
- Compute limits one component at a time with the tools of Chapter 2. If one component limit fails to exist, the vector limit fails to exist.
- $\mathbf{r}$ is continuous exactly when all components are continuous. Limits pass through sums, scalar multiples, dot products, and cross products.
🧪 Practice
- Evaluate $\displaystyle\lim_{t \to 1}\left\langle \frac{t^2 - 1}{t - 1},\ \sqrt{t + 3},\ \frac{\ln t}{t - 1} \right\rangle$ and $\displaystyle\lim_{t \to \infty}\left\langle e^{-t},\ \frac{2t^2 + 1}{t^2 - t},\ t\sin\frac{1}{t} \right\rangle$.
- Where is $\mathbf{r}(t) = \left\langle \lfloor t \rfloor,\ \frac{1}{t^2 + 1},\ \tan t \right\rangle$ discontinuous? Classify each discontinuity as a jump or an infinite discontinuity.
- Prove that if $\mathbf{u}(t) \to \mathbf{A}$ and $\mathbf{v}(t) \to \mathbf{B}$, then $\mathbf{u}(t) \cdot \mathbf{v}(t) \to \mathbf{A} \cdot \mathbf{B}$, using the componentwise rule and the limit laws.
- Interview-style: A test asserts that a simulated trajectory converges to a target point by checking that each coordinate is within $10^{-6}$ of the target. A teammate proposes checking $\lvert \mathbf{r} - \mathbf{L} \rvert < 10^{-6}$ instead. Are the two checks equivalent, and which tolerance in one corresponds to a given tolerance in the other? Hint: use $\max_i \lvert r_i - L_i \rvert \le \lvert \mathbf{r} - \mathbf{L} \rvert \le \sqrt{3}\max_i \lvert r_i - L_i \rvert$.
Derivatives of Vector Functions
Theory
For a real function, the derivative is the limit of the slopes of secant lines. For a vector function, the analogous object is the limit of secant vectors. As $h$ shrinks, the vector $\mathbf{r}(t + h) - \mathbf{r}(t)$ points from one point on the curve to a nearby one. Dividing by $h$ keeps it from shrinking to $\mathbf{0}$, and in the limit, it points along the curve:
$$\mathbf{r}'(t) = \lim_{h \to 0}\frac{\mathbf{r}(t + h) - \mathbf{r}(t)}{h}.$$
* r(t + h)
.'
secant .' (r(t+h) - r(t)) / h
vectors .' rotates toward the tangent
.' as h -> 0
--------- *==============> r'(t): tangent vector
.' r(t)
.'
.' the curve
O (origin; r(t) and r(t+h) are arrows from here)
By the componentwise rule for limits, derivatives are computed componentwise:
$$\mathbf{r}'(t) = \langle f'(t), g'(t), h'(t) \rangle.$$
Interpretation. If $\mathbf{r}(t)$ is a position at time $t$, then $\mathbf{r}'(t)$ is the velocity: its direction is the direction of motion, and its length $\lvert \mathbf{r}'(t) \rvert$ is the speed (Section 11.4). Geometrically, $\mathbf{r}'(t)$ is the tangent vector at the point $\mathbf{r}(t)$, provided it is not $\mathbf{0}$. The tangent line at $t_0$ is $\mathbf{L}(s) = \mathbf{r}(t_0) + s,\mathbf{r}'(t_0)$.
A curve is smooth on an interval if $\mathbf{r}'$ is continuous and never $\mathbf{0}$ there. As in the plane, a zero derivative means the point stops, and it can leave in a new direction, producing a cusp or corner. For example, $\mathbf{c}(t) = \langle t^3, t^2, 0 \rangle$ has the perfectly smooth component functions $t^3$ and $t^2$, but $\mathbf{c}'(0) = \mathbf{0}$, and the curve $y = x^{2/3}$ it traces has a cusp at the origin.
Differentiation rules. Let $\mathbf{u}, \mathbf{v}$ be differentiable vector functions, $c$ a scalar, and $f$ a real function.
| Rule | Formula |
|---|---|
| Sum | $(\mathbf{u} + \mathbf{v})' = \mathbf{u}' + \mathbf{v}'$ |
| Constant multiple | $(c\mathbf{u})' = c\mathbf{u}'$ |
| Scalar function times vector | $(f\mathbf{u})' = f'\mathbf{u} + f\mathbf{u}'$ |
| Dot product | $(\mathbf{u} \cdot \mathbf{v})' = \mathbf{u}' \cdot \mathbf{v} + \mathbf{u} \cdot \mathbf{v}'$ |
| Cross product | $(\mathbf{u} \times \mathbf{v})' = \mathbf{u}' \times \mathbf{v} + \mathbf{u} \times \mathbf{v}'$ (keep the order) |
| Chain rule | $\frac{d}{dt}\mathbf{u}(f(t)) = f'(t),\mathbf{u}'(f(t))$ |
Each rule follows from the ordinary product rule (Section 3.3) applied to components. The three product rules all have the familiar "derivative of the first times the second, plus the first times derivative of the second" shape. The only trap is that the cross product is anticommutative, so the factors must stay in their original order.
A key consequence: constant length means perpendicular derivative. If $\lvert \mathbf{r}(t) \rvert = c$ for all $t$, then $\mathbf{r}(t) \cdot \mathbf{r}'(t) = 0$. Proof: differentiate $\mathbf{r} \cdot \mathbf{r} = c^2$ with the dot product rule to get $2,\mathbf{r} \cdot \mathbf{r}' = 0$.
Geometrically, a point moving on a sphere centered at the origin always moves tangent to the sphere, perpendicular to the radius. This small fact is used repeatedly: it shows that the unit tangent $\mathbf{T}$ is perpendicular to its own derivative (Section 11.3), and that uniform circular motion has acceleration pointing to the center (Section 11.4).
Higher derivatives are also componentwise: $\mathbf{r}''(t) = \langle f''(t), g''(t), h''(t) \rangle$ is the acceleration.
Examples
For the helix $\mathbf{r}(t) = \langle \cos t, \sin t, t \rangle$,
$$\mathbf{r}'(t) = \langle -\sin t, \cos t, 1 \rangle, \qquad \lvert \mathbf{r}'(t) \rvert = \sqrt{\sin^2 t + \cos^2 t + 1} = \sqrt{2}.$$
The speed is constant, and since $\mathbf{r}'$ is never zero, the helix is smooth. At $t_0 = \frac{\pi}{2}$, the point is $(0, 1, \frac{\pi}{2})$ and the tangent vector is $\langle -1, 0, 1 \rangle$, so the tangent line is $\mathbf{L}(s) = \langle -s,\ 1,\ \frac{\pi}{2} + s \rangle$.
import sympy as sp
t = sp.symbols('t', real=True)
r = sp.Matrix([sp.cos(t), sp.sin(t), t]) # circular helix
rp = r.diff(t) # differentiates each component
print(rp.T, sp.simplify(rp.norm())) # Matrix([[-sin(t), cos(t), 1]]) sqrt(2)
# Tangent line at t = pi/2: point + s * tangent vector
t0, s = sp.pi / 2, sp.symbols('s')
print((r.subs(t, t0) + s * rp.subs(t, t0)).T) # Matrix([[-s, 1, s + pi/2]])
# Check the product rules on two arbitrary vector functions
u = sp.Matrix([t**2, sp.exp(t), sp.sin(t)])
v = sp.Matrix([1 / (1 + t**2), t, sp.cos(t)])
dot_rule = sp.diff(u.dot(v), t) - (u.diff(t).dot(v) + u.dot(v.diff(t)))
cross_rule = u.cross(v).diff(t) - (u.diff(t).cross(v) + u.cross(v.diff(t)))
swapped = u.cross(v).diff(t) - (v.cross(u.diff(t)) + u.cross(v.diff(t))) # wrong order in the first term
print(sp.simplify(dot_rule), sp.simplify(cross_rule).T) # 0 Matrix([[0, 0, 0]]) both rules hold
print(swapped.subs(t, 1).evalf(4).T) # Matrix([[1.857, -1.621, 1.282]]) not zero
# Constant length implies r . r' = 0
w = sp.Matrix([sp.cos(t) * sp.sin(t), sp.sin(t)**2, sp.cos(t)])
print(sp.simplify(w.dot(w)), sp.simplify(w.dot(w.diff(t)))) # 1 0 on the unit sphere, so w is perpendicular to w'
# Smooth components, but not a smooth curve: c'(0) = 0 (a cusp)
c = sp.Matrix([t**3, t**2, 0])
print(c.diff(t).subs(t, 0).T) # Matrix([[0, 0, 0]])
Key Takeaways
- $\mathbf{r}'(t)$ is the limit of secant vectors and is computed componentwise. It is the tangent vector, or the velocity.
- A curve is smooth where $\mathbf{r}'$ is continuous and nonzero. Smooth components alone do not guarantee a smooth curve.
- Sum, scalar, dot, cross, and chain rules all hold. Keep the factor order in the cross product rule.
- If $\lvert \mathbf{r}(t) \rvert$ is constant, then $\mathbf{r} \cdot \mathbf{r}' = 0$.
🧪 Practice
- Find $\mathbf{r}'(t)$ and $\mathbf{r}''(t)$ for $\mathbf{r}(t) = \langle t e^{-t},\ \arctan t,\ \ln(1 + t^2) \rangle$, and find the tangent line at $t = 0$.
- Find the points where $\mathbf{r}(t) = \langle t^3 - 3t,\ t^2,\ t^3 \rangle$ fails to be smooth, if any, and the points where its tangent vector is parallel to the $xy$-plane.
- Show that $\frac{d}{dt}\big(\mathbf{r} \times \mathbf{r}'\big) = \mathbf{r} \times \mathbf{r}''$ and $\frac{d}{dt}\big[\mathbf{u} \cdot (\mathbf{v} \times \mathbf{w})\big] = \mathbf{u}' \cdot (\mathbf{v} \times \mathbf{w}) + \mathbf{u} \cdot (\mathbf{v}' \times \mathbf{w}) + \mathbf{u} \cdot (\mathbf{v} \times \mathbf{w}')$.
- Interview-style: A simulation rotates a vector with explicit Euler steps, $\mathbf{r} \leftarrow \mathbf{r} + \Delta t,(\boldsymbol{\omega} \times \mathbf{r})$. After one revolution with $\Delta t = 0.01$, the vector's length has grown from $1$ to about $1.032$. Explain why, and suggest a fix. Hint: the exact flow keeps $\mathbf{r} \cdot \mathbf{r}' = 0$, but a straight step along the tangent always lands outside the circle: each step multiplies the length by $\sqrt{1 + (\lvert\boldsymbol{\omega}\rvert\Delta t)^2}$. Renormalize after each step, or apply an exact rotation (rotation matrix or quaternion) instead.
Integrals of Vector Functions
Theory
Differentiation turns position into velocity. The reverse question comes up constantly: a navigation system measures velocity (or acceleration) and must reconstruct position. That is integration, and like everything else in this subchapter, it is done componentwise.
Definite integral. Define $\int_a^b \mathbf{r}(t),dt$ as a limit of Riemann sums $\sum \mathbf{r}(t_i^*),\Delta t$, exactly as in Section 5.2. Each Riemann sum is a vector whose components are Riemann sums of $f$, $g$, $h$, so
$$\int_a^b \mathbf{r}(t),dt = \left\langle \int_a^b f(t),dt,\ \int_a^b g(t),dt,\ \int_a^b h(t),dt \right\rangle.$$
Fundamental Theorem of Calculus. If $\mathbf{R}' = \mathbf{r}$ on $[a, b]$ (that is, $\mathbf{R}$ is an antiderivative of $\mathbf{r}$), then
$$\int_a^b \mathbf{r}(t),dt = \mathbf{R}(b) - \mathbf{R}(a).$$
Indefinite integral. Two antiderivatives of the same vector function differ by a constant vector, so
$$\int \mathbf{r}(t),dt = \mathbf{R}(t) + \mathbf{C}, \qquad \mathbf{C} = \langle C_1, C_2, C_3 \rangle.$$
There is one arbitrary constant per component, and one initial condition such as $\mathbf{r}(0) = \mathbf{r}_0$ fixes all three at once. Physically: velocity determines the shape of a path, and the starting point determines where it is.
Displacement versus distance. If $\mathbf{v}(t)$ is a velocity, then $\int_a^b \mathbf{v}(t),dt = \mathbf{r}(b) - \mathbf{r}(a)$ is the displacement: the net change in position, regardless of the route. The distance traveled is $\int_a^b \lvert \mathbf{v}(t) \rvert,dt$, the integral of the speed, which is the arc length of Section 11.3. A runner who completes a lap has displacement $\mathbf{0}$ but a distance of $400$ m. In general,
$$\left\lvert \int_a^b \mathbf{v}(t),dt \right\rvert \le \int_a^b \lvert \mathbf{v}(t) \rvert,dt,$$
the integral form of the triangle inequality.
Properties such as linearity, splitting the interval, and $\int_a^b \mathbf{c} \cdot \mathbf{r}(t),dt = \mathbf{c} \cdot \int_a^b \mathbf{r}(t),dt$ for a constant vector $\mathbf{c}$ all carry over.
Examples
A particle has velocity $\mathbf{r}'(t) = \langle 2t,\ 3\sqrt{t},\ e^{-t} \rangle$ for $t \ge 0$ and starts at $\mathbf{r}(0) = \langle 1, 0, 2 \rangle$. Integrating each component,
$$\mathbf{r}(t) = \langle t^2,\ 2t^{3/2},\ -e^{-t} \rangle + \mathbf{C}.$$
At $t = 0$, this gives $\langle 0, 0, -1 \rangle + \mathbf{C} = \langle 1, 0, 2 \rangle$, so $\mathbf{C} = \langle 1, 0, 3 \rangle$ and
$$\mathbf{r}(t) = \langle t^2 + 1,\ 2t^{3/2},\ 3 - e^{-t} \rangle.$$
A definite integral: $\displaystyle\int_0^{\pi/2}\langle \cos t, \sin t, t^2 \rangle,dt = \left\langle 1,\ 1,\ \frac{\pi^3}{24} \right\rangle$.
import numpy as np
import sympy as sp
from scipy.integrate import quad
t = sp.symbols('t', real=True)
rp = sp.Matrix([2 * t, 3 * sp.sqrt(t), sp.exp(-t)]) # given velocity r'(t)
R = rp.integrate((t, 0, t)) # integral from 0 to t, component by component
r = sp.Matrix([1, 0, 2]) + R # r(t) = r(0) + integral_0^t r'(u) du
print(sp.simplify(r).T) # Matrix([[t**2 + 1, 2*t**(3/2), 3 - exp(-t)]])
print([float(x) for x in r.subs(t, 1)]) # [2.0, 2.0, 2.6321205588285577] position at t = 1
# Definite integral of a vector function, three ways
f = lambda s: np.array([np.cos(s), np.sin(s), s**2])
exact = np.array([1.0, 1.0, (np.pi / 2)**3 / 3])
by_quad = np.array([quad(lambda s, k=k: f(s)[k], 0, np.pi / 2)[0] for k in range(3)]) # one call per component
s = np.linspace(0, np.pi / 2, 101)
by_trapz = np.trapezoid(f(s), s, axis=1) # all components at once
print(exact, by_quad, by_trapz, sep="\n")
# [1. 1. 1.2919282]
# [1. 1. 1.2919282]
# [0.99997944 0.99997944 1.29199279] trapezoid rule, 100 panels: error about 2e-5
Key Takeaways
- Integrate vector functions componentwise. The Fundamental Theorem of Calculus holds unchanged.
- The constant of integration is a constant vector, and one initial condition determines it.
- $\int \mathbf{v},dt$ gives displacement, and $\int \lvert \mathbf{v} \rvert,dt$ gives distance traveled. The first is never longer than the second.
🧪 Practice
- Evaluate $\displaystyle\int_0^1 \left\langle \frac{1}{1 + t^2},\ t e^{t^2},\ \sqrt{t} \right\rangle dt$.
- Find $\mathbf{r}(t)$ if $\mathbf{r}'(t) = \langle \sin 2t,\ \sec^2 t,\ 6t \rangle$ and $\mathbf{r}(0) = \langle 0, 1, -1 \rangle$.
- A particle moves with velocity $\mathbf{v}(t) = \langle \cos t, \sin t, 0 \rangle$ for $0 \le t \le 2\pi$. Compute its displacement and the distance it travels, and explain the difference.
- Interview-style: An inertial navigation unit integrates measured acceleration twice to estimate position. A tiny constant bias $\mathbf{b}$ in the accelerometer turns into a large position error after a few minutes. How fast does the error grow, and how do real systems keep it bounded? Hint: integrate $\mathbf{b}$ twice from rest to get an error of $\frac{1}{2}\mathbf{b}t^2$. Systems fuse in absolute measurements (GPS, landmarks) to reset the drift.
Space Curves and Their Graphs
Theory
The set of points traced by a continuous vector function $\mathbf{r}(t)$ is a space curve. As in Section 10.1, the curve is the picture, and the vector function is one particular trip along it, with an orientation and timing. Drawing space curves is harder than drawing plane curves, because a 2D sketch flattens out the depth. Three techniques help.
1. Recognize the surface the curve lives on. If the components satisfy an equation, the curve lies on that surface. For the circular helix $\mathbf{r}(t) = \langle a\cos t, a\sin t, ct \rangle$, $x^2 + y^2 = a^2$, so the curve winds around a cylinder of radius $a$, rising $2\pi c$ per turn (the pitch), like the thread of a screw or a spiral staircase. If $c > 0$, it is a right-handed helix, which is the shape of DNA and of most screws.
z
| .-------. one turn rises by 2 pi c
| .' ^ '.
| ( | ) <- t = 2 pi
| '. | .'
| '---|---'
| .-----|-.
| .' | '.
| ( | ) <- t = pi, on the far side
| '. | .'
| '-----|-'
| .--------|.
|.' | '.
( * ) <- t = 0 at (a, 0, 0)
'. .'
'--------' shadow on the xy-plane: the circle x^2 + y^2 = a^2
2. Look at the shadows. Dropping one coordinate projects the curve onto a coordinate plane, giving a plane curve that is usually easy to identify. The twisted cubic $\mathbf{r}(t) = \langle t, t^2, t^3 \rangle$ has shadows $y = x^2$ (a parabola) in the $xy$-plane, $z = x^3$ (a cubic) in the $xz$-plane, and $z^2 = y^3$ (a cusp) in the $yz$-plane. Putting the three views together, as in an engineering drawing, reveals the 3D shape.
3. Parametrize an intersection of surfaces. A curve defined as the intersection of two surfaces (two equations) can often be parametrized by solving one equation in a natural way and substituting into the other.
- Cylinder and plane. On $x^2 + y^2 = 1$, set $x = \cos t$, $y = \sin t$. The plane $y + z = 2$ then forces $z = 2 - \sin t$. The intersection is $\mathbf{r}(t) = \langle \cos t, \sin t, 2 - \sin t \rangle$, $0 \le t \le 2\pi$: an ellipse, tilted with the plane.
- Sphere and cylinder (Viviani's curve). The sphere $x^2 + y^2 + z^2 = 4$ meets the cylinder $(x - 1)^2 + y^2 = 1$. Parametrize the cylinder by $x = 1 + \cos t$, $y = \sin t$. Then $z^2 = 4 - x^2 - y^2 = 2 - 2\cos t = 4\sin^2\frac{t}{2}$, so $z = 2\sin\frac{t}{2}$ for $0 \le t \le 4\pi$ covers the whole figure-eight-shaped curve.
- Graph-type surfaces. On the paraboloid $z = x^2 + y^2$, any plane curve $\langle x(t), y(t) \rangle$ lifts to the space curve $\langle x(t), y(t), x(t)^2 + y(t)^2 \rangle$.
Some notable curves:
| Curve | Vector function | Lies on |
|---|---|---|
| Circular helix | $\langle a\cos t, a\sin t, ct \rangle$ | Cylinder $x^2 + y^2 = a^2$ |
| Conical helix | $\langle t\cos t, t\sin t, t \rangle$ | Cone $z^2 = x^2 + y^2$ |
| Twisted cubic | $\langle t, t^2, t^3 \rangle$ | Parabolic cylinder $y = x^2$ |
| Toroidal spiral | $\langle (R + r\cos nt)\cos t,\ (R + r\cos nt)\sin t,\ r\sin nt \rangle$ | Torus (a donut) |
| Viviani's curve | $\langle 1 + \cos t, \sin t, 2\sin\frac{t}{2} \rangle$ | Sphere of radius $2$ and cylinder |
Is a point on the curve? Unlike a surface equation, a parametrization does not give a direct membership test. Find the $t$ that the point would need from one component (often the simplest), then check the other components with that $t$.
Examples
The code checks each parametrization against the surfaces it should lie on. Residuals at the level of $10^{-16}$ are rounding error, which confirms the derivations.
import numpy as np
t = np.linspace(0, 4 * np.pi, 2001)
# 1) Ellipse cut from the cylinder x^2 + y^2 = 1 by the plane y + z = 2
x, y, z = np.cos(t), np.sin(t), 2 - np.sin(t)
print(np.abs(x**2 + y**2 - 1).max(), np.abs(y + z - 2).max())
# 2.220446049250313e-16 2.220446049250313e-16 on both surfaces, up to rounding
# 2) Viviani's curve: the sphere x^2 + y^2 + z^2 = 4 meets the cylinder (x - 1)^2 + y^2 = 1
x, y, z = 1 + np.cos(t), np.sin(t), 2 * np.sin(t / 2)
print(np.abs(x**2 + y**2 + z**2 - 4).max(), np.abs((x - 1)**2 + y**2 - 1).max())
# 8.881784197001252e-16 4.440892098500626e-16
# 3) Twisted cubic <t, t^2, t^3>: its shadows on the coordinate planes
s = np.array([-1.0, 0.0, 1.0, 2.0])
xs, ys, zs = s, s**2, s**3
print("xy-plane: y = x^2 ->", np.allclose(ys, xs**2)) # True
print("xz-plane: z = x^3 ->", np.allclose(zs, xs**3)) # True
print("yz-plane: z^2 = y^3 ->", np.allclose(zs**2, ys**3)) # True
# 4) Is a point on the helix <cos t, sin t, t>? Recover t from z, then compare all components
def on_helix(P, tol=1e-9):
t0 = P[2] # z = t gives the only candidate parameter
return np.linalg.norm(np.array([np.cos(t0), np.sin(t0), t0]) - P) < tol
print(on_helix(np.array([-1.0, 0.0, np.pi])), on_helix(np.array([1.0, 0.0, np.pi]))) # True False
Key Takeaways
- A space curve is the set of points traced by a continuous vector function. The function adds orientation and timing.
- Identify a curve by finding a surface it lies on and by projecting it onto the coordinate planes.
- Parametrize an intersection of two surfaces by parametrizing one surface naturally and solving the other equation for the remaining coordinate.
- To test whether a point is on a parametrized curve, solve one component for $t$ and check the rest.
🧪 Practice
- Sketch $\mathbf{r}(t) = \langle 2\cos t,\ 2\sin t,\ 4 - t \rangle$ for $0 \le t \le 4\pi$, and describe its projections onto the three coordinate planes. Is $(-2, 0, 4 - \pi)$ on the curve? Is $(0, 2, 0)$?
- Find a vector function for the curve of intersection of (a) the paraboloid $z = 4x^2 + y^2$ and the parabolic cylinder $y = x^2$; (b) the cone $z = \sqrt{x^2 + y^2}$ and the plane $z = 1 + y$.
- Show that the conical helix $\langle t\cos t, t\sin t, t \rangle$ lies on the cone $z^2 = x^2 + y^2$, and find the points where it meets the sphere $x^2 + y^2 + z^2 = 8$.
- Interview-style: Two drones fly along the paths $\mathbf{r}_1(t)$ and $\mathbf{r}_2(t)$. A colleague checks for a collision by testing whether the two space curves intersect. Why is that the wrong test, and what is the right one? Hint: the curves can cross while the drones pass the crossing point at different times. A collision requires $\mathbf{r}_1(t) = \mathbf{r}_2(t)$ at the same $t$, or in practice $\min_t \lvert \mathbf{r}_1(t) - \mathbf{r}_2(t) \rvert$ below a safety radius.
<a id="113-geometry-of-curves"></a>
11.3 Geometry of Curves
A vector function describes both a curve and a particular way of moving along it. This subchapter separates the two. It measures properties that belong to the curve itself, such as length, bending (curvature), and twisting (torsion), and it builds the moving frame $\mathbf{T}, \mathbf{N}, \mathbf{B}$ that encodes them.
Arc Length in Space
Theory
How long is a helical spring, or the cable wound around a drum? Section 7.2 answered this question for plane curves. The same idea works in space without change: approximate the curve by a polygon, add up the lengths of its segments, and let the segments shrink.
Split $[a, b]$ into pieces of width $\Delta t$. Over one piece, the curve moves from $\mathbf{r}(t_{i-1})$ to $\mathbf{r}(t_i)$, a chord of length $\lvert \mathbf{r}(t_i) - \mathbf{r}(t_{i-1}) \rvert \approx \lvert \mathbf{r}'(t_i) \rvert,\Delta t$. Summing and taking the limit gives a Riemann sum for the speed:
$$L = \int_a^b \lvert \mathbf{r}'(t) \rvert,dt = \int_a^b \sqrt{f'(t)^2 + g'(t)^2 + h'(t)^2},dt.$$
This assumes the curve is smooth (or piecewise smooth, in which case the pieces are added) and is traced exactly once as $t$ goes from $a$ to $b$. If the parametrization retraces part of the curve, the integral counts that part twice: it measures distance traveled, which equals the curve's length only for a one-time traversal.
In words, length = integral of speed over time. If you drive for an hour at a varying speed, the odometer reading is the integral of your speedometer reading, whatever the shape of the road.
Independence of parametrization. Two parametrizations that trace the same curve once give the same length. If $t = \varphi(u)$ with $\varphi' > 0$, the chain rule gives $\left\lvert \frac{d}{du}\mathbf{r}(\varphi(u)) \right\rvert = \lvert \mathbf{r}'(\varphi(u)) \rvert,\varphi'(u)$, and the substitution rule turns one integral into the other. Faster travel means less time spent, and the two effects cancel exactly. Length is a property of the curve, not of the trip.
Unrolling the helix. For the helix $\langle a\cos t, a\sin t, ct \rangle$, the speed is $\sqrt{a^2 + c^2}$, constant, so one turn has length $2\pi\sqrt{a^2 + c^2}$. There is a nice way to see this. Cut the cylinder along a vertical line and flatten it into a rectangle of width $2\pi a$. One turn of the helix becomes a straight diagonal that rises $2\pi c$, so its length is the hypotenuse $\sqrt{(2\pi a)^2 + (2\pi c)^2}$.
cylinder cut open and flattened
2 pi c +--------------------------* one turn of the helix becomes the
| .' | straight diagonal of the rectangle:
| .' |
| .' | L = sqrt((2 pi a)^2 + (2 pi c)^2)
| .' | = 2 pi sqrt(a^2 + c^2)
| .' |
| .' |
0 *--------------------------+
0 2 pi a (circumference of the cylinder)
Most arc length integrals have no elementary antiderivative (the square root rarely simplifies), so in practice they are evaluated numerically with the methods of Section 6.4. The exercises designed to be done by hand are ones where the expression under the root is a perfect square.
Examples
Find the length of $\mathbf{r}(t) = \langle 2t, t^2, \ln t \rangle$ for $1 \le t \le e$. Here $\mathbf{r}'(t) = \left\langle 2, 2t, \frac{1}{t} \right\rangle$, and
$$\lvert \mathbf{r}'(t) \rvert^2 = 4 + 4t^2 + \frac{1}{t^2} = \left(2t + \frac{1}{t}\right)^2,$$
a perfect square. So
$$L = \int_1^e \left(2t + \frac{1}{t}\right)dt = \Big[t^2 + \ln t\Big]_1^e = (e^2 + 1) - 1 = e^2 \approx 7.389.$$
import numpy as np
from scipy.integrate import quad
def arc_length(rp, a, b):
"""Integrate the speed |r'(t)| over [a, b]."""
return quad(lambda t: np.linalg.norm(rp(t)), a, b)[0]
# r(t) = <2t, t^2, ln t> on [1, e]: exact length e^2
rp = lambda t: np.array([2.0, 2 * t, 1 / t])
print(arc_length(rp, 1, np.e), np.e**2) # 7.38905609893065 7.3890560989306495
# Helix <3 cos t, 3 sin t, 4t>, one turn: speed 5, length 10 pi
hp = lambda t: np.array([-3 * np.sin(t), 3 * np.cos(t), 4.0])
print(arc_length(hp, 0, 2 * np.pi), 10 * np.pi) # 31.41592653589793 31.41592653589793
# The definition itself: inscribed polygons approach the length from below
r = lambda t: np.array([2 * t, t**2, np.log(t)])
for n in [2, 8, 32, 128]:
pts = np.array([r(t) for t in np.linspace(1, np.e, n + 1)]) # n + 1 vertices
print(n, np.sum(np.linalg.norm(np.diff(pts, axis=0), axis=1))) # sum of chord lengths
# 2 7.369630178224763
# 8 7.387755623588834
# 32 7.388974434250622
# 128 7.3890509933708 error shrinks by about 16x per 4x more segments (O(1/n^2))
# Twisted cubic <t, t^2, t^3> on [0, 1]: sqrt(1 + 4t^2 + 9t^4) has no elementary antiderivative
cp = lambda t: np.array([1.0, 2 * t, 3 * t**2])
print(arc_length(cp, 0, 1)) # 1.8630229825122513
Key Takeaways
- $L = \int_a^b \lvert \mathbf{r}'(t) \rvert,dt$: the integral of speed, the limit of inscribed polygon lengths.
- The curve must be traced once; otherwise, the integral gives distance traveled, not the curve's length.
- Length does not depend on the parametrization.
- Exact answers usually require $\lvert \mathbf{r}' \rvert^2$ to be a perfect square. Otherwise, integrate numerically.
🧪 Practice
- Find the length of $\mathbf{r}(t) = \langle \cos^3 t, \sin^3 t, \cos 2t \rangle$ for $0 \le t \le \frac{\pi}{2}$.
- Show that $\mathbf{r}(t) = \left\langle \sqrt{2},t,\ e^t,\ e^{-t} \right\rangle$ has speed $e^t + e^{-t}$, and use it to show that the length from $t = 0$ to $t = 1$ is $e - e^{-1}$.
- A wire is wound around a cylindrical drum of radius $0.5$ m in a helix that makes $40$ turns while rising $2$ m. Find the length of the wire, and compare it with $40$ times the circumference.
- Interview-style: A GPS logger records a hiker's position once per minute, and the app reports distance as the sum of straight-line segments. Is the reported distance biased high or low, and what makes the bias worse? Is there an opposite effect? Hint: chords are never longer than arcs, so sharp, winding trails are underestimated. But GPS noise adds random zigzags, which inflate the sum, especially at high sampling rates when the hiker barely moves.
Arc Length Parametrization
Theory
A parametrization describes the curve together with a timing. For geometry, the timing is a distraction: a curve's bending should not depend on how fast someone drives along it. The cure is to use distance along the curve as the parameter, the way mile markers label a highway.
The arc length function. Fix a starting time $a$ and define
$$s(t) = \int_a^t \lvert \mathbf{r}'(u) \rvert,du,$$
the distance traveled along the curve from $\mathbf{r}(a)$ to $\mathbf{r}(t)$. By the Fundamental Theorem of Calculus,
$$\frac{ds}{dt} = \lvert \mathbf{r}'(t) \rvert,$$
which says that the rate at which arc length accumulates is the speed. For a smooth curve, $\frac{ds}{dt} > 0$, so $s$ is strictly increasing and has an inverse $t = t(s)$.
Reparametrizing. Substituting $t = t(s)$ gives $\mathbf{r}(t(s))$, the same curve with arc length as the parameter. Its key property is unit speed: by the chain rule,
$$\left\lvert \frac{d\mathbf{r}}{ds} \right\rvert = \lvert \mathbf{r}'(t) \rvert,\left\lvert \frac{dt}{ds} \right\rvert = \lvert \mathbf{r}'(t) \rvert \cdot \frac{1}{\lvert \mathbf{r}'(t) \rvert} = 1.$$
Conversely, a parametrization has unit speed exactly when its parameter measures arc length (from some starting point). The steps are:
- Compute the speed $\lvert \mathbf{r}'(t) \rvert$.
- Integrate to get $s(t) = \int_a^t \lvert \mathbf{r}'(u) \rvert,du$.
- Solve for $t$ in terms of $s$.
- Substitute into $\mathbf{r}(t)$.
Step 3 is the bottleneck. It works by hand only in special cases, which is why arc length parametrizations are mainly a theoretical tool. Definitions such as curvature (next topics) are stated in terms of $s$ because that makes them geometric. Formulas are then translated back to an arbitrary $t$ with $\frac{d}{ds} = \frac{1}{\lvert \mathbf{r}'(t) \rvert}\frac{d}{dt}$.
In computation, $s(t)$ is tabulated numerically, and the table is inverted by interpolation. This is how animation systems move objects along splines at constant speed.
Examples
The helix $\mathbf{r}(t) = \langle \cos t, \sin t, t \rangle$ has speed $\sqrt{2}$, so measuring from $t = 0$, $s = \sqrt{2},t$ and $t = \frac{s}{\sqrt{2}}$:
$$\mathbf{r}(s) = \left\langle \cos\frac{s}{\sqrt{2}},\ \sin\frac{s}{\sqrt{2}},\ \frac{s}{\sqrt{2}} \right\rangle.$$
A less trivial case: $\mathbf{r}(t) = e^t\langle \cos t, \sin t, 1 \rangle$, a spiral on a cone. Its derivative is $e^t\langle \cos t - \sin t,\ \sin t + \cos t,\ 1 \rangle$, and since $(\cos t - \sin t)^2 + (\sin t + \cos t)^2 = 2$, the speed is $\sqrt{3},e^t$. Then
$$s(t) = \int_0^t \sqrt{3},e^u,du = \sqrt{3}(e^t - 1) \quad\Longrightarrow\quad t = \ln\left(1 + \frac{s}{\sqrt{3}}\right),$$
$$\mathbf{r}(s) = \left(1 + \frac{s}{\sqrt{3}}\right)\left\langle \cos\ln\left(1 + \tfrac{s}{\sqrt{3}}\right),\ \sin\ln\left(1 + \tfrac{s}{\sqrt{3}}\right),\ 1 \right\rangle.$$
import numpy as np
# r(t) = e^t <cos t, sin t, 1>: speed sqrt(3) e^t, so s(t) = sqrt(3)(e^t - 1)
def r(t):
return np.exp(t) * np.array([np.cos(t), np.sin(t), 1.0])
def r_of_s(s):
t = np.log(1 + s / np.sqrt(3)) # invert s(t) in closed form
return r(t)
# Unit speed check with a central difference: |dr/ds| should be 1
h = 1e-6
for s in [0.0, 1.0, 5.0]:
print(s, np.linalg.norm((r_of_s(s + h) - r_of_s(s - h)) / (2 * h)))
# 0.0 0.9999999998968969
# 1.0 1.0000000000987614
# 5.0 0.9999999999670178
# When s(t) has no closed-form inverse: tabulate s(t), then interpolate t(s)
def rp(t): # twisted cubic <t, t^2, t^3>: velocity
return np.array([np.ones_like(t), 2 * t, 3 * t**2])
t = np.linspace(0, 1, 2001)
speed = np.linalg.norm(rp(t), axis=0)
s = np.concatenate([[0], np.cumsum(0.5 * (speed[1:] + speed[:-1]) * np.diff(t))]) # cumulative trapezoid
L = s[-1]
s_targets = np.linspace(0, L, 5) # 5 positions equally spaced along the curve
t_at = np.interp(s_targets, s, t) # inverse lookup: swap the roles of s and t
print(round(L, 6), np.round(t_at, 4)) # 1.863023 [0. 0.4137 0.6797 0.8611 1. ]
pts = np.array([t_at, t_at**2, t_at**3]).T
print(np.round(np.linalg.norm(np.diff(pts, axis=0), axis=1), 4))
# [0.4533 0.4631 0.4652 0.4656] chords slightly under L/4 = 0.4658, most where the curve bends most
The equally spaced arc lengths correspond to unequal steps in $t$ (larger near $t = 0$, where the curve moves slowly), which is the whole point of reparametrizing.
Key Takeaways
- $s(t) = \int_a^t \lvert \mathbf{r}'(u) \rvert,du$ measures distance along the curve, and $\frac{ds}{dt} = \lvert \mathbf{r}'(t) \rvert$.
- Reparametrizing by $s$ gives a unit-speed curve: $\left\lvert \frac{d\mathbf{r}}{ds} \right\rvert = 1$.
- Inverting $s(t)$ in closed form is rare. Arc length is mainly used to define geometric quantities, and $\frac{d}{ds} = \frac{1}{\lvert \mathbf{r}' \rvert}\frac{d}{dt}$ converts them to any parameter.
- Numerically, tabulate $s(t)$ and invert by interpolation.
🧪 Practice
- Reparametrize $\mathbf{r}(t) = \langle 3 - 2t,\ 1 + 4t,\ 4t \rangle$ by arc length from the point where $t = 0$. Where is the point $6$ units along the line in the direction of increasing $t$?
- Reparametrize the helix $\langle 2\cos t, 2\sin t, t \rangle$ by arc length measured from $(2, 0, 0)$, and find the point reached after traveling $5\pi$ units.
- Show that $\mathbf{r}(s) = \left\langle \frac{s}{\sqrt{2}},\ \frac{1}{\sqrt{2}}\cos s,\ \frac{1}{\sqrt{2}}\sin s \right\rangle$ is already parametrized by arc length, and explain how you know without computing $s(t)$.
- Interview-style: Your tabulated $s(t)$ for a camera path has 2001 samples, but the camera must be positioned 60 times per second for a 30-second shot. How do you make each lookup cheap, and what accuracy issues arise if the path has a nearly stationary section ($\lvert \mathbf{r}' \rvert \approx 0$)? Hint: the table is sorted, so binary search plus linear interpolation is $O(\log n)$. Near zero speed, $t(s)$ is very steep, so small errors in $s$ cause large jumps in $t$; refine the samples there.
Unit Tangent and Normal Vectors
Theory
To describe how a curve bends, we need reference directions attached to the curve itself, not to the fixed $x, y, z$ axes. Think of a pilot. The useful directions are "forward", "toward the inside of the turn", and "the wing direction", and they rotate with the plane. This topic builds the first two.
Unit tangent vector. At a point where the curve is smooth ($\mathbf{r}' \ne \mathbf{0}$), the direction of motion is
$$\mathbf{T}(t) = \frac{\mathbf{r}'(t)}{\lvert \mathbf{r}'(t) \rvert}.$$
Dividing by the speed removes the timing information, so $\mathbf{T}$ depends only on the curve and its orientation. In terms of arc length, $\mathbf{T} = \frac{d\mathbf{r}}{ds}$.
Principal unit normal vector. Since $\lvert \mathbf{T} \rvert = 1$ for all $t$, the constant-length rule (Section 11.2) gives $\mathbf{T} \cdot \mathbf{T}' = 0$. A unit vector can only change by turning, so its derivative is perpendicular to it. The direction in which $\mathbf{T}$ turns is the direction in which the curve bends. Where $\mathbf{T}' \ne \mathbf{0}$, define
$$\mathbf{N}(t) = \frac{\mathbf{T}'(t)}{\lvert \mathbf{T}'(t) \rvert}.$$
$\mathbf{N}$ is a unit vector, perpendicular to $\mathbf{T}$, pointing toward the concave side of the curve: the inside of the turn. On a straight line, $\mathbf{T}$ never changes, $\mathbf{T}' = \mathbf{0}$, and $\mathbf{N}$ is undefined. A line does not turn, so there is no "inside".
*-------> T the point * moves to the right along
.''' | '''. the top of the circle (clockwise)
.' | '.
: v N : T: unit vector in the direction of motion
: o : N: unit vector toward the inside of the bend,
: center : here straight down toward the center o
'. .'
'.........'
A practical note: computing $\mathbf{T}'$ by differentiating $\frac{\mathbf{r}'}{\lvert \mathbf{r}' \rvert}$ is often messy because of the quotient. Two alternatives:
- For curves with constant speed, $\mathbf{T}'$ is just $\frac{\mathbf{r}''}{\lvert \mathbf{r}' \rvert}$, so $\mathbf{N}$ is the unit vector of $\mathbf{r}''$.
- In general, $\mathbf{N}$ is the unit vector in the direction of the part of $\mathbf{r}''$ perpendicular to $\mathbf{r}'$: $\mathbf{r}'' - \operatorname{proj}_{\mathbf{r}'}\mathbf{r}''$. This follows from the acceleration decomposition in Section 11.4, and it is how $\mathbf{N}$ is computed in numerical code.
Examples
For the helix $\mathbf{r}(t) = \langle a\cos t, a\sin t, ct \rangle$, the speed is $\sqrt{a^2 + c^2}$ and
$$\mathbf{T}(t) = \frac{1}{\sqrt{a^2 + c^2}}\langle -a\sin t,\ a\cos t,\ c \rangle, \qquad \mathbf{T}'(t) = \frac{1}{\sqrt{a^2 + c^2}}\langle -a\cos t,\ -a\sin t,\ 0 \rangle,$$
$$\mathbf{N}(t) = \langle -\cos t,\ -\sin t,\ 0 \rangle.$$
$\mathbf{N}$ is horizontal and points from the point straight toward the $z$-axis, the axis of the cylinder. Climbing a spiral staircase, you are always turning toward the central column.
For the parabola $\mathbf{r}(t) = \langle t, t^2 \rangle$ at $(1, 1)$: $\mathbf{T} = \frac{1}{\sqrt{5}}\langle 1, 2 \rangle$ and $\mathbf{N} = \frac{1}{\sqrt{5}}\langle -2, 1 \rangle$, which points up and to the left, into the cup of the parabola.
import sympy as sp
t = sp.symbols('t', real=True)
a, c = sp.symbols('a c', positive=True)
def unit(v):
return sp.simplify(v / sp.sqrt(sp.simplify(v.dot(v))))
r = sp.Matrix([a * sp.cos(t), a * sp.sin(t), c * t]) # helix
T = unit(r.diff(t)) # T = r' / |r'|
N = unit(T.diff(t)) # N = T' / |T'|
print(T.T) # Matrix([[-a*sin(t)/sqrt(a**2 + c**2), a*cos(t)/sqrt(a**2 + c**2), c/sqrt(a**2 + c**2)]])
print(N.T) # Matrix([[-cos(t), -sin(t), 0]]) horizontal, toward the z-axis
print(sp.simplify(T.dot(N))) # 0
# Plane curve y = x^2 (as a space curve with z = 0) at the point (1, 1)
q = sp.Matrix([t, t**2, 0])
Tq = unit(q.diff(t))
Nq = unit(Tq.diff(t))
print(Tq.subs(t, 1).T, Nq.subs(t, 1).T)
# Matrix([[sqrt(5)/5, 2*sqrt(5)/5, 0]]) Matrix([[-2*sqrt(5)/5, sqrt(5)/5, 0]])
Key Takeaways
- $\mathbf{T} = \frac{\mathbf{r}'}{\lvert \mathbf{r}' \rvert}$ is the direction of motion. It depends on the curve and its orientation, not on the speed.
- $\lvert \mathbf{T} \rvert = 1$ forces $\mathbf{T}' \perp \mathbf{T}$. The unit vector $\mathbf{N} = \frac{\mathbf{T}'}{\lvert \mathbf{T}' \rvert}$ points toward the inside of the bend.
- $\mathbf{N}$ is undefined where the curve is locally straight ($\mathbf{T}' = \mathbf{0}$).
- For constant-speed curves, $\mathbf{N}$ is the direction of $\mathbf{r}''$.
🧪 Practice
- Find $\mathbf{T}(t)$ and $\mathbf{N}(t)$ for $\mathbf{r}(t) = \langle 3\sin t,\ 4t,\ 3\cos t \rangle$.
- Find $\mathbf{T}$ and $\mathbf{N}$ for $\mathbf{r}(t) = \langle t^2, \frac{2}{3}t^3, t \rangle$ at $t = 1$. (Use the projection method: $\mathbf{N}$ is the unit vector along $\mathbf{r}'' - \operatorname{proj}_{\mathbf{r}'}\mathbf{r}''$.)
- For the plane curve $y = \sin x$, find $\mathbf{N}$ at $x = \frac{\pi}{2}$ and at $x = \frac{3\pi}{2}$. What happens at $x = \pi$, and why?
- Interview-style: A car-racing game aligns the car's model with the track using $\mathbf{T}$ and $\mathbf{N}$ computed from a spline. On straight sections, the car flickers between orientations. What is going on, and how would you fix it? Hint: on straight sections, $\mathbf{T}' \approx \mathbf{0}$, so $\mathbf{N}$ is the normalization of a tiny, noisy vector. Use a fixed "up" vector, or a frame carried along without twisting (see the Frenet-Serret Frame topic).
Curvature
Theory
Some roads are gentle sweeps and others are hairpin turns. Curvature puts a number on this: how sharply a curve bends at a point. To make it a property of the curve and not of the driving, measure how fast the direction changes per unit distance traveled, not per unit time.
Definition.
$$\kappa = \left\lvert \frac{d\mathbf{T}}{ds} \right\rvert.$$
A line has $\kappa = 0$. On a circle of radius $a$, going a distance $s$ turns the tangent through an angle $\frac{s}{a}$, so $\kappa = \frac{1}{a}$: small circles bend sharply, and large circles are nearly straight. This motivates the radius of curvature $\rho = \frac{1}{\kappa}$.
Formulas in terms of $t$. Using $\frac{d}{ds} = \frac{1}{\lvert \mathbf{r}' \rvert}\frac{d}{dt}$:
$$\kappa(t) = \frac{\lvert \mathbf{T}'(t) \rvert}{\lvert \mathbf{r}'(t) \rvert} = \frac{\lvert \mathbf{r}'(t) \times \mathbf{r}''(t) \rvert}{\lvert \mathbf{r}'(t) \rvert^3}.$$
The second form avoids computing $\mathbf{T}'$ and is usually the fastest. Proof: write $v = \lvert \mathbf{r}' \rvert$, so $\mathbf{r}' = v\mathbf{T}$. Differentiating, $\mathbf{r}'' = v'\mathbf{T} + v\mathbf{T}'$. Crossing with $\mathbf{r}'$ and using $\mathbf{T} \times \mathbf{T} = \mathbf{0}$ gives $\mathbf{r}' \times \mathbf{r}'' = v^2(\mathbf{T} \times \mathbf{T}')$. Since $\mathbf{T}$ and $\mathbf{T}'$ are perpendicular and $\lvert \mathbf{T} \rvert = 1$, $\lvert \mathbf{T} \times \mathbf{T}' \rvert = \lvert \mathbf{T}' \rvert = v\kappa$. So $\lvert \mathbf{r}' \times \mathbf{r}'' \rvert = v^3\kappa$.
Plane curves. For a graph $y = f(x)$, use $\mathbf{r}(x) = \langle x, f(x), 0 \rangle$. Then $\mathbf{r}' \times \mathbf{r}'' = \langle 0, 0, f''(x) \rangle$, and
$$\kappa(x) = \frac{\lvert f''(x) \rvert}{\left(1 + f'(x)^2\right)^{3/2}}.$$
Curvature is not the same as $\lvert f'' \rvert$, although the two agree where the graph is horizontal ($f' = 0$). On steep parts of the graph, the denominator corrects for the stretching. For a plane parametric curve $\langle x(t), y(t) \rangle$, the formula is $\kappa = \frac{\lvert x'y'' - y'x'' \rvert}{(x'^2 + y'^2)^{3/2}}$.
Osculating circle. At a point where $\kappa > 0$, the circle that best fits the curve has radius $\rho = \frac{1}{\kappa}$, lies in the plane of $\mathbf{T}$ and $\mathbf{N}$, and has its center at $\mathbf{r} + \rho\mathbf{N}$ (the center of curvature). "Osculating" means "kissing". This circle matches the curve's position, tangent direction, and curvature, so it is the curve's second-order approximation, just as the tangent line is the first-order approximation.
y = x^2 osculating circle at the vertex:
radius 1/kappa = 1/2, center (0, 1/2)
\ /
\ /
\ .---. /
\ .' '. /
/ * \ * = center of curvature (0, 1/2)
\ /
'._____.'
^
vertex (0, 0): kappa = 2, the maximum
| Curve | Curvature |
|---|---|
| Line | $0$ |
| Circle of radius $a$ | $\frac{1}{a}$ |
| Helix $\langle a\cos t, a\sin t, ct \rangle$ | $\frac{a}{a^2 + c^2}$ (constant; less than $\frac{1}{a}$, as stretching the coil straightens it) |
| Parabola $y = x^2$ | $\frac{2}{(1 + 4x^2)^{3/2}}$, maximum $2$ at the vertex |
| Twisted cubic $\langle t, t^2, t^3 \rangle$ at $t = 0$ | $2$ |
Applications. For a car on a curve at speed $v$, the sideways acceleration is $\kappa v^2$ (Section 11.4). Road and railway engineers therefore limit curvature, and they avoid sudden jumps in curvature. Optimization of beams, lenses, and computer-vision edge detectors also use it.
Examples
For the helix, $\mathbf{r}' = \langle -a\sin t, a\cos t, c \rangle$ and $\mathbf{r}'' = \langle -a\cos t, -a\sin t, 0 \rangle$. Then $\mathbf{r}' \times \mathbf{r}'' = \langle ac\sin t,\ -ac\cos t,\ a^2 \rangle$, with length $a\sqrt{a^2 + c^2}$, so
$$\kappa = \frac{a\sqrt{a^2 + c^2}}{(a^2 + c^2)^{3/2}} = \frac{a}{a^2 + c^2}.$$
The code below checks these formulas and also estimates curvature from three sample points alone. The circle through three points $P, Q, R$ has curvature $\frac{4 \cdot \text{area}(PQR)}{\lvert PQ \rvert\lvert QR \rvert\lvert PR \rvert}$, and as the points merge, this circle becomes the osculating circle. This is how curvature is estimated from sensor data or polylines.
import numpy as np
import sympy as sp
t = sp.symbols('t', real=True)
def curvature(r):
r1, r2 = r.diff(t), r.diff(t, 2)
return sp.simplify(r1.cross(r2).norm() / r1.norm()**3) # |r' x r''| / |r'|^3
a, c = sp.symbols('a c', positive=True)
print(curvature(sp.Matrix([a * sp.cos(t), a * sp.sin(t), c * t]))) # a/(a**2 + c**2)
print(curvature(sp.Matrix([t, t**2, 0]))) # 2/(4*t**2 + 1)**(3/2)
print(curvature(sp.Matrix([t, t**2, t**3])).subs(t, 0)) # 2
# Numerical curvature from three nearby samples: 1 / radius of the circle through them
def menger_curvature(P, Q, R):
a, b, c = np.linalg.norm(Q - P), np.linalg.norm(R - Q), np.linalg.norm(R - P)
area2 = np.linalg.norm(np.cross(Q - P, R - P)) # twice the triangle's area
return 2 * area2 / (a * b * c) # 4 * area / (a b c)
f = lambda x: np.array([x, x**2, 0.0])
for x0 in [0.0, 0.5, 1.0, 2.0]:
h = 1e-3
k_num = menger_curvature(f(x0 - h), f(x0), f(x0 + h))
k_exact = 2 / (1 + 4 * x0**2)**1.5
print(x0, round(k_num, 6), round(k_exact, 6))
# 0.0 1.999998 2.0
# 0.5 0.707107 0.707107
# 1.0 0.178885 0.178885
# 2.0 0.028534 0.028534 far from the vertex, the parabola is nearly straight
Key Takeaways
- $\kappa = \left\lvert \frac{d\mathbf{T}}{ds} \right\rvert$: the rate at which the direction turns per unit length. A line has $\kappa = 0$; a circle of radius $a$ has $\kappa = \frac{1}{a}$.
- In practice: $\kappa = \frac{\lvert \mathbf{r}' \times \mathbf{r}'' \rvert}{\lvert \mathbf{r}' \rvert^3}$, or $\frac{\lvert f'' \rvert}{(1 + f'^2)^{3/2}}$ for a graph.
- The osculating circle has radius $\frac{1}{\kappa}$ and center $\mathbf{r} + \frac{1}{\kappa}\mathbf{N}$.
- Curvature can be estimated from three nearby points via the circle through them.
🧪 Practice
- Find the curvature of $\mathbf{r}(t) = \langle t, \frac{1}{2}t^2, \frac{1}{3}t^3 \rangle$ at $t = 1$ and of $y = \ln x$ at $x = 1$. Where is the curvature of $y = \ln x$ largest?
- Find the curvature of the ellipse $\langle 3\cos t, 2\sin t \rangle$ and the points where it is largest and smallest. Find the osculating circles at the ends of the major axis.
- Show that the curvature of a polar curve $r = f(\theta)$ is $\kappa = \frac{\lvert r^2 + 2r'^2 - rr'' \rvert}{(r^2 + r'^2)^{3/2}}$, and use it to find the curvature of the cardioid $r = 1 + \cos\theta$ at $\theta = 0$.
- Interview-style: When designing a highway or railway, engineers do not connect a straight section directly to a circular arc. Instead, they insert a transition curve called a clothoid (or Euler spiral). Why, and what property defines a clothoid? Hint: a straight line has $\kappa = 0$ and an arc has $\kappa = \frac{1}{R}$, so a direct joint makes the lateral acceleration $\kappa v^2$ jump instantly, which drivers feel as a jolt. A clothoid's curvature grows linearly with arc length, $\kappa(s) = cs$.
Binormal Vector and Torsion
Theory
Curvature tells you how sharply a curve bends, but in space, there is a second question: is the bending staying within one plane, or is the curve twisting out of it? A circle and a helix both have constant curvature, yet the circle stays flat while the helix climbs. The third frame vector and a second number, the torsion, capture this.
Binormal vector. At points where $\mathbf{N}$ is defined, complete $\mathbf{T}$ and $\mathbf{N}$ to a right-handed frame with
$$\mathbf{B} = \mathbf{T} \times \mathbf{N}.$$
Since $\mathbf{T}$ and $\mathbf{N}$ are orthogonal unit vectors, $\mathbf{B}$ is a unit vector orthogonal to both. In practice, $\mathbf{B}$ is easiest to compute as the unit vector along $\mathbf{r}' \times \mathbf{r}''$, because the curvature proof above showed that $\mathbf{r}' \times \mathbf{r}''$ is a positive multiple of $\mathbf{T} \times \mathbf{T}'$, and hence of $\mathbf{T} \times \mathbf{N}$.
The three vectors define three planes through each point of the curve:
| Plane | Spanned by | Normal vector | Meaning |
|---|---|---|---|
| Osculating plane | $\mathbf{T}, \mathbf{N}$ | $\mathbf{B}$ | The plane the curve is momentarily bending in; contains the osculating circle |
| Normal plane | $\mathbf{N}, \mathbf{B}$ | $\mathbf{T}$ | All directions perpendicular to the motion; a cross-section of a tube around the curve |
| Rectifying plane | $\mathbf{T}, \mathbf{B}$ | $\mathbf{N}$ | The plane the curve "does not bend toward" |
Torsion. A plane curve stays in its osculating plane, so $\mathbf{B}$ is constant. For a twisting curve, the osculating plane rotates, and $\mathbf{B}$ changes. Two facts determine how:
- $\mathbf{B}$ has constant length, so $\frac{d\mathbf{B}}{ds} \perp \mathbf{B}$.
- $\frac{d\mathbf{B}}{ds} = \frac{d\mathbf{T}}{ds} \times \mathbf{N} + \mathbf{T} \times \frac{d\mathbf{N}}{ds} = \kappa\mathbf{N} \times \mathbf{N} + \mathbf{T} \times \frac{d\mathbf{N}}{ds} = \mathbf{T} \times \frac{d\mathbf{N}}{ds}$, which is perpendicular to $\mathbf{T}$.
Perpendicular to both $\mathbf{B}$ and $\mathbf{T}$, it must be a multiple of $\mathbf{N}$. The torsion $\tau$ is defined by
$$\frac{d\mathbf{B}}{ds} = -\tau\mathbf{N}.$$
The minus sign is a convention that makes a right-handed helix (like a standard screw) have positive torsion. Torsion can be positive or negative, unlike curvature: its sign records the handedness of the twist. A curve with $\kappa > 0$ everywhere lies in a plane exactly when $\tau = 0$ everywhere.
Formula in terms of $t$. Differentiating $\mathbf{r}'$, $\mathbf{r}''$, $\mathbf{r}'''$ and expressing them in the $\mathbf{T}, \mathbf{N}, \mathbf{B}$ frame (Section 11.4 does the first two) leads to
$$\tau = \frac{(\mathbf{r}' \times \mathbf{r}'') \cdot \mathbf{r}'''}{\lvert \mathbf{r}' \times \mathbf{r}'' \rvert^2}.$$
The numerator is a scalar triple product, so it measures how far $\mathbf{r}'''$ sticks out of the osculating plane spanned by $\mathbf{r}'$ and $\mathbf{r}''$. If the third derivative stays in that plane, the curve does not twist.
B (binormal: normal to the osculating plane)
^
|
| osculating plane (T, N): the curve bends here
| +-----------------------+
|/ /
*-----------> T / torsion = rate at which this
/| / plane rotates about T as you
/ | / move along the curve
v +--------------------+
N (toward the center of curvature)
Examples
The twisted cubic $\mathbf{r}(t) = \langle t, t^2, t^3 \rangle$ at $t = 0$ has $\mathbf{r}' = \langle 1, 0, 0 \rangle$, $\mathbf{r}'' = \langle 0, 2, 0 \rangle$, and $\mathbf{r}''' = \langle 0, 0, 6 \rangle$. So $\mathbf{r}' \times \mathbf{r}'' = \langle 0, 0, 2 \rangle$ and
$$\mathbf{T} = \mathbf{i}, \quad \mathbf{N} = \mathbf{j}, \quad \mathbf{B} = \mathbf{k}, \quad \kappa = \frac{2}{1} = 2, \quad \tau = \frac{\langle 0, 0, 2 \rangle \cdot \langle 0, 0, 6 \rangle}{4} = 3.$$
At the origin, the osculating plane is the $xy$-plane, $z = 0$. The curve's $z = t^3$ term is what carries it out of that plane.
For the helix $\langle a\cos t, a\sin t, ct \rangle$, the torsion is $\tau = \frac{c}{a^2 + c^2}$, positive for a right-handed helix ($c > 0$) and negative for a left-handed one. As $c \to 0$, the helix flattens into a circle and $\tau \to 0$.
import sympy as sp
t = sp.symbols('t', real=True)
def frenet_data(r):
r1, r2, r3 = r.diff(t), r.diff(t, 2), r.diff(t, 3)
c = r1.cross(r2)
kappa = sp.simplify(c.norm() / r1.norm()**3)
tau = sp.simplify(c.dot(r3) / c.dot(c)) # (r' x r'') . r''' / |r' x r''|^2
B = sp.simplify(c / c.norm()) # B is the unit vector along r' x r''
return kappa, tau, B
a, c = sp.symbols('a c', positive=True)
print(frenet_data(sp.Matrix([a * sp.cos(t), a * sp.sin(t), c * t]))[:2])
# (a/(a**2 + c**2), c/(a**2 + c**2)) right-handed helix: tau > 0
print(frenet_data(sp.Matrix([a * sp.cos(t), a * sp.sin(t), -c * t]))[:2])
# (a/(a**2 + c**2), -c/(a**2 + c**2)) mirror image: same kappa, opposite tau
k, tau, B = frenet_data(sp.Matrix([t, t**2, t**3]))
print(k.subs(t, 0), tau, B.subs(t, 0).T)
# 2 3/(9*t**4 + 9*t**2 + 1) Matrix([[0, 0, 1]]) twisted cubic: tau = 3 at t = 0
print(frenet_data(sp.Matrix([t, t**2, 0]))[1]) # 0 a plane curve does not twist
Key Takeaways
- $\mathbf{B} = \mathbf{T} \times \mathbf{N}$, the unit vector along $\mathbf{r}' \times \mathbf{r}''$, is normal to the osculating plane.
- Torsion is defined by $\frac{d\mathbf{B}}{ds} = -\tau\mathbf{N}$. It measures how fast the osculating plane rotates, and its sign gives the handedness.
- $\tau = \frac{(\mathbf{r}' \times \mathbf{r}'') \cdot \mathbf{r}'''}{\lvert \mathbf{r}' \times \mathbf{r}'' \rvert^2}$.
- A curve with nonzero curvature is planar exactly when $\tau \equiv 0$. Mirror images have the same $\kappa$ and opposite $\tau$.
🧪 Practice
- For $\mathbf{r}(t) = \langle \cos t, \sin t, t \rangle$ at $t = \frac{\pi}{2}$, find $\mathbf{T}$, $\mathbf{N}$, $\mathbf{B}$, and the equations of the osculating and normal planes.
- Compute the torsion of $\mathbf{r}(t) = \langle e^t, e^{-t}, \sqrt{2},t \rangle$ as a function of $t$.
- Show that $\mathbf{r}(t) = \langle 1 + t,\ 2t - t^2,\ 3 + t + t^2 \rangle$ has zero torsion, and find the plane that contains it. (Hint: $\mathbf{B}$ is constant, and it is the plane's normal.)
- Interview-style: You compute $\tau$ numerically from noisy samples of a 3D curve using finite differences for $\mathbf{r}'$, $\mathbf{r}''$, $\mathbf{r}'''$, and the results are wildly unstable, especially on nearly straight sections. Why, and what would you do? Hint: third differences amplify noise like $\frac{\epsilon}{h^3}$, and the denominator $\lvert \mathbf{r}' \times \mathbf{r}'' \rvert^2$ goes to $0$ where $\kappa \to 0$. Smooth first (fit a spline or local polynomial and differentiate that), and do not report $\tau$ where $\kappa$ is below a threshold.
Frenet-Serret Frame
Theory
The vectors $\mathbf{T}, \mathbf{N}, \mathbf{B}$ form a right-handed orthonormal frame that travels with the point along the curve, like the axes of a roller coaster car: forward, toward the inside of the turn, and the axis around which the car is turning. The Frenet-Serret formulas describe how this frame rotates as it moves, and they show that curvature and torsion are all there is to know about the shape of a curve.
The formulas. With respect to arc length,
$$\frac{d\mathbf{T}}{ds} = \kappa\mathbf{N}, \qquad \frac{d\mathbf{N}}{ds} = -\kappa\mathbf{T} + \tau\mathbf{B}, \qquad \frac{d\mathbf{B}}{ds} = -\tau\mathbf{N}.$$
The first equation is the definition of $\kappa$ and $\mathbf{N}$, and the third is the definition of $\tau$. For the middle one, differentiate $\mathbf{N} = \mathbf{B} \times \mathbf{T}$:
$$\frac{d\mathbf{N}}{ds} = \frac{d\mathbf{B}}{ds} \times \mathbf{T} + \mathbf{B} \times \frac{d\mathbf{T}}{ds} = -\tau,\mathbf{N} \times \mathbf{T} + \kappa,\mathbf{B} \times \mathbf{N} = \tau\mathbf{B} - \kappa\mathbf{T}.$$
In matrix form, the coefficient matrix is skew-symmetric:
$$\frac{d}{ds}\begin{pmatrix} \mathbf{T} \ \mathbf{N} \ \mathbf{B} \end{pmatrix} = \begin{pmatrix} 0 & \kappa & 0 \ -\kappa & 0 & \tau \ 0 & -\tau & 0 \end{pmatrix}\begin{pmatrix} \mathbf{T} \ \mathbf{N} \ \mathbf{B} \end{pmatrix}.$$
Skew-symmetry is exactly the condition that keeps a frame orthonormal as it evolves (differentiate $\mathbf{T} \cdot \mathbf{N} = 0$ and the other dot products to see it). In other words, the frame is never stretched or sheared, only rotated. The rotation is described by the Darboux vector $\boldsymbol{\omega} = \tau\mathbf{T} + \kappa\mathbf{B}$: each frame vector satisfies $\frac{d\mathbf{X}}{ds} = \boldsymbol{\omega} \times \mathbf{X}$. Curvature is the rate of rotation about $\mathbf{B}$ (turning, or yaw), and torsion is the rate about $\mathbf{T}$ (twisting, or roll). There is no rotation about $\mathbf{N}$; that is what makes this frame special.
With an arbitrary parameter $t$, multiply each right side by the speed $v = \lvert \mathbf{r}' \rvert$, since $\frac{d}{dt} = v\frac{d}{ds}$. For example, $\mathbf{T}'(t) = \kappa v\mathbf{N}$.
Fundamental theorem of space curves. Given continuous functions $\kappa(s) > 0$ and $\tau(s)$, there is a curve, parametrized by arc length, with exactly this curvature and torsion, and it is unique up to a rigid motion (a rotation plus a translation). Why: the Frenet-Serret formulas are a linear system of differential equations for $\mathbf{T}, \mathbf{N}, \mathbf{B}$. With an initial position and frame, it has a unique solution (Section 8.1), and integrating $\frac{d\mathbf{r}}{ds} = \mathbf{T}$ gives the curve. A different initial position and frame just moves the curve rigidly. So $\kappa$ and $\tau$ are a complete, coordinate-free description of shape.
Consequences:
| $\kappa(s)$ and $\tau(s)$ | Curve |
|---|---|
| $\kappa = 0$ | Straight line |
| $\kappa$ constant $> 0$, $\tau = 0$ | Circle of radius $\frac{1}{\kappa}$ |
| $\kappa$, $\tau$ both constant and nonzero | Circular helix with radius $\frac{\kappa}{\kappa^2 + \tau^2}$ and pitch $\frac{2\pi\tau}{\kappa^2 + \tau^2}$ |
| $\tau = 0$ | Plane curve |
| $\kappa(s) = cs$, $\tau = 0$ | Clothoid (Euler spiral) |
Limitations in practice. The Frenet frame needs $\kappa > 0$. Where a curve is straight or has an inflection point, $\mathbf{N}$ is undefined, and it can flip direction suddenly as the curve passes from bending one way to bending the other. Applications that need a smooth frame everywhere (camera rails, extruding a tube mesh along a path, robot tool paths) often use a rotation-minimizing or parallel-transport frame instead: it keeps $\mathbf{T}$ but rotates the other two vectors as little as possible, avoiding the roll that torsion introduces.
Examples
The fundamental theorem can be tested numerically. The code below takes constant $\kappa = 2$, $\tau = 1$, integrates the Frenet-Serret system from a starting frame, and checks that (a) the frame stays orthonormal, and (b) the curve is a helix with the predicted radius $\frac{\kappa}{\kappa^2 + \tau^2} = 0.4$, wound around an axis along the Darboux vector.
import numpy as np
from scipy.integrate import solve_ivp
kappa, tau = 2.0, 1.0 # constant curvature and torsion
def frenet_rhs(s, y):
r, T, N, B = y.reshape(4, 3)
return np.concatenate([T, # dr/ds = T
kappa * N, # dT/ds = kappa N
-kappa * T + tau * B, # dN/ds = -kappa T + tau B
-tau * N]) # dB/ds = -tau N
y0 = np.concatenate([np.zeros(3), np.eye(3).ravel()]) # start at the origin with T, N, B = i, j, k
sol = solve_ivp(frenet_rhs, (0, 20), y0, rtol=1e-10, atol=1e-12, dense_output=True)
s = np.linspace(0, 20, 400)
r, T, N, B = sol.sol(s).reshape(4, 3, -1) # each has shape (3, 400)
# (a) The frame stays orthonormal: all pairwise dot products match the identity matrix
F = np.stack([T, N, B])
G = np.einsum('ais,bis->sab', F, F) # Gram matrix at each s
print(np.abs(G - np.eye(3)).max()) # 1.1557962364960872e-09 (integration error only)
# (b) Distance from the predicted axis is constant: a helix of radius kappa / (kappa^2 + tau^2)
R_pred = kappa / (kappa**2 + tau**2)
axis = tau * np.array([1, 0, 0]) + kappa * np.array([0, 0, 1]) # Darboux vector tau T(0) + kappa B(0)
axis /= np.linalg.norm(axis)
center0 = N[:, 0] * R_pred # the axis passes through r(0) + R N(0)
d = r.T - center0
radial = d - np.outer(d @ axis, axis) # remove the component along the axis
dist = np.linalg.norm(radial, axis=1)
print(R_pred, round(dist.min(), 8), round(dist.max(), 8)) # 0.4 0.4 0.4
print(2 * np.pi * tau / (kappa**2 + tau**2)) # 1.2566370614359172 predicted pitch
Key Takeaways
- $\mathbf{T}' = \kappa\mathbf{N}$, $\mathbf{N}' = -\kappa\mathbf{T} + \tau\mathbf{B}$, $\mathbf{B}' = -\tau\mathbf{N}$ (derivatives with respect to $s$; multiply by the speed for any other parameter).
- The coefficient matrix is skew-symmetric, so the frame only rotates. The rotation vector is $\boldsymbol{\omega} = \tau\mathbf{T} + \kappa\mathbf{B}$.
- Curvature and torsion as functions of arc length determine a curve up to a rigid motion. Constant $\kappa$ and $\tau$ give a helix.
- The frame breaks down where $\kappa = 0$. Graphics and robotics often use rotation-minimizing frames instead.
🧪 Practice
- Verify all three Frenet-Serret formulas for the helix $\langle \cos t, \sin t, t \rangle$, remembering to divide $t$-derivatives by the speed $\sqrt{2}$.
- Show that a helix with radius $a$ and pitch $2\pi c$ has $\kappa = \frac{a}{a^2 + c^2}$ and $\tau = \frac{c}{a^2 + c^2}$, then invert these relations to recover $a = \frac{\kappa}{\kappa^2 + \tau^2}$ and $c = \frac{\tau}{\kappa^2 + \tau^2}$.
- Prove that if all the normal lines of a curve (the lines through $\mathbf{r}(s)$ in the direction of $\mathbf{N}(s)$) pass through one fixed point, then the curve lies on a circle. (Hint: write the fixed point as $\mathbf{r} + \lambda\mathbf{N}$ and differentiate.)
- Interview-style: You are extruding a circular tube mesh along a 3D spline by placing a ring of vertices in the plane of $\mathbf{N}$ and $\mathbf{B}$ at each sample. The tube looks twisted and, at some samples, collapses or flips. What causes this, and what frame would you use instead? Hint: torsion rotates $\mathbf{N}$ and $\mathbf{B}$ around $\mathbf{T}$, and $\mathbf{N}$ is undefined at inflection points. A parallel-transport frame (for example, the double-reflection method) carries the previous ring's orientation forward, rotating it only enough to stay perpendicular to the new $\mathbf{T}$.
<a id="114-motion-in-space"></a>
11.4 Motion in Space
When $\mathbf{r}(t)$ is the position of a moving object, its derivatives are velocity and acceleration, and Newton's second law connects acceleration to force. This subchapter splits acceleration into a part that changes speed and a part that changes direction, solves projectile problems, and derives Kepler's first law from the inverse-square law with vector calculus.
Velocity and Acceleration Vectors
Theory
Section 4.6 studied motion along a line, where velocity is a signed number. In space, the direction of motion can change too, so velocity and acceleration become vectors. If $\mathbf{r}(t)$ is the position at time $t$:
| Quantity | Definition | Meaning |
|---|---|---|
| Velocity | $\mathbf{v}(t) = \mathbf{r}'(t)$ | Direction of motion and rate of travel; tangent to the path |
| Speed | $v(t) = \lvert \mathbf{v}(t) \rvert = \frac{ds}{dt}$ | How fast, regardless of direction; what a speedometer shows |
| Acceleration | $\mathbf{a}(t) = \mathbf{v}'(t) = \mathbf{r}''(t)$ | Rate of change of velocity, in both size and direction |
The key conceptual point is in the last row. Acceleration is not just "speeding up". A car driving around a bend at a constant 60 km/h is accelerating, because its velocity vector is turning. You feel this acceleration as a sideways push. In everyday language, the gas pedal, the brake, and the steering wheel are all accelerators.
Newton's second law. For an object of constant mass $m$ acted on by a total force $\mathbf{F}$,
$$\mathbf{F}(t) = m,\mathbf{a}(t).$$
This is three equations at once, one per component. It works in both directions:
- Given the motion, find the force. Differentiate $\mathbf{r}$ twice and multiply by $m$.
- Given the force, find the motion. Integrate $\mathbf{a} = \frac{\mathbf{F}}{m}$ twice, using the initial velocity $\mathbf{v}(0) = \mathbf{v}_0$ and initial position $\mathbf{r}(0) = \mathbf{r}_0$ to fix the two constant vectors (Section 11.2):
$$\mathbf{v}(t) = \mathbf{v}_0 + \int_0^t \mathbf{a}(u),du, \qquad \mathbf{r}(t) = \mathbf{r}_0 + \int_0^t \mathbf{v}(u),du.$$
When the force depends on position or velocity (springs, gravity toward a center, air resistance), the second direction becomes a system of differential equations (Chapter 8), usually solved numerically.
From data. When positions are measured rather than given by a formula, velocity and acceleration are estimated with finite differences: $\mathbf{v}(t) \approx \frac{\mathbf{r}(t + h) - \mathbf{r}(t - h)}{2h}$ and $\mathbf{a}(t) \approx \frac{\mathbf{r}(t + h) - 2\mathbf{r}(t) + \mathbf{r}(t - h)}{h^2}$, applied to each component.
Examples
For $\mathbf{r}(t) = \langle t^2,\ e^t,\ te^t \rangle$:
$$\mathbf{v}(t) = \langle 2t,\ e^t,\ (1 + t)e^t \rangle, \qquad \mathbf{a}(t) = \langle 2,\ e^t,\ (2 + t)e^t \rangle.$$
At $t = 0$, the object is at $(0, 1, 0)$ with velocity $\langle 0, 1, 1 \rangle$, speed $\sqrt{2}$, and acceleration $\langle 2, 1, 2 \rangle$.
Going the other way: a ball thrown from $\mathbf{r}_0 = \langle 0, 0, 5 \rangle$ with $\mathbf{v}_0 = \langle 1, 2, 3 \rangle$ has acceleration $\langle 0, 0, -10 \rangle$ (taking $g \approx 10$ m/s²). Integrating once gives $\mathbf{v}(t) = \langle 1, 2, 3 - 10t \rangle$, and again gives $\mathbf{r}(t) = \langle t,\ 2t,\ 5 + 3t - 5t^2 \rangle$.
import numpy as np
import sympy as sp
t = sp.symbols('t', real=True)
r = sp.Matrix([t**2, sp.exp(t), t * sp.exp(t)])
v, acc = r.diff(t), r.diff(t, 2)
print(v.T, acc.T)
# Matrix([[2*t, exp(t), t*exp(t) + exp(t)]]) Matrix([[2, exp(t), (t + 2)*exp(t)]])
print(v.subs(t, 0).T, sp.sqrt(v.subs(t, 0).dot(v.subs(t, 0))), acc.subs(t, 0).T)
# Matrix([[0, 1, 1]]) sqrt(2) Matrix([[2, 1, 2]])
# From acceleration back to position: a = <0, 0, -10>, v(0) = <1, 2, 3>, r(0) = <0, 0, 5>
a = sp.Matrix([0, 0, -10])
vel = sp.Matrix([1, 2, 3]) + a.integrate((t, 0, t))
pos = sp.Matrix([0, 0, 5]) + vel.integrate((t, 0, t))
print(vel.T, pos.T)
# Matrix([[1, 2, 3 - 10*t]]) Matrix([[t, 2*t, -5*t**2 + 3*t + 5]])
# Finite differences recover v and a from sampled positions (here at t = 0.5)
f = sp.lambdify(t, r, 'numpy')
h, t0 = 1e-4, 0.5
r0, rp, rm = (np.array(f(x), float).ravel() for x in (t0, t0 + h, t0 - h))
print(np.round((rp - rm) / (2 * h), 6), np.round(np.array(v.subs(t, t0).evalf(), float).ravel(), 6))
# [1. 1.648721 2.473082] [1. 1.648721 2.473082] central difference vs exact v
print(np.round((rp - 2 * r0 + rm) / h**2, 6), np.round(np.array(acc.subs(t, t0).evalf(), float).ravel(), 6))
# [2. 1.648721 4.121803] [2. 1.648721 4.121803] second difference vs exact a
Key Takeaways
- $\mathbf{v} = \mathbf{r}'$ is tangent to the path, the speed is $\lvert \mathbf{v} \rvert$, and $\mathbf{a} = \mathbf{r}''$.
- Acceleration measures any change in velocity: speeding up, slowing down, or turning.
- $\mathbf{F} = m\mathbf{a}$ is one equation per component. Integrate twice with $\mathbf{r}_0$ and $\mathbf{v}_0$ to recover the motion from the force.
- From sampled data, estimate $\mathbf{v}$ and $\mathbf{a}$ with central differences, component by component.
🧪 Practice
- Find the velocity, speed, and acceleration of $\mathbf{r}(t) = \langle 3\cos t,\ 3\sin t,\ t^2 \rangle$ at $t = \pi$. When is the speed smallest?
- A 2 kg object starts at the origin with velocity $\langle 1, -1, 3 \rangle$ and is acted on by the force $\mathbf{F}(t) = \langle 6t,\ 4,\ -2e^{-t} \rangle$ newtons. Find its position at time $t$.
- Show that if the speed of an object is constant, then its acceleration is always perpendicular to its velocity. Is the converse true?
- Interview-style: You estimate a drone's acceleration from GPS positions with the second difference $\frac{\mathbf{r}(t + h) - 2\mathbf{r}(t) + \mathbf{r}(t - h)}{h^2}$. Increasing the sampling rate (smaller $h$) makes the estimate worse, not better. Why? Hint: if each position has a measurement error of size $\epsilon$, the error in the estimate is about $\frac{4\epsilon}{h^2}$, while the truncation error shrinks only like $h^2$. Balance the two, or smooth the positions (or fuse in an accelerometer) before differentiating.
Tangential and Normal Components of Acceleration
Theory
When you ride in a car, you feel two different kinds of acceleration: the push into your seat when the driver speeds up (or the lurch forward when braking), and the sideways push in a turn. These are the components of $\mathbf{a}$ along $\mathbf{T}$ and along $\mathbf{N}$, and separating them shows exactly how much of the acceleration changes speed and how much changes direction.
The decomposition. Write the velocity as speed times direction, $\mathbf{v} = v\mathbf{T}$. Differentiating with the product rule,
$$\mathbf{a} = v'\mathbf{T} + v\mathbf{T}'.$$
From the Frenet-Serret frame (Section 11.3), $\mathbf{T}' = \frac{d\mathbf{T}}{ds}\frac{ds}{dt} = \kappa v\mathbf{N}$. So
$$\boxed{\mathbf{a} = \underbrace{v'}{a_T},\mathbf{T} + \underbrace{\kappa v^2}{a_N},\mathbf{N}}$$
| Component | Formula | Effect | Felt as |
|---|---|---|---|
| Tangential $a_T$ | $\frac{dv}{dt}$ | Changes the speed. Positive when speeding up, negative when slowing down | Pushed into the seat, or thrown forward |
| Normal $a_N$ | $\kappa v^2 = \frac{v^2}{\rho}$ | Changes the direction. Always $\ge 0$, pointing into the turn | Pushed sideways, toward the outside of the turn |
Three consequences:
- There is no $\mathbf{B}$ component. Acceleration always lies in the osculating plane, whatever the motion.
- $a_N$ grows with the square of speed. Taking a curve twice as fast needs four times the sideways force (friction from the tires). This is why speed limits drop sharply on tight curves.
- $a_N$ grows with curvature. At the same speed, a tighter curve (larger $\kappa$, smaller radius $\rho$) needs more sideways force.
Computing the components without finding $\mathbf{T}$, $\mathbf{N}$, or $\kappa$ first. Dot and cross the decomposition with $\mathbf{v} = v\mathbf{T}$:
- $\mathbf{v} \cdot \mathbf{a} = v,a_T$ (since $\mathbf{T} \cdot \mathbf{T} = 1$ and $\mathbf{T} \cdot \mathbf{N} = 0$), so $a_T = \dfrac{\mathbf{v} \cdot \mathbf{a}}{\lvert \mathbf{v} \rvert}$, the scalar projection of $\mathbf{a}$ onto $\mathbf{v}$ (Section 11.1).
- $\mathbf{v} \times \mathbf{a} = v,a_N,(\mathbf{T} \times \mathbf{N}) = v,a_N\mathbf{B}$, so $a_N = \dfrac{\lvert \mathbf{v} \times \mathbf{a} \rvert}{\lvert \mathbf{v} \rvert}$.
- Since $\mathbf{T}$ and $\mathbf{N}$ are orthonormal, $\lvert \mathbf{a} \rvert^2 = a_T^2 + a_N^2$, which gives $a_N = \sqrt{\lvert \mathbf{a} \rvert^2 - a_T^2}$ as a shortcut.
The second formula, combined with $a_N = \kappa v^2$, gives $\kappa = \frac{\lvert \mathbf{v} \times \mathbf{a} \rvert}{v^3}$, which is the curvature formula of Section 11.3 in physical language.
a (total acceleration)
^
/|
/ |
/ | a_N N (turning: points into the curve)
/ |
*----+------> a_T T (speeding up: along the motion)
.' path
.' a = a_T T + a_N N, |a|^2 = a_T^2 + a_N^2
.'
Examples
For $\mathbf{r}(t) = \langle t, t^2, t^3 \rangle$ at $t = 1$: $\mathbf{v} = \langle 1, 2, 3 \rangle$ and $\mathbf{a} = \langle 0, 2, 6 \rangle$. Then
$$a_T = \frac{\mathbf{v} \cdot \mathbf{a}}{\lvert \mathbf{v} \rvert} = \frac{0 + 4 + 18}{\sqrt{14}} = \frac{22}{\sqrt{14}} \approx 5.880, \qquad \mathbf{v} \times \mathbf{a} = \langle 6, -6, 2 \rangle, \quad a_N = \frac{\sqrt{76}}{\sqrt{14}} \approx 2.330.$$
The object is speeding up ($a_T > 0$), and most of its acceleration goes into that rather than into turning. The curvature there is $\kappa = \frac{a_N}{v^2} = \frac{2.330}{14} \approx 0.166$.
import numpy as np
def accel_components(v, a):
speed = np.linalg.norm(v)
aT = v @ a / speed # signed: > 0 speeding up, < 0 slowing down
aN = np.linalg.norm(np.cross(v, a)) / speed # always >= 0
kappa = aN / speed**2 # from a_N = kappa v^2
return aT, aN, kappa
# r(t) = <t, t^2, t^3> at t = 1
v, a = np.array([1.0, 2.0, 3.0]), np.array([0.0, 2.0, 6.0])
aT, aN, kappa = accel_components(v, a)
print(aT, aN, kappa) # 5.879747322073337 2.32992949004287 0.16642353500306215
print(np.hypot(aT, aN), np.linalg.norm(a)) # 6.324555320336759 6.324555320336759 Pythagoras holds
# Highway curve design: largest speed keeping lateral acceleration below 0.3 g
g = 9.81
for radius in [50, 200, 800]: # curve radius in meters
v_max = np.sqrt(0.3 * g * radius) # a_N = v^2 / R <= 0.3 g
print(radius, round(v_max, 1), round(v_max * 3.6)) # m/s and km/h
# 50 12.1 44
# 200 24.3 87
# 800 48.5 175 quadrupling the radius only doubles the safe speed
Key Takeaways
- $\mathbf{a} = a_T\mathbf{T} + a_N\mathbf{N}$, with $a_T = v'$ (changes speed) and $a_N = \kappa v^2$ (changes direction). There is never a $\mathbf{B}$ component.
- $a_T = \frac{\mathbf{v} \cdot \mathbf{a}}{\lvert \mathbf{v} \rvert}$, $a_N = \frac{\lvert \mathbf{v} \times \mathbf{a} \rvert}{\lvert \mathbf{v} \rvert} = \sqrt{\lvert \mathbf{a} \rvert^2 - a_T^2}$.
- Sideways acceleration grows with the square of the speed and linearly with the curvature: $v_{\max} = \sqrt{a_{\max}/\kappa}$.
🧪 Practice
- Find $a_T$ and $a_N$ for $\mathbf{r}(t) = \langle e^t,\ \sqrt{2},t,\ e^{-t} \rangle$ as functions of $t$.
- Find $a_T$ and $a_N$ for the helix $\langle \cos t, \sin t, t \rangle$, and explain why $a_T = 0$ without computing it.
- A car travels around a circular track of radius $100$ m, speeding up at a constant rate of $2$ m/s². At the moment its speed is $20$ m/s, find the magnitude of its acceleration and the angle between the acceleration and the velocity.
- Interview-style: A self-driving car's planner must slow down for curves so that passengers never feel more than $2$ m/s² of sideways acceleration, and it must also brake no harder than $3$ m/s². Given a planned path with known curvature $\kappa(s)$, how would you compute a speed profile $v(s)$? Hint: the curvature limit gives a ceiling $v(s) \le \sqrt{2/\kappa(s)}$. The braking limit means $v$ cannot drop faster than $\frac{d(v^2)}{ds} \ge -6$, so sweep backward along the path from each tight curve to start slowing in time. (Some planners also cap the combined $\sqrt{a_T^2 + a_N^2}$.)
Projectile Motion
Theory
A thrown ball, a fired shell, a jet of water from a fountain: once launched, each moves under gravity alone (if we ignore air). Projectile motion is the cleanest example of finding a trajectory from a force by integration, and it reveals some surprisingly simple geometry.
Setup. Put the launch point at the origin, with the $z$-axis pointing up. Near the Earth's surface, gravity gives a constant acceleration $\mathbf{a} = -g\mathbf{k}$, with $g \approx 9.81$ m/s². With initial velocity $\mathbf{v}_0$, integrating twice gives
$$\mathbf{v}(t) = \mathbf{v}_0 - gt,\mathbf{k}, \qquad \mathbf{r}(t) = \mathbf{v}_0 t - \frac{1}{2}gt^2,\mathbf{k}.$$
The trajectory always lies in the vertical plane containing $\mathbf{v}_0$, since $\mathbf{r}(t)$ is a combination of $\mathbf{v}_0$ and $\mathbf{k}$. So the problem is really two-dimensional. Use a horizontal axis $x$ in the launch direction. If the launch speed is $v_0$ and the angle of elevation is $\alpha$, then $\mathbf{v}_0 = \langle v_0\cos\alpha,\ v_0\sin\alpha \rangle$ and
$$x(t) = (v_0\cos\alpha),t, \qquad z(t) = (v_0\sin\alpha),t - \frac{1}{2}gt^2.$$
The horizontal velocity never changes (no horizontal force), and the vertical motion is the falling-body motion of Section 4.6.
Shape of the path. Eliminating $t = \frac{x}{v_0\cos\alpha}$ gives
$$z = (\tan\alpha),x - \frac{g}{2v_0^2\cos^2\alpha},x^2,$$
a downward parabola, as Galileo found.
Key quantities (launch and landing at the same height):
| Quantity | Formula | How to get it |
|---|---|---|
| Time of flight | $T = \frac{2v_0\sin\alpha}{g}$ | Solve $z(t) = 0$ for $t > 0$ |
| Range | $R = \frac{v_0^2\sin 2\alpha}{g}$ | $x(T)$, using $2\sin\alpha\cos\alpha = \sin 2\alpha$ |
| Maximum height | $H = \frac{(v_0\sin\alpha)^2}{2g}$ | $z$ at $t = \frac{T}{2}$, where $z'(t) = 0$ |
| Speed at landing | $v_0$ | By symmetry (and conservation of energy) |
The range formula has two famous consequences. The maximum range, $\frac{v_0^2}{g}$, is reached at $\alpha = 45°$, where $\sin 2\alpha = 1$. And complementary angles such as $30°$ and $60°$ give the same range, because $\sin 2\alpha = \sin(180° - 2\alpha)$: the higher shot stays in the air longer but moves forward more slowly.
z v0 = 50 m/s, no air resistance
|
| . ' ' . alpha = 60 deg: H = 95.6 m, T = 8.83 s
| .' '.
| .' . ' . '. alpha = 45 deg: H = 63.7 m, R = 254.8 m (maximum)
| / .' '. \
| / .' .' ' '. '. \ alpha = 30 deg: H = 31.9 m, T = 5.10 s
| / .' .' '. '. \
| /.' .' '. '.\
| /'.' '. '\
+-------*--------------------------*---*------- x
0 220.7 254.8
30 and 60 degrees land at the same point
Launching from a height $h_0$: solve $h_0 + (v_0\sin\alpha)t - \frac{1}{2}gt^2 = 0$ with the quadratic formula and take the positive root. The optimal angle is then a little below $45°$.
Air resistance. Real projectiles are slowed by drag. The simplest model, linear drag, adds a force opposite to the velocity: $\mathbf{a} = -g\mathbf{k} - b\mathbf{v}$, with drag coefficient $b > 0$ (per second). This is a linear differential equation for $\mathbf{v}$ (Section 8.2), with solution
$$\mathbf{v}(t) = \mathbf{v}\infty + (\mathbf{v}0 - \mathbf{v}\infty)e^{-bt}, \qquad \mathbf{v}\infty = -\frac{g}{b}\mathbf{k}.$$
The velocity relaxes exponentially toward the terminal velocity $\mathbf{v}_\infty$: the horizontal motion dies out, and the projectile ends up falling straight down at speed $\frac{g}{b}$. The path is no longer a parabola, it comes down more steeply than it went up, and the best launch angle drops below $45°$. Fast projectiles experience drag closer to $\lvert \mathbf{v} \rvert\mathbf{v}$ (quadratic drag), which has no closed-form solution and is integrated numerically.
Examples
A ball is kicked at $50$ m/s at $45°$. With $g = 9.81$: $T = \frac{2(50)\sin 45°}{9.81} \approx 7.21$ s, $R = \frac{2500}{9.81} \approx 254.8$ m, and $H = \frac{(50\sin 45°)^2}{19.62} \approx 63.7$ m. The code checks these, then adds linear drag with $b = 0.1$ s$^{-1}$ by integrating numerically and stopping when the projectile lands.
import numpy as np
from scipy.integrate import solve_ivp
g = 9.81
def ideal(v0, alpha_deg, h0=0.0):
"""No air resistance. Returns (flight time, range, max height)."""
a = np.radians(alpha_deg)
vx, vz = v0 * np.cos(a), v0 * np.sin(a)
T = (vz + np.sqrt(vz**2 + 2 * g * h0)) / g # positive root of h0 + vz t - g t^2 / 2 = 0
return T, vx * T, h0 + vz**2 / (2 * g)
print(np.round(ideal(50, 45), 2)) # [ 7.21 254.84 63.71]
print(np.round(ideal(50, 30), 2), np.round(ideal(50, 60), 2))
# [ 5.1 220.7 31.86] [ 8.83 220.7 95.57] complementary angles: same range
print(np.round(ideal(50, 45, h0=20), 2)) # [ 7.74 273.48 83.71] launched from 20 m up
def range_with_drag(v0, alpha_deg, b):
"""Linear drag: a = -g k - b v. Integrates until the height returns to 0."""
a = np.radians(alpha_deg)
def rhs(t, y):
x, z, vx, vz = y
return [vx, vz, -b * vx, -g - b * vz] # r' = v, v' = a
landed = lambda t, y: y[1] # event function: height z
landed.terminal, landed.direction = True, -1 # stop when z crosses 0 going down
sol = solve_ivp(rhs, (0, 100), [0, 0, v0 * np.cos(a), v0 * np.sin(a)], events=landed, rtol=1e-9)
return sol.y_events[0][0][0] # x at the landing event
print(round(range_with_drag(50, 45, 0.0), 2), round(range_with_drag(50, 45, 0.1), 2))
# 254.84 169.12 drag costs a third of the range
angles = np.arange(30, 46)
ranges = [range_with_drag(50, A, 0.1) for A in angles]
print(angles[int(np.argmax(ranges))], round(max(ranges), 2))
# 39 172.06 with drag, the best angle drops below 45 degrees
Key Takeaways
- With $\mathbf{a} = -g\mathbf{k}$: $\mathbf{r}(t) = \mathbf{r}_0 + \mathbf{v}_0 t - \frac{1}{2}gt^2\mathbf{k}$. The path is a parabola in the vertical plane of $\mathbf{v}_0$.
- Horizontal velocity is constant, and the vertical motion is free fall.
- On level ground: $T = \frac{2v_0\sin\alpha}{g}$, $R = \frac{v_0^2\sin 2\alpha}{g}$ (maximum at $45°$; complementary angles tie), $H = \frac{(v_0\sin\alpha)^2}{2g}$.
- Linear drag gives exponential approach to a terminal velocity and lowers the optimal angle. Realistic drag needs numerical integration.
🧪 Practice
- A ball is thrown at $20$ m/s at $35°$ above the horizontal from a height of $1.8$ m. Find when and where it lands, and its speed at landing.
- A basketball player releases the ball from $2$ m above the floor, $6$ m from the hoop (which is $3.05$ m high), at an angle of $50°$. What launch speed puts the ball through the hoop?
- Show that a projectile launched with speed $v_0$ can reach a point $(x, z)$ only if $z \le \frac{v_0^2}{2g} - \frac{gx^2}{2v_0^2}$. (This parabola, the envelope of all trajectories, is called the safety parabola. Hint: treat $\tan\alpha$ as the unknown in the trajectory equation and require real roots.)
- Interview-style: In a game, an enemy turret must hit a target moving at constant velocity $\mathbf{u}$, firing projectiles with a fixed speed $v_0$ under gravity. How would you compute the firing direction? Hint: work in the frame of the target, or require $\mathbf{p} + \mathbf{u}t = \mathbf{v}t - \frac{1}{2}gt^2\mathbf{k}$ with $\lvert \mathbf{v} \rvert = v_0$. Taking the length of $\mathbf{v}t$ gives a quartic in $t$; solve it numerically and pick the smallest positive root (or the larger one for a lob).
Circular and Planetary Motion
Theory
Rotation is everywhere: wheels, centrifuges, satellites, planets. This topic first analyzes the simplest case, uniform circular motion, and then uses the full vector toolkit of this chapter to derive Kepler's laws from Newton's law of gravitation. Section 10.3 stated those laws and derived the first one in polar coordinates; here the vector proof shows why the orbit is a conic, using only cross products and the product rule.
Uniform circular motion. An object moving counterclockwise on a circle of radius $R$ with constant angular speed $\omega$ (radians per second) has
$$\mathbf{r}(t) = R\langle \cos\omega t,\ \sin\omega t \rangle, \qquad \mathbf{v}(t) = R\omega\langle -\sin\omega t,\ \cos\omega t \rangle, \qquad \mathbf{a}(t) = -\omega^2,\mathbf{r}(t).$$
- The velocity is perpendicular to the radius ($\mathbf{r} \cdot \mathbf{v} = 0$), as it must be for any motion at constant distance from a center (Section 11.2).
- The speed is $v = R\omega$, and the period is $T = \frac{2\pi}{\omega}$.
- The acceleration points toward the center (centripetal) with magnitude $R\omega^2 = \frac{v^2}{R}$. This is consistent with the previous topics: the speed is constant, so $a_T = 0$, and $a_N = \kappa v^2 = \frac{v^2}{R}$.
By Newton's second law, keeping an object of mass $m$ on the circle takes a net force $\frac{mv^2}{R}$ directed toward the center. It is supplied by string tension for a whirled stone, by friction for a car in a turn, and by gravity for a satellite. The outward "centrifugal force" felt by a passenger is the absence of that force on their body: by inertia, they tend to continue in a straight line.
Circular orbits. If gravity supplies the centripetal force, then $\frac{GMm}{R^2} = \frac{mv^2}{R}$, so $v = \sqrt{\frac{GM}{R}}$ and $T = 2\pi\sqrt{\frac{R^3}{GM}}$, which is Kepler's third law for a circle. A geostationary satellite has a period of one sidereal day, which fixes $R \approx 42{,}164$ km from the Earth's center, about $35{,}786$ km above the surface.
Newton's law of gravitation in vector form. A planet at position $\mathbf{r}$ relative to the Sun (mass $M$) is pulled toward the Sun with
$$\mathbf{F} = -\frac{GMm}{r^3},\mathbf{r} = -\frac{GMm}{r^2},\mathbf{u}, \qquad r = \lvert \mathbf{r} \rvert, \quad \mathbf{u} = \frac{\mathbf{r}}{r}.$$
So $\mathbf{a} = -\frac{GM}{r^2}\mathbf{u}$: an inverse-square, central acceleration. (Treating the Sun as fixed is accurate because it is much heavier than any planet.)
Step 1: angular momentum is conserved, so the orbit is planar. Let $\mathbf{h} = \mathbf{r} \times \mathbf{v}$, the angular momentum per unit mass. By the cross product rule,
$$\frac{d\mathbf{h}}{dt} = \mathbf{v} \times \mathbf{v} + \mathbf{r} \times \mathbf{a} = \mathbf{0} + \mathbf{0},$$
because $\mathbf{a}$ is parallel to $\mathbf{r}$. So $\mathbf{h}$ is a constant vector. Since $\mathbf{r} \perp \mathbf{h}$ at every moment, the planet stays in the fixed plane through the Sun with normal $\mathbf{h}$ (assuming $\mathbf{h} \ne \mathbf{0}$, which excludes falling straight into the Sun). This holds for any central force.
Step 2: Kepler's second law. In polar coordinates in that plane, $\lvert \mathbf{h} \rvert = r^2\frac{d\theta}{dt}$, so the area swept per unit time, $\frac{dA}{dt} = \frac{1}{2}r^2\frac{d\theta}{dt} = \frac{h}{2}$, is constant: equal areas in equal times. (This is the computation of Section 10.3, now justified by Step 1.)
Step 3: Kepler's first law. Write $\mathbf{r} = r\mathbf{u}$, so $\mathbf{v} = r'\mathbf{u} + r\mathbf{u}'$, and
$$\mathbf{h} = \mathbf{r} \times \mathbf{v} = r\mathbf{u} \times (r'\mathbf{u} + r\mathbf{u}') = r^2(\mathbf{u} \times \mathbf{u}').$$
Now compute $\mathbf{a} \times \mathbf{h}$ using the "BAC minus CAB" rule (Section 11.1) and the facts $\mathbf{u} \cdot \mathbf{u} = 1$ and $\mathbf{u} \cdot \mathbf{u}' = 0$ (since $\mathbf{u}$ has constant length):
$$\mathbf{a} \times \mathbf{h} = -\frac{GM}{r^2}\mathbf{u} \times r^2(\mathbf{u} \times \mathbf{u}') = -GM\big[(\mathbf{u} \cdot \mathbf{u}')\mathbf{u} - (\mathbf{u} \cdot \mathbf{u})\mathbf{u}'\big] = GM,\mathbf{u}'.$$
Since $\mathbf{h}$ is constant, $\frac{d}{dt}(\mathbf{v} \times \mathbf{h}) = \mathbf{a} \times \mathbf{h} = GM,\mathbf{u}'$. Integrating both sides,
$$\mathbf{v} \times \mathbf{h} = GM,\mathbf{u} + \mathbf{c}$$
for a constant vector $\mathbf{c}$ lying in the orbital plane. Finally, dot both sides with $\mathbf{r} = r\mathbf{u}$. On the left, the scalar triple product can be cycled: $\mathbf{r} \cdot (\mathbf{v} \times \mathbf{h}) = (\mathbf{r} \times \mathbf{v}) \cdot \mathbf{h} = h^2$. On the right, $\mathbf{r} \cdot (GM\mathbf{u} + \mathbf{c}) = GMr + rc\cos\theta$, where $\theta$ is the angle between $\mathbf{c}$ and $\mathbf{r}$. Solving for $r$,
$$r = \frac{h^2}{GM + c\cos\theta} = \frac{h^2/GM}{1 + e\cos\theta}, \qquad e = \frac{c}{GM}.$$
This is the polar equation of a conic with a focus at the Sun and eccentricity $e$ (Section 10.3). The orbit is an ellipse when $e < 1$, a parabola when $e = 1$, and a hyperbola when $e > 1$. The constant vector
$$\mathbf{e} = \frac{\mathbf{v} \times \mathbf{h}}{GM} - \mathbf{u},$$
called the eccentricity vector (or Laplace-Runge-Lenz vector), points toward perihelion, and its length is the eccentricity.
Step 4: Kepler's third law follows from the first two as in Section 10.3: $T^2 = \frac{4\pi^2}{GM}a^3$.
The whole derivation uses only ideas from this chapter: the cross product, the product rules, the constant-length rule, and the scalar triple product.
Examples
Starting at $\mathbf{r}_0 = \langle 1, 0, 0 \rangle$ with $\mathbf{v}_0 = \langle 0, 1.2, 0 \rangle$ (in units where $GM = 1$), we get $\mathbf{h} = \langle 0, 0, 1.2 \rangle$ and $\mathbf{e} = \langle 0, 1.2, 0 \rangle \times \langle 0, 0, 1.2 \rangle - \langle 1, 0, 0 \rangle = \langle 0.44, 0, 0 \rangle$: an ellipse with $e = 0.44$, with perihelion at the starting point. The code simulates one orbital period with two integrators and checks the conserved quantities. The velocity Verlet method respects the geometry of the problem and returns to the starting point with $\mathbf{h}$ and $\mathbf{e}$ intact. Explicit Euler drifts: the orbit does not close, and the "conserved" quantities change.
import numpy as np
GM = 1.0 # units where GM = 1
def accel(r):
return -GM * r / np.linalg.norm(r)**3 # inverse-square, pointing at the origin
def invariants(r, v):
h = np.cross(r, v) # angular momentum per unit mass
e = np.cross(v, h) / GM - r / np.linalg.norm(r) # eccentricity (Laplace-Runge-Lenz) vector
E = 0.5 * v @ v - GM / np.linalg.norm(r) # energy per unit mass
return h, e, E
def verlet(r, v, dt, steps):
a = accel(r)
for _ in range(steps):
v_half = v + 0.5 * dt * a # half kick
r = r + dt * v_half # drift
a = accel(r)
v = v_half + 0.5 * dt * a # half kick with the new acceleration
return r, v
def euler(r, v, dt, steps):
for _ in range(steps):
r, v = r + dt * v, v + dt * accel(r)
return r, v
r0, v0 = np.array([1.0, 0.0, 0.0]), np.array([0.0, 1.2, 0.0])
h0, e0, E0 = invariants(r0, v0)
a_semi = -GM / (2 * E0) # semi-major axis from the energy
T = 2 * np.pi * np.sqrt(a_semi**3 / GM) # Kepler's third law
print(h0, np.round(e0, 4) + 0.0, round(E0, 4), round(a_semi, 4), round(T, 4))
# [0. 0. 1.2] [0.44 0. 0. ] -0.28 1.7857 14.9933
steps = 20000
dt = T / steps
for name, method in [("verlet", verlet), ("euler", euler)]:
r1, v1 = method(r0, v0, dt, steps) # one full period
h1, e1, E1 = invariants(r1, v1)
print(name, np.round(r1, 4) + 0.0, np.round(h1, 6) + 0.0, np.round(e1, 4) + 0.0, round(E1, 6))
# verlet [1. 0. 0.] [0. 0. 1.2] [0.44 0. 0. ] -0.28 back at perihelion
# euler [ 1.0007 -0.1001 0. ] [0. 0. 1.203188] [0.4417 0. 0. ] -0.278012 drifted
# Uniform circular motion: a = -omega^2 r, r perpendicular to v, |a| = v^2 / R
R, w = 2.0, 3.0
t = 0.7
r = R * np.array([np.cos(w * t), np.sin(w * t), 0])
v = R * w * np.array([-np.sin(w * t), np.cos(w * t), 0])
a = -R * w**2 * np.array([np.cos(w * t), np.sin(w * t), 0])
print(np.allclose(a, -w**2 * r), r @ v, np.linalg.norm(a), (v @ v) / R)
# True 8.881784197001252e-16 18.0 18.0
# Geostationary orbit: period = one sidereal day
GM_earth, T_day = 3.986004418e14, 86164.0905 # m^3/s^2, seconds
r_geo = (GM_earth * T_day**2 / (4 * np.pi**2)) ** (1 / 3) # from T^2 = 4 pi^2 r^3 / GM
print(round(r_geo / 1000), round((r_geo - 6378137) / 1000)) # 42164 35786 km from center, km altitude
Key Takeaways
- Uniform circular motion: $\mathbf{a} = -\omega^2\mathbf{r}$, centripetal, with magnitude $\frac{v^2}{R} = R\omega^2$.
- For any central force, $\mathbf{h} = \mathbf{r} \times \mathbf{v}$ is constant, so the motion is planar and sweeps equal areas in equal times.
- For an inverse-square force, $\mathbf{v} \times \mathbf{h} - GM\mathbf{u}$ is also constant. Dotting with $\mathbf{r}$ gives the conic $r = \frac{h^2/GM}{1 + e\cos\theta}$.
- In simulations, check conserved quantities ($\mathbf{h}$, $\mathbf{e}$, energy) to detect integration error. Symplectic methods such as velocity Verlet preserve them far better than explicit Euler.
🧪 Practice
- A centrifuge spins a sample at $3000$ revolutions per minute at a radius of $10$ cm. Find the sample's speed and its acceleration as a multiple of $g$.
- The International Space Station orbits about $420$ km above the Earth's surface (radius $6371$ km). Assuming a circular orbit, find its speed and period, using $GM = 3.986 \times 10^{14}$ m³/s².
- Using $\mathbf{v} \times \mathbf{h} = GM(\mathbf{u} + \mathbf{e})$, show that the speed at perihelion is $\frac{GM}{h}(1 + e)$ and at aphelion is $\frac{GM}{h}(1 - e)$. Check the ratio against the result of the Kepler's Laws exercise in Section 10.3.
- Interview-style: A colleague's orbital simulation uses explicit Euler, and after many orbits, the planet spirals outward. They propose halving the time step. What does that buy, and what would you recommend instead? Hint: Euler's error per orbit shrinks only in proportion to $\Delta t$, but the drift is systematic, so it still accumulates without bound over long runs. A symplectic integrator (velocity Verlet, leapfrog) keeps the energy error bounded and oscillating for all time, at the same cost per step. Monitor $\mathbf{h}$, $\mathbf{e}$, and the energy as diagnostics.
<a id="12-multivariable-differential-calculus"></a>