55  Vector-valued functions, \(f:R \rightarrow R^n\)

This section uses these add-on packages:

using CalculusWithJulia
using Plots
plotly()
using SymPy
using Roots
using LinearAlgebra
using QuadGK

and

import DifferentialEquations
import DifferentialEquations: ODEProblem, Tsit5

We discuss functions of a single variable that return a vector in \(R^n\). There are many parallels to univariate functions (when \(n=1\)) and differences.

55.1 Definition

A function \(\vec{f}: R \rightarrow R^n\), \(n > 1\) is called a vector-valued function. Some examples:

\[ \vec{f}(t) = \langle \sin(t), 2\cos(t) \rangle, \quad \vec{g}(t) = \langle \sin(t), \cos(t), t \rangle, \quad \vec{h}(t) = \langle 2, 3 \rangle + t \cdot \langle 1, 2 \rangle. \]

The components themselves are also functions of \(t\), in this case univariate functions. Depending on the context, it can be useful to view vector-valued functions as a function that returns a vector, or a vector of the component functions.

The above example functions have \(n\) equal \(2\), \(3\), and \(2\) respectively. We will see that many concepts of calculus for univariate functions (\(n=1\)) have direct counterparts.

(We use \(\vec{f}\) above to emphasize the return value is a vector, but will quickly drop that notation and let context determine if \(f\) refers to a scalar- or vector-valued function.)

55.2 Representation in Julia

In Julia, the representation of a vector-valued function is straightforward: we define a function of a single variable that returns a vector. For example, the three functions above would be represented by:

f(t) = [sin(t), 2*cos(t)]
g(t) = [sin(t), cos(t), t]
h(t) = [2, 3] + t * [1, 2]
h (generic function with 1 method)

For a given t, these evaluate to a vector. For example:

h(2)
2-element Vector{Int64}:
 4
 7

We can create a vector of functions, e.g., F = [cos, sin, identity], but calling this object, as in F(t), would require some work, such as t = 1; [f(t) for f in F] or 1 .|> F.

F = [cos, sin, identity]
[f(1) for f in F]
3-element Vector{Real}:
 0.5403023058681398
 0.8414709848078965
 1

or

1 .|> F
3-element Vector{Real}:
 0.5403023058681398
 0.8414709848078965
 1

55.3 Space curves

A vector-valued function is typically visualized as a curve. That is, for some range, \(a \leq t \leq b\) the set of points \(\{\vec{f}(t): a \leq t \leq b\}\) are plotted. If, say in \(n=2\), we have \(x(t)\) and \(y(t)\) as the component functions, then the graph would also be the parametric plot of \(x\) and \(y\). The term planar curve is common for the \(n=2\) case and space curve for the \(n \geq 3\) case.

This plot represents the vectors with their tails at the origin.

There is a convention for plotting the component functions to yield a parametric plot within the Plots package (e.g., plot(x, y, a, b)). This can be used to make polar plots, where x is t -> r(t)*cos(t) and y is t -> r(t)*sin(t).

However, we will use a different approach, as the component functions are not naturally produced from the vector-valued function.

In Plots, the command plot(xs, ys), where, say, xs=[x1, x2, ..., xn] and ys=[y1, y2, ..., yn], will make a connect-the-dot plot between corresponding pairs of points. As previously discussed, this can be used as an alternative to plotting a function through plot(f, a, b): first make a set of \(x\) values, say xs=range(a, b, length=100); then the corresponding \(y\) values, say ys = f.(xs); and then plotting through plot(xs, ys).

Similarly, were a third vector, zs, for \(z\) components used, plot(xs, ys, zs) will make a \(3\)-dimensional connect the dot plot

However, our representation of vector-valued functions naturally generates a vector of points: [[x1,y1], [x2, y2], ..., [xn, yn]], as this comes from broadcasting f over some time values. That is, for a collection of time values, ts the command f.(ts) will produce a vector of points. (Technically a vector of vectors, but points if you identify the \(2\)-\(d\) vectors as points.)

To get the xs and ys from this is conceptually easy: just iterate over all the points and extract the corresponding component. For example, to get xs we would have a command like [p[1] for p in f.(ts)]. Similarly, the ys would use p[2] in place of p[1]. The unzip function from the CalculusWithJulia package does this for us. The name comes from how the zip function in base Julia takes two vectors and returns a vector of the values paired off. This is the reverse. As previously mentioned, unzip uses the invert function of the SplitApplyCombine package to invert the indexing (the \(j\)th component of the \(i\)th point can be referenced by vs[i][j] or invert(vs)[j][i]).

Visually, we have unzip performing this reassociation:

[[x1, y1, z1],         (⌈x1⌉,  ⌈y1⌉, ⌈z1⌉,
 [x2, y2, z2],          |x2|, |y2|, |z2|,
 [x3, y3, z3],   -->    |x3|, |y3|, |z3|,
     ⋮                         ⋮
 [xn, yn, zn]]          ⌊xn⌋,  ⌊yn⌋, ⌊zn⌋ )

To turn a collection of vectors into separate arguments for a function, splatting (the ...) is used.

Finally, with these definitions, we can visualize the three functions we have defined.

Here we show the plot of f over the values between \(0\) and \(2\pi\) and also add a vector anchored at the origin defined by f(1).

ts = range(0, 2pi, length=200)
xs, ys = unzip(f.(ts))
plot(xs, ys)
arrow!([0, 0], f(1))

The trace of the plot is an ellipse. If we describe the components as \(\vec{f}(t) = \langle x(t), y(t) \rangle\), then we have \(x(t)^2 + y(t)^2/4 = 1\). That is, for any value of \(t\), the resulting point satisfies the equation \(x^2 + y^2/4 =1\) for an ellipse.

The plot of \(g\) needs \(3\)-dimensions to render. For most plotting backends, the following should work with no differences, save the additional vector is anchored in \(3\) dimensions now:

ts = range(0, 6pi, length=200)
plot(unzip(g.(ts))...) # use splatting to avoid xs,ys,zs = unzip(g.(ts))
arrow!([0, 0, 0], g(2pi))

Here the graph is a helix; three turns are plotted. If we write \(g(t) = \langle x(t), y(t), z(t) \rangle\), as the \(x\) and \(y\) values trace out a circle, the \(z\) value increases. When the graph is viewed from above, as below, we see only \(x\) and \(y\) components, and the view is circular.

ts = range(0, 6pi, length=200)
plot(unzip(g.(ts))..., camera=(0, 90))

The graph of \(h\) shows that this function parameterizes a line in space. The line segment for \(-2 \leq t \leq 2\) is shown below:

ts = range(-2, 2, length=200)
plot(unzip(h.(ts))...)

55.3.1 The plot_parametric function

While the unzip function is easy to understand as a function that reshapes data from one format into one that plot can use, its usage is a bit cumbersome. The CalculusWithJulia package provides a function plot_parametric which hides the use of unzip and the splatting within a function definition.

The function borrows a calling style for Makie. The interval to plot over is specified first using a..b notation (which specifies a closed interval in the IntervalSets package), then the function is specified. Additional keyword arguments are passed along to plot.

plot_parametric(-2..2, h)
Note

Defining plotting functions in Julia for Plots is facilitated by the RecipesBase package. There are two common choices: creating a new function for plotting, as is done with plot_parametric and plot_polar; or creating a new type so that plot can dispatch to an appropriate plotting method. The latter would also be a reasonable choice, but wasn’t taken here. In any case, each can be avoided by creating the appropriate values for xs and ys (and possibly zs).

Example

Familiarity with equations for lines, circles, and ellipses is important, as these fundamental geometric shapes are often building blocks in the description of other more complicated things.

The point-slope equation of a line, \(y = y_0 + m \cdot (x - x_0)\) finds an analog. The slope, \(m\), is replaced with a vector \(\vec{v}\) and the point, \((x_0, y_0)\) is replaced with a vector \(\vec{p}\) identified with a point in the plane. A parameterization would then be \(\vec{f}(t) = \vec{p} + (t - t_0) \vec{v}\). From this, we have \(\vec{f}(t_0) = \vec{p}\).

The unit circle is instrumental in introducing the trigonometric functions though the identification of an angle \(t\) with a point on the unit circle \((x,y)\) through \(y = \sin(t)\) and \(x=\cos(t)\). With this identification certain properties of the trigonometric functions are immediately seen, such as the period of \(\sin\) and \(\cos\) being \(2\pi\), or the angles for which \(\sin\) and \(\cos\) are positive or even increasing. Further, this gives a natural parameterization for a vector-valued function whose plot yields the unit circle, namely \(\vec{f}(t) = \langle \cos(t), \sin(t) \rangle\). This parameterization starts (at \(t=0\)) at the point \((1, 0)\). More generally, we might have additional parameters \(\vec{f}(t) = \vec{p} + R \cdot \langle \cos(\omega(t-t_0)), \sin(\omega(t-t_0)) \rangle\) to change the origin, \(\vec{p}\); the radius, \(R\); the starting angle, \(t_0\); and the rotational frequency, \(\omega\).

An ellipse has a slightly more general equation than a circle and in simplest forms may satisfy the equation \(x^2/a^2 + y^2/b^2 = 1\), where when \(a=b\) a circle is being described. A vector-valued function of the form \(\vec{f}(t) = \langle a\cdot\cos(t), b\cdot\sin(t) \rangle\) will trace out an ellipse.

The above description of an ellipse is useful, but it can also be useful to re-express the ellipse so that one of the foci is at the origin. With this, the ellipse can be given in polar coordinates through a description of the radius:

\[ r(\theta) = \frac{a (1 - e^2)}{1 + e \cos(\theta)}. \]

Here, \(a\) is the semi-major axis (\(a > b\)); \(e\) is the eccentricity given by \(b = a \sqrt{1 - e^2}\); and \(\theta\) a polar angle.

Using the conversion to Cartesian equations, we have \(\vec{f}(\theta) = \langle r(\theta) \cos(\theta), r(\theta) \sin(\theta)\rangle\).

For example:

a, ecc = 20, 3/4
f(t) = a*(1-ecc^2)/(1 + ecc*cos(t)) * [cos(t), sin(t)]
plot_parametric(0..2pi, f, legend=false)
scatter!([0],[0], markersize=4)
Example

The Spirograph is “… a geometric drawing toy that produces mathematical roulette curves of the variety technically known as hypotrochoids and epitrochoids. It was developed by British engineer Denys Fisher and first sold in \(1965\).” These can be used to make interesting geometrical curves.

Following Wikipedia: Consider a fixed outer circle \(C_o\) of radius \(R\) centered at the origin. A smaller inner circle \(C_i\) of radius \(r < R\) rolling inside \(C_o\) and is continuously tangent to it. \(C_i\) will be assumed never to slip on \(C_o\) (in a real Spirograph, teeth on both circles prevent such slippage). Now assume that a point \(A\) lying somewhere inside \(C_{i}\) is located a distance \(\rho < r\) from \(C_i\)’s center.

The center of the inner circle will move in a circular manner with radius \(R-r\). The fixed point on the inner circle will rotate about this center. The accumulated angle may be described by the angle the point of contact of the inner circle with the outer circle. Call this angle \(t\).

Suppose the outer circle is centered at the origin and the inner circle starts (\(t=0\)) with center \((R-r, 0)\) and rotates around counterclockwise. Then if the point of contact makes angle \(t\), the arc length along the outer circle is \(Rt\). The inner circle will have moved a distance \(r t'\) in the opposite direction, so \(Rt =-r t'\) and solving the angle will be \(t' = -(R/r)t\).

If the initial position of the fixed point is at \((\rho, 0)\) relative to the origin, then the following function will describe the motion:

\[ \vec{s}(t) = (R-r) \cdot \langle \cos(t), \sin(t) \rangle + \rho \cdot \langle \cos(-\frac{R}{r}t), \sin(-\frac{R}{r}t) \rangle. \]

To visualize this we first define a helper function to draw a circle at point \(P\) with radius \(R\):

circle!(P, R; kwargs...) = plot_parametric!(0..2pi, t -> P + R * [cos(t), sin(t)]; kwargs...)
circle! (generic function with 1 method)

Then we have this function to visualize the spirograph for different \(t\) values:

function spiro(t; r=2, R=5, rho=0.8*r)

    cent(t) = (R-r) * [cos(t), sin(t)]

    p = plot(legend=false, aspect_ratio=:equal)
    circle!([0,0], R, color=:blue)
    circle!(cent(t), r, color=:black)

    tp(t) = -R/r * t

    s(t) = cent(t) + rho * [cos(tp(t)), sin(tp(t))]
    plot_parametric!(0..t, s, color=:red)

    p
end
spiro (generic function with 1 method)

And we can see the trace for \(t=\pi\):

spiro(pi)

The point of contact is at \((-R, 0)\), as expected. Carrying this forward to a full circle’s worth is done through:

spiro(2pi)

The curve does not match up at the start. For that, a second time around the outer circle is needed:

spiro(4pi)

Whether the curve will have a period or not is decided by the ratio of \(R/r\) being rational or irrational.

Example

In 1935 Marcel Duchamp showed a collection of “Rotorelief” discs at a French fair for inventors. Disk number 10 is comprised of several nested, off-center circles on disk that would be rotated to give a sense of movement. To mimic the effect:

  • for each circle, \(3\) points where selected using a mouse from an image and their pixels recorded;
  • as \(3\) points determine a circle, the center and radius of each circle can be solved for
  • the exterior of the disc is drawn (the last specified circle below);
  • each nested circle is drawn after its center is rotated by \(\theta\) radian;
  • an animation captures the movement for display.
A Figure

Duchamp rotorelief

Example

Ivars Peterson described the carnival ride “tilt-a-whirl” as a chaotic system, whose equations of motion are presented in American Journal of Physics by Kautz and Huggard. The tilt-a-whirl has a platform that moves in a circle that also moves up and down. To describe the motion of a point on the platform assuming it has radius \(R\) and period \(T\) and rises twice in that period could be done with the function:

\[ \vec{u}(t) = \langle R \sin(2\pi t/T), R \cos(2\pi t/T), h + h \cdot \sin(2\pi t/ T) \rangle. \]

A passenger sits on a circular platform with radius \(r\) attached at some point on the larger platform. The dynamics of the person on the tilt-a-whirl depend on physics, but for simplicity, let’s assume the platform moves at a constant rate with period \(S\) and has no relative \(z\) component. The motion of the platform in relation to the point it is attached would be modeled by:

\[ \vec{v}(t) = \langle r \sin(2\pi t/S), r \cos(2\pi t/S), 0 \rangle. \]

And the motion relative to the origin would be the vector sum, or superposition:

\[ \vec{f}(t) = \vec{u}(t) + \vec{v}(t). \]

To visualize for some parameters, we have:

M, m = 25, 5
height = 5
S, T = 8, 2
outer(t) = [M * sin(2pi*t/T), M * cos(2pi*t/T), height*(1 +sin(2pi * (t-pi/2)/T))]
inner(t) = [m * sin(2pi*t/S), m * cos(2pi*t/S), 0]
f(t) = outer(t) + inner(t)
plot_parametric(0..8, f)

55.4 Limits and continuity

The definition of a limit for a univariate function is: For every \(\epsilon > 0\) there exists a \(\delta > 0\) such that if \(0 < |x-c| < \delta\) then \(|f(x) - L | < \epsilon\).

If the notion of “\(\vec{f}\) is close to \(L\)” is replaced by close in the sense of a norm, or vector distance, then the same limit definition can be used, with the new wording “… \(\| \vec{f}(x) - L \| < \epsilon\)”.

The notion of continuity is identical: \(\vec{f}(t)\) is continuous at \(t_0\) if \(\lim_{t \rightarrow t_0}\vec{f}(t) = \vec{f}(t_0)\). More informally \(\| \vec{f}(t) - \vec{f}(t_0)\| \rightarrow 0\).

A consequence of the triangle inequality is that a vector-valued function is continuous or has a limit if and only if its component functions do.

55.4.1 Derivatives

If \(\vec{f}(t)\) is vector valued, and \(\Delta t > 0\) then we can consider the vector:

\[ \vec{f}(t + \Delta t) - \vec{f}(t) \]

For example, if \(\vec{f}(t) = \langle 3\cos(t), 2\sin(t) \rangle\) and \(t=\pi/4\) and \(\Delta t = \pi/16\) we have this picture:

f(t) = [3cos(t), 2sin(t)]
t, Δt = pi/4, pi/16
df = f(t + Δt) - f(t)

plot(legend=false)
arrow!([0,0], f(t))
arrow!([0,0], f(t + Δt))
arrow!(f(t), df)

The length of the difference appears to be related to the length of \(\Delta t\), in a similar manner as the univariate derivative. The following limit defines the derivative of a vector-valued function:

\[ \vec{f}'(t) = \lim_{\Delta t \rightarrow 0} \frac{f(t + \Delta t) - f(t)}{\Delta t}. \]

The limit exists if the component limits do. The component limits are just the derivatives of the component functions. So, if \(\vec{f}(t) = \langle x(t), y(t) \rangle\), then \(\vec{f}'(t) = \langle x'(t), y'(t) \rangle\).

If the derivative is never \(\vec{0}\), the curve is called regular. For a regular curve the derivative is a tangent vector to the parameterized curve, akin to the case for a univariate function. We can use ForwardDiff to compute the derivative in the exact same manner as was done for univariate functions:

using ForwardDiff
D(f,n=1) = n > 1 ? D(D(f),n-1) : x -> ForwardDiff.derivative(f, float(x))
Base.adjoint(f::Function) = D(f)         # allow f' to compute derivative

(This is already done by the CalculusWithJulia package.)

We can visualize the tangential property through a graph:

f(t) = [3cos(t), 2sin(t)]
p = plot_parametric(0..2pi, f, legend=false, aspect_ratio=:equal)
for t in [1,2,3]
    arrow!(f(t), f'(t))   # add arrow with tail on curve, in direction of derivative
end
p

55.4.2 Symbolic representation

Were symbolic expressions used in place of functions, the vector-valued function would naturally be represented as a vector of expressions:

@syms 𝒕
𝒗vf = [cos(𝒕), sin(𝒕), 𝒕]

\(\left[\begin{smallmatrix}\cos{\left(𝒕 \right)}\\\sin{\left(𝒕 \right)}\\𝒕\end{smallmatrix}\right]\)

We will see working with these expressions is not identical to working with a vector-valued function.

To plot, we can avail ourselves of the the parametric plot syntax. The following expands to plot(cos(t), sin(t), t, 0, 2pi):

plot(𝒗vf..., 0, 2pi)

The unzip usage, as was done above, could be used, but it would be more trouble in this case.

To evaluate the function at a given value, say \(t=2\), we can use subs with broadcasting to substitute into each component:

subs.(𝒗vf, 𝒕=>2)

\(\left[\begin{smallmatrix}\cos{\left(2 \right)}\\\sin{\left(2 \right)}\\2\end{smallmatrix}\right]\)

Limits are performed component by component, and can also be defined by broadcasting, again with the need to adjust the values:

@syms Δ
limit.((subs.(𝒗vf, 𝒕 => 𝒕 + Δ) - 𝒗vf) / Δ, Δ => 0)

\(\left[\begin{smallmatrix}- \sin{\left(𝒕 \right)}\\\cos{\left(𝒕 \right)}\\1\end{smallmatrix}\right]\)

Derivatives, as was just done through a limit, are a bit more straightforward than evaluation or limit taking, as we won’t bump into the shape mismatch when broadcasting:

diff.(𝒗vf, 𝒕)

\(\left[\begin{smallmatrix}- \sin{\left(𝒕 \right)}\\\cos{\left(𝒕 \right)}\\1\end{smallmatrix}\right]\)

The second derivative, can be found through:

diff.(𝒗vf, 𝒕, 𝒕)

\(\left[\begin{smallmatrix}- \cos{\left(𝒕 \right)}\\- \sin{\left(𝒕 \right)}\\0\end{smallmatrix}\right]\)

55.4.3 Applications of the derivative

Here are some sample applications of the derivative.

Example: equation of the tangent line

The derivative of a vector-valued function is similar to that of a univariate function, in that it indicates a direction tangent to a curve. The point-slope form offers a straightforward parameterization. We have a point given through the vector-valued function and a direction given by its derivative. (After identifying a vector with its tail at the origin with the point that is the head of the vector.)

With this, the equation is simply \(\vec{tl}(t) = \vec{f}(t_0) + \vec{f}'(t_0) \cdot (t - t_0)\), where the dot indicates scalar multiplication.

Example: parabolic motion

In physics, we learn that the equation \(F=ma\) can be used to derive a formula for position, when acceleration, \(a\), is a constant. The resulting equation of motion is \(x = x_0 + v_0t + (1/2) at^2\). Similarly, if \(x(t)\) is a vector-valued position vector, and the second derivative, \(x''(t) =\vec{a}\), a constant, then we have: \(x(t) = \vec{x_0} + \vec{v_0}t + (1/2) \vec{a} t^2\).

For two dimensions, we have the force due to gravity acts downward, only in the \(y\) direction. The acceleration is then \(\vec{a} = \langle 0, -g \rangle\). If we start at the origin, with initial velocity \(\vec{v_0} = \langle 2, 3\rangle\), then we can plot the trajectory until the object returns to ground (\(y=0\)) as follows:

gravity = 9.8
x0, v0, a = [0,0], [2, 3], [0, -gravity]
xpos(t) = x0 + v0*t + (1/2)*a*t^2

t_0 = find_zero(t -> xpos(t)[2], (1/10, 100))  # find when y=0

plot_parametric(0..t_0, xpos)
Example: a tractrix

A tractrix, studied by Perrault, Newton, Huygens, and many others, is the curve along which an object moves when pulled in a horizontal plane by a line segment attached to a pulling point (Wikipedia). If the object is placed at \((a,0)\) and the puller at the origin, and the puller moves along the positive \(y\) axis, then the line will always be tangent to the curve and of fixed length, so determinable from the motion of the puller. In this example \(dy/dx = -\sqrt{a^2-x^2}/x\).

This is the key property: “Due to the geometrical way it was defined, the tractrix has the property that the segment of its tangent, between the asymptote and the point of tangency, has constant length \(a\).”

The tracks made by the front and rear bicycle wheels also have this same property and similarly afford a mathematical description. We follow Dunbar, Bosman, and Nooij from The Track of a Bicycle Back Tire below, though Levi and Tabachnikov and Foote, Levi, and Tabachnikov were also consulted. Let \(a\) be the distance between the front and back wheels, whose positions are parameterized by \(\vec{F}(t)\) and \(\vec{B}(t)\), respectively. The key property is the distance between the two is always \(a\), and, as the back wheel is always moving in the direction of the front wheel, we have \(\vec{B}'(t)\) is in the direction of \(\vec{F}(t) - \vec{B}(t)\), that is the vector \((\vec{F}(t)-\vec{B}(t))/a\) is a unit vector in the direction of the derivative of \(\vec{B}\). How long is the derivative vector? That would be answered by the speed of the back wheel, which is related to the velocity of the front wheel. But only the component of the velocity in the direction of \(\vec{F}(t)-\vec{B}(t)\), so the speed of the back wheel is the length of the projection of \(\vec{F}'(t)\) onto the unit vector \((\vec{F}(t)-\vec{B}(t))/a\), which is identified through the dot product.

Combined, this gives the following equations relating \(\vec{F}(t)\) to \(\vec{B}(t)\):

\[ s_B(t) = \vec{F}'(t) \cdot \frac{\vec{F}(t)-\vec{B}(t)}{a}, \quad \vec{B}'(t) = s_B(t) \frac{\vec{F}(t)-\vec{B}(t)}{a}. \]

This is a differential equation describing the motion of the back wheel in terms of the front wheel.

If the back wheel trajectory is known, the relationship is much easier, as the two differ by a vector of length \(a\) in the direction of \(\vec{B}'(t)\), or:

\[ F(t) = \vec{B}(t) + a \frac{\vec{B'(t)}}{\|\vec{B}'(t)\|}. \]

We don’t discuss when a differential equation has a solution, or if it is unique when it does, but note that the differential equation above may be solved numerically, in a manner somewhat similar to what was discussed in ODEs. Here we will use the DifferentialEquations package for finding numeric solutions.

We can define our equation as follows, using p to pass in the two parameters: the wheel-base length \(a\), and \(F(t)\), the parameterization of the front wheel in time:

function bicycle(dB, B, p, t)

  a, F = p   # unpack parameters

  speed =  F'(t) ⋅ (F(t) - B) / a
  dB[1], dB[2] = speed * (F(t) - B) / a

end
bicycle (generic function with 1 method)

Let’s consider a few simple cases first. We suppose \(a=1\) and the front wheel moves in a circle of radius \(3\). Here is how we can plot two loops:

t₀, t₁ = 0.0, 4pi

tspan₁ = (t₀, t₁)  # time span to consider

a₁ = 1
F₁(t) = 3 * [cos(t), sin(t)]
p₁ = (a₁, F₁)      # combine parameters

B₁0 = F₁(0) - [0, a₁]  # some initial position for the back
prob₁ = ODEProblem(bicycle, B₁0, tspan₁, p₁)

out₁ = solve(prob₁, reltol=1e-6, Tsit5())
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 141-element Vector{Float64}:
  0.0
  0.024824742235786373
  0.06751397058437489
  0.12257757397808876
  0.19299151337065923
  0.27111297357572545
  0.3857318027804938
  0.4674013218599202
  0.5643259834822881
  0.6554327060040652
  0.7518286045591882
  0.8474585608640102
  0.9451683751433828
  ⋮
 11.626710774978088
 11.715969045813116
 11.806663049962983
 11.898926531478901
 11.99240798013197
 12.086816057538613
 12.181697903789814
 12.27665249477995
 12.371291594129387
 12.465297210621507
 12.558395336376323
 12.566370614359172
u: 141-element Vector{Vector{Float64}}:
 [3.0, -1.0]
 [2.9999774774962056, -0.9255330157507721]
 [2.999561441536238, -0.7975914672002802]
 [2.997483559850816, -0.6329875818495893]
 [2.9906910543525567, -0.42353639734088794]
 [2.975659153342371, -0.19294877767647486]
 [2.9355038293214006, 0.14094436759864082]
 [2.89166595475666, 0.374747334436287]
 [2.821569496286955, 0.6465254432715674]
 [2.7369054701865503, 0.8949976424867925]
 [2.6270391498651806, 1.14884250900696]
 [2.49741327194199, 1.389663176947978]
 [2.3441433681310797, 1.622435087597618]
 ⋮
 [0.8123167712462607, -2.7092695725718965]
 [1.0505867865985716, -2.626074550668984]
 [1.2841118515417536, -2.5201303309780667]
 [1.5108364609525202, -2.391102949230913]
 [1.7274381999578818, -2.239633315522586]
 [1.930871223974605, -2.066817950409659]
 [2.1179957505055906, -1.874591723393362]
 [2.2861883302367114, -1.6653357183999136]
 [2.4333284917614986, -1.4418434716748734]
 [2.5579265192097687, -1.2070675530143982]
 [2.659062879683303, -0.964046039251787]
 [2.666666766573337, -0.9428088491842747]

The object out holds the answer. This object is callable, in that out(t) will return the numerically computed value for the answer to our equation at time point t.

To plot the two trajectories, we could use that out.u holds the \(x\) and \(y\) components of the computed trajectory, but more simply, we can just call out like a function.

plt₁ = plot_parametric(t₀..t₁, F₁, legend=false)
plot_parametric!(t₀..t₁, out₁,   linewidth=3)

## add the bicycle as a line segment at a few times along the path
for t in range(t₀, t₁, length=11)
    plot!(unzip([out₁(t), F₁(t)])..., linewidth=3, color=:black)
end
plt₁

That the rear wheel track appears shorter, despite the rear wheel starting outside the circle, is typical of bicycle tracks and also a reason to rotate tires on car, as the front ones move a bit more than the rear, so presumably wear faster.

Let’s look what happens if the front wheel wobbles back and forth following a sine curve. Repeating the above, only with \(F\) redefined, we have:

a₂ = 1
F₂(t) = [t, 2sin(t)]
p₂ = (a₂, F₂)

B₂0 = F₂(0) - [0, a₂]  # some initial position for the back
prob₂ = ODEProblem(bicycle, B₂0, tspan₁, p₂)

out₂ = solve(prob₂, reltol=1e-6, Tsit5())

plot_parametric(t₀..t₁, F₂, legend=false)
plot_parametric!(t₀..t₁, t -> out₂(t),  linewidth=3)

Again, the back wheel moves less than the front.

The motion of the back wheel need not be smooth, even if the motion of the front wheel is, as this curve illustrates:

a₃ = 1
F₃(t) = [cos(t), sin(t)] + [cos(2t), sin(2t)]
p₃ = (a₃, F₃)

B₃0 = F₃(0) - [0,a₃]
prob₃ = ODEProblem(bicycle, B₃0, tspan₁, p₃)

out₃ = solve(prob₃, reltol=1e-6, Tsit5())
plot_parametric(t₀..t₁, F₃, legend=false)
plot_parametric!(t₀..t₁, t -> out₃(t), linewidth=3)

The back wheel is moving backwards for part of the above trajectory.

This effect can happen even for a front wheel motion as simple as a circle when the front wheel radius is less than the wheelbase:

a₄ = 1
F₄(t) = a₄/3 * [cos(t), sin(t)]
p₄ = (a₄, F₄)

t₀₄, t₁₄ = 0.0, 25pi
tspan₄ = (t₀₄, t₁₄)

B₄0 = F₄(0) - [0, a₄]
prob₄ = ODEProblem(bicycle, B₄0, tspan₄, p₄)

out₄ = solve(prob₄, reltol=1e-6, Tsit5())
plot_parametric(t₀₄..t₁₄, F₄, legend=false, aspect_ratio=:equal)
plot_parametric!(t₀₄..t₁₄, t -> out₄(t), linewidth=3)

Later we will characterize when there are cusps in the rear-wheel trajectory.

55.5 Derivative rules

From the definition, as it is for univariate functions, for vector-valued functions \(\vec{f}, \vec{g}: R \rightarrow R^n\):

\[ [\vec{f} + \vec{g}]'(t) = \vec{f}'(t) + \vec{g}'(t), \quad\text{and } [a\vec{f}]'(t) = a \vec{f}'(t). \]

If \(a(t)\) is a univariate (scalar) function of \(t\), then a product rule holds:

\[ [a(t) \vec{f}(t)]' = a'(t)\vec{f}(t) + a(t)\vec{f}'(t). \]

If \(s\) is a univariate function, then the composition \(\vec{f}(s(t))\) can be differentiated. Each component would satisfy the chain rule, and consequently:

\[ \frac{d}{dt}\left(\vec{f}(s(t))\right) = \vec{f}'(s(t)) \cdot s'(t), \]

The dot being scalar multiplication by the derivative of the univariate function \(s\).

Vector-valued functions do not have multiplication or division defined for them, so there are no ready analogues of the product and quotient rule. However, the dot product and the cross product produce new functions that may have derivative rules available.

For the dot product, the combination \(\vec{f}(t) \cdot \vec{g}(t)\) we have a univariate function of \(t\), so we know a derivative is well defined. Can it be represented in terms of the vector-valued functions? In terms of the component functions, we have this calculation specific to \(n=2\), but that which can be generalized:

\[ \begin{align*} \frac{d}{dt}(\vec{f}(t) \cdot \vec{g}(t)) &= \frac{d}{dt}(f_1(t) g_1(t) + f_2(t) g_2(t))\\ &= f_1'(t) g_1(t) + f_1(t) g_1'(t) + f_2'(t) g_2(t) + f_2(t) g_2'(t)\\ &= f_1'(t) g_1(t) + f_2'(t) g_2(t) + f_1(t) g_1'(t) + f_2(t) g_2'(t)\\ &= \vec{f}'(t)\cdot \vec{g}(t) + \vec{f}(t) \cdot \vec{g}'(t). \end{align*} \]

Suggesting that a product rule like formula applies for dot products.

For the cross product, we let SymPy derive a formula for us.

@syms tₛ us()[1:3] vs()[1:3]
uₛ = tₛ .|> us  # evaluate each of us at t
vₛ = tₛ .|> vs

\(\left[\begin{smallmatrix}\operatorname{vs₁}{\left(tₛ \right)}\\\operatorname{vs₂}{\left(tₛ \right)}\\\operatorname{vs₃}{\left(tₛ \right)}\end{smallmatrix}\right]\)

Then the cross product has a derivative:

diff.(uₛ × vₛ, tₛ)

\(\left[\begin{smallmatrix}\operatorname{us₂}{\left(tₛ \right)} \frac{d}{d tₛ} \operatorname{vs₃}{\left(tₛ \right)} - \operatorname{us₃}{\left(tₛ \right)} \frac{d}{d tₛ} \operatorname{vs₂}{\left(tₛ \right)} - \operatorname{vs₂}{\left(tₛ \right)} \frac{d}{d tₛ} \operatorname{us₃}{\left(tₛ \right)} + \operatorname{vs₃}{\left(tₛ \right)} \frac{d}{d tₛ} \operatorname{us₂}{\left(tₛ \right)}\\- \operatorname{us₁}{\left(tₛ \right)} \frac{d}{d tₛ} \operatorname{vs₃}{\left(tₛ \right)} + \operatorname{us₃}{\left(tₛ \right)} \frac{d}{d tₛ} \operatorname{vs₁}{\left(tₛ \right)} + \operatorname{vs₁}{\left(tₛ \right)} \frac{d}{d tₛ} \operatorname{us₃}{\left(tₛ \right)} - \operatorname{vs₃}{\left(tₛ \right)} \frac{d}{d tₛ} \operatorname{us₁}{\left(tₛ \right)}\\\operatorname{us₁}{\left(tₛ \right)} \frac{d}{d tₛ} \operatorname{vs₂}{\left(tₛ \right)} - \operatorname{us₂}{\left(tₛ \right)} \frac{d}{d tₛ} \operatorname{vs₁}{\left(tₛ \right)} - \operatorname{vs₁}{\left(tₛ \right)} \frac{d}{d tₛ} \operatorname{us₂}{\left(tₛ \right)} + \operatorname{vs₂}{\left(tₛ \right)} \frac{d}{d tₛ} \operatorname{us₁}{\left(tₛ \right)}\end{smallmatrix}\right]\)

Admittedly, that isn’t very clear. With a peek at the answer, we show that the derivative is the same as the product rule would suggest (\(\vec{u}' \times \vec{v} + \vec{u} \times \vec{v}'\)):

diff.(uₛ × vₛ, tₛ) - (diff.(uₛ, tₛ) × vₛ + uₛ × diff.(vₛ, tₛ))

\(\left[\begin{smallmatrix}0\\0\\0\end{smallmatrix}\right]\)

In summary, these two derivative formulas hold for vector-valued functions \(R \rightarrow R^n\):

\[ \begin{align*} (\vec{u} \cdot \vec{v})' &= \vec{u}' \cdot \vec{v} + \vec{u} \cdot \vec{v}',\\ (\vec{u} \times \vec{v})' &= \vec{u}' \times \vec{v} + \vec{u} \times \vec{v}'. \end{align*} \]

Application. Circular motion and the tangent vector.

The parameterization \(\vec{r}(t) = \langle \cos(t), \sin(t) \rangle\) describes a circle. Characteristic of this motion is a constant radius, or in terms of a norm: \(\| \vec{r}(t) \| = c\). The norm squared, can be expressed in terms of the dot product:

\[ \| \vec{r}(t) \|^2 = \vec{r}(t) \cdot \vec{r}(t). \]

Differentiating this for the case of a constant radius yields the equation \(0 = [\vec{r}\cdot\vec{r}]'(t)\), which simplifies through the product rule and commutativity of the dot product to \(0 = 2 \vec{r}(t) \cdot \vec{r}'(t)\). That is, the two vectors are orthogonal to each other. This observation proves to be very useful, as will be seen.

Example: Kepler’s laws

Kepler’s laws of planetary motion are summarized by:

  • The orbit of a planet is an ellipse with the Sun at one of the two foci.
  • A line segment joining a planet and the Sun sweeps out equal areas during equal intervals of time.
  • The square of the orbital period of a planet is directly proportional to the cube of the semi-major axis of its orbit.

Kepler was a careful astronomer, and derived these laws empirically. We show next how to derive these laws using vector calculus assuming some facts on Newtonian motion, as postulated by Newton. This approach is borrowed from Joyce.

We adopt a sun-centered view of the universe, placing the sun at the origin and letting \(\vec{x}(t)\) be the position of a planet relative to this origin. We can express this in terms of a magnitude and direction through \(r(t) \hat{x}(t)\).

Newton’s law of gravitational force between the sun and this planet is then expressed by:

\[ \vec{F} = -\frac{G M m}{r^2} \hat{x}(t). \]

Newton’s famous law relating force and acceleration is

\[ \vec{F} = m \vec{a} = m \ddot{\vec{x}}. \]

Combining, Newton states \(\vec{a} = -(GM/r^2) \hat{x}\).

Now to show the first law. Consider \(\vec{x} \times \vec{v}\). It is constant, as:

\[ \begin{align*} (\vec{x} \times \vec{v})' &= \vec{x}' \times \vec{v} + \vec{x} \times \vec{v}'\\ &= \vec{v} \times \vec{v} + \vec{x} \times \vec{a}. \end{align*} \]

Both terms are \(\vec{0}\), as \(\vec{a}\) is parallel to \(\vec{x}\) by the above, and clearly \(\vec{v}\) is parallel to itself.

This says, \(\vec{x} \times \vec{v} = \vec{c}\) is a constant vector, meaning, the motion of \(\vec{x}\) must lie in a plane, as \(\vec{x}\) is always orthogonal to the fixed vector \(\vec{c}\).

Now, by differentiating \(\vec{x} = r \hat{x}\) we have:

\[ \begin{align*} \vec{v} &= \vec{x}'\\ &= (r\hat{x})'\\ &= r' \hat{x} + r \hat{x}', \end{align*} \]

and so

\[ \begin{align*} \vec{c} &= \vec{x} \times \vec{v}\\ &= (r\hat{x}) \times (r'\hat{x} + r \hat{x}')\\ &= r^2 (\hat{x} \times \hat{x}'). \end{align*} \]

From this, we can compute \(\vec{a} \times \vec{c}\):

\[ \begin{align*} \vec{a} \times \vec{c} &= (-\frac{GM}{r^2})\hat{x} \times r^2(\hat{x} \times \hat{x}')\\ &= -GM \hat{x} \times (\hat{x} \times \hat{x}') \\ &= GM (\hat{x} \times \hat{x}')\times \hat{x}. \end{align*} \]

The last line by anti-commutativity.

But, the triple cross product can be simplified through the identify \((\vec{u}\times\vec{v})\times\vec{w} = (\vec{u}\cdot\vec{w})\vec{v} - (\vec{v}\cdot\vec{w})\vec{u}\). So, the above becomes:

\[ \begin{align*} \vec{a} \times \vec{c} &= GM ((\hat{x}\cdot\hat{x})\hat{x}' - (\hat{x} \cdot \hat{x}')\hat{x})\\ &= GM (1 \hat{x}' - 0 \hat{x}). \end{align*} \]

Now, since \(\vec{c}\) is constant, we have:

\[ \begin{align*} (\vec{v} \times \vec{c})' &= (\vec{a} \times \vec{c})\\ &= GM \hat{x}'\\ &= (GM\hat{x})'. \end{align*} \]

The two sides have the same derivative, hence differ by a constant:

\[ \vec{v} \times \vec{c} = GM \hat{x} + \vec{d}. \]

As \(\vec{x}\) and \(\vec{v}\times\vec{c}\) lie in the same plane - orthogonal to \(\vec{c}\) - so does \(\vec{d}\). With a suitable re-orientation, so that \(\vec{d}\) is along the \(x\) axis, \(\vec{c}\) is along the \(z\)-axis, then we have \(\vec{c} = \langle 0,0,c\rangle\) and \(\vec{d} = \langle d ,0,0 \rangle\), and \(\vec{x} = \langle x, y, 0 \rangle\). Set \(\theta\) to be the angle, then \(\hat{x} = \langle \cos(\theta), \sin(\theta), 0\rangle\).

Now

\[ \begin{align*} c^2 &= \|\vec{c}\|^2 \\ &= \vec{c} \cdot \vec{c}\\ &= (\vec{x} \times \vec{v}) \cdot \vec{c}\\ &= \vec{x} \cdot (\vec{v} \times \vec{c})\\ &= r\hat{x} \cdot (GM\hat{x} + \vec{d})\\ &= GMr + r \hat{x} \cdot \vec{d}\\ &= GMr + rd \cos(\theta). \end{align*} \]

Solving, this gives the first law. That is, the radial distance is in the form of an ellipse:

\[ r = \frac{c^2}{GM + d\cos(\theta)} = \frac{c^2/(GM)}{1 + (d/GM) \cos(\theta)}. \]

Kepler’s second law can also be derived from vector calculus. This derivation follows that given at MIT OpenCourseWare and OpenCourseWare.

The second law states that the area being swept out during a time duration only depends on the duration of time, not the time. Let \(\Delta t\) be this duration. Then if \(\vec{x}(t)\) is the position vector, as above, we have the area swept out between \(t\) and \(t + \Delta t\) is visualized along the lines of:

x1(t) = [cos(t), 2 * sin(t)]
t0, t1, Delta = 1.0, 2.0, 1/10
plot_parametric(0..pi/2, x1)

arrow!([0,0], x1(t0)); arrow!([0,0], x1(t0 + Delta))
arrow!(x1(t0), x1(t0+Delta)- x1(t0), linewidth=5)

The area swept out, is basically the half the area of the parallelogram formed by \(\vec{x}(t)\) and \(\Delta \vec{x}(t) = \vec{x}(t + \Delta t) - \vec{x}(t)\). This area is \((1/2) (\vec{x} \times \Delta\vec{x}(t))\).

If we divide through by \(\Delta t\), and take a limit we have:

\[ \frac{dA}{dt} = \| \frac{1}{2}\lim_{\Delta t \rightarrow 0} (\vec{x} \times \frac{\vec{x}(t + \Delta t) - \vec{x}(t)}{\Delta t})\| = \frac{1}{2}\|\vec{x} \times \vec{v}\|. \]

But we saw above, that for the motion of a planet, that \(\vec{x} \times \vec{v} = \vec{c}\), a constant. This says, \(dA\) is a constant independent of \(t\), and consequently, the area swept out over a duration of time will not depend on the particular times involved, just the duration.

The third law relates the period to a parameter of the ellipse. We have from the above a strong suggestion that area of the ellipse can be found by integrating \(dA\) over the period, say \(T\). Assuming that is the case and letting \(a\) be the semi-major axis length, and \(b\) the semi-minor axis length, then

\[ \pi a b = \int_0^T dA = \int_0^T (1/2) \|\vec{x} \times \vec{v}\| dt = \| \vec{x} \times \vec{v}\| \frac{T}{2}. \]

As \(c = \|\vec{x} \times \vec{v}\|\) is a constant, this allows us to express \(c\) by: \(2\pi a b/T\).

But, we have

\[ r(\theta) = \frac{c^2}{GM + d\cos(\theta)} = \frac{c^2/(GM)}{1 + d/(GM) \cos(\theta)}. \]

So, \(e = d/(GM)\) and \(a (1 - e^2) = c^2/(GM)\). Using \(b = a \sqrt{1-e^2}\) we have:

\[ a(1-e^2) = c^2/(GM) = (\frac{2\pi a b}{T})^2 \frac{1}{GM} = \frac{(2\pi)^2}{GM} \frac{a^2 (a^2(1-e^2))}{T^2}, \]

or after cancelling \((1-e^2)\) from each side:

\[ T^2 = \frac{(2\pi)^2}{GM} \frac{a^4}{a} = \frac{(2\pi)^2}{GM} a^3. \]

The above shows how Newton might have derived Kepler’s observational facts. Next we show, that assuming the laws of Kepler can anticipate Newton’s equation for gravitational force. This follows Wikipedia.

Now let \(\vec{r}(t)\) be the position of the planet relative to the Sun at the origin, in two dimensions (we used \(\vec{x}(t)\) above). Assume \(\vec{r}(0)\) points in the \(x\) direction. Write \(\vec{r} = r \hat{r}\). Define \(\hat{\theta}(t)\) to be the mapping from time \(t\) to the angle defined by \(\hat{r}\) through the unit circle.

Then we express the velocity (\(\dot{\vec{r}}\)) and acceleration (\(\ddot{\vec{r}}\)) in terms of the orthogonal vectors \(\hat{r}\) and \(\hat{\theta}\), as follows:

\[ \frac{d}{dt}(r \hat{r}) = \dot{r} \hat{r} + r \dot{\hat{r}} = \dot{r} \hat{r} + r \dot{\theta}\hat{\theta}. \]

The last equality from expressing \(\hat{r}(t) = \hat{r}(\theta(t))\) and using the chain rule, noting \(d(\hat{r}(\theta))/d\theta = \hat{\theta}\).

Continuing,

\[ \frac{d^2}{dt^2}(r \hat{r}) = (\ddot{r} \hat{r} + \dot{r} \dot{\hat{r}}) + (\dot{r} \dot{\theta}\hat{\theta} + r \ddot{\theta}\hat{\theta} + r \dot{\theta}\dot{\hat{\theta}}). \]

Noting, similar to above, \(\dot{\hat{\theta}} = d\hat{\theta}/dt = d\hat{\theta}/d\theta \cdot d\theta/dt = -\dot{\theta} \hat{r}\) we can express the above in terms of \(\hat{r}\) and \(\hat{\theta}\) as:

\[ \vec{a} = \frac{d^2}{dt^2}(r \hat{r}) = (\ddot{r} - r (\dot{\theta})^2) \hat{r} + (r\ddot{\theta} + 2\dot{r}\dot{\theta}) \hat{\theta}. \]

That is, in general, the acceleration has a radial component and a transversal component.

Kepler’s second law says that the area increment over time is constant (\(dA/dt\)), but this area increment is approximated by the following wedge in polar coordinates: \(dA = (1/2) r \cdot rd\theta\). We have then \(dA/dt = r^2 \dot{\theta}\) is constant.

Differentiating, we have:

\[ 0 = \frac{d(r^2 \dot{\theta})}{dt} = 2r\dot{r}\dot{\theta} + r^2 \ddot{\theta}, \]

which is the transversal component of the acceleration times \(r\), as decomposed above. This means, that the acceleration of the planet is completely towards the Sun at the origin.

Kepler’s first law, relates \(r\) and \(\theta\) through the polar equation of an ellipse:

\[ r = \frac{p}{1 + \epsilon \cos(\theta)}. \]

Expressing in terms of \(p/r\) and differentiating in \(t\) gives:

\[ -\frac{p \dot{r}}{r^2} = -\epsilon\sin(\theta) \dot{\theta}. \]

Or

\[ p\dot{r} = \epsilon\sin(\theta) r^2 \dot{\theta} = \epsilon \sin(\theta) C, \]

For a constant \(C\), used above, as the second law implies \(r^2 \dot{\theta}\) is constant. (This constant can be expressed in terms of parameters describing the ellipse.)

Differentiating again in \(t\), gives:

\[ p \ddot{r} = C\epsilon \cos(\theta) \dot{\theta} = C\epsilon \cos(\theta)\frac{C}{r^2}. \]

So \(\ddot{r} = (C^2 \epsilon / p) \cos{\theta} (1/r^2)\).

The radial acceleration from above is:

\[ \ddot{r} - r (\dot{\theta})^2 = (C^2 \epsilon/p) \cos{\theta} \frac{1}{r^2} - r\frac{C^2}{r^4} = \frac{C^2}{pr^2}(\epsilon \cos(\theta) - \frac{p}{r}). \]

Using \(p/r = 1 + \epsilon\cos(\theta)\), we have the radial acceleration is \(C^2/p \cdot (1/r^2)\). That is the acceleration, is proportional to the inverse square of the position, and using the relation between \(F\), force, and acceleration, we see the force on the planet follows the inverse-square law of Newton.

55.6 Moving frames of reference

In the last example, it proved useful to represent vectors in terms of other unit vectors, in that case \(\hat{r}\) and \(\hat{\theta}\). Here we discuss a coordinate system defined intrinsically by the motion along the trajectory of a curve.

Let \(r(t)\) be a smooth vector-valued function in \(R^3\). It gives rise to a space curve, through its graph. This curve has tangent vector \(\vec{r}'(t)\), indicating the direction of travel along \(\vec{r}\) as \(t\) increases. The length of \(\vec{r}'(t)\) depends on the parameterization of \(\vec{r}\), as for any increasing, differentiable function \(s(t)\), the composition \(\vec{r}(s(t))\) will have derivative, \(\vec{r}'(s(t)) s'(t)\), having the same direction as \(\vec{r}'(t)\) (at suitably calibrated points), but not the same magnitude, the factor of \(s(t)\) being involved.

To discuss properties intrinsic to the curve, the unit vector is considered:

\[ \hat{T}(t) = \frac{\vec{r}'(t)}{\|\vec{r}'(t)\|}. \]

The function \(\hat{T}(t)\) is the unit tangent vector. An assumption of regularity ensures the denominator is never \(0\).

Now define the unit normal, \(\hat{N}(t)\), by:

\[ \hat{N}(t) = \frac{\hat{T}'(t)}{\| \hat{T}'(t) \|}. \]

Since \(\|\hat{T}(t)\| = 1\), a constant, it must be that \(\hat{T}'(t) \cdot \hat{T}(t) = 0\), that is, the \(\hat{N}\) and \(\hat{T}\) are orthogonal.

Finally, define the binormal, \(\hat{B}(t) = \hat{T}(t) \times \hat{N}(t)\). At each time \(t\), the three unit vectors are orthogonal to each other. They form a moving coordinate system for the motion along the curve that does not depend on the parameterization.

We can visualize this, for example along a Viviani curve, as is done in a Wikipedia animation:

function viviani(t, a=1)
    [a*(1-cos(t)), a*sin(t), 2a*sin(t/2)]
end


Tangent(t) = viviani'(t)/norm(viviani'(t))
Normal(t) = Tangent'(t)/norm(Tangent'(t))
Binormal(t) = Tangent(t) × Normal(t)

p = plot(legend=false)
plot_parametric!(-2pi..2pi, viviani)

t0, t1 = -pi/3, pi/2 + 2pi/5
r0, r1 = viviani(t0), viviani(t1)
arrow!(r0, Tangent(t0)); arrow!(r0, Binormal(t0)); arrow!(r0, Normal(t0))
arrow!(r1, Tangent(t1)); arrow!(r1, Binormal(t1)); arrow!(r1, Normal(t1))
p

The curvature of a \(3\)-dimensional space curve is defined by:

The curvature: For a \(3-D\) curve the curvature is defined by:

\(\kappa = \frac{\| r'(t) \times r''(t) \|}{\| r'(t) \|^3}.\)

For \(2\)-dimensional space curves, the same formula applies after embedding a \(0\) third component. It can also be expressed directly as

\[ \kappa = (x'y''-x''y')/\|r'\|^3. \quad (r(t) =\langle x(t), y(t) \rangle) \]

Curvature can also be defined as derivative of the tangent vector, \(\hat{T}\), when the curve is parameterized by arc length, a topic still to be taken up. The vector \(\vec{r}'(t)\) is the direction of motion, whereas \(\vec{r}''(t)\) indicates how fast and in what direction this is changing. For curves with little curve in them, the two will be nearly parallel and the cross product small (reflecting the presence of \(\cos(\theta)\) in the definition). For “curvy” curves, \(\vec{r}''\) will be in a direction opposite of \(\vec{r}'\) to the \(\cos(\theta)\) term in the cross product will be closer to \(1\).

Let \(\vec{r}(t) = k \cdot \langle \cos(t), \sin(t), 0 \rangle\). This will have curvature:

@syms k::positive t::real
r1 = k * [cos(t), sin(t), 0]
norm(diff.(r1,t) × diff.(r1,t,t)) / norm(diff.(r1,t))^3 |> simplify

\(\frac{1}{k}\)

For larger circles (bigger \(\|k\|\)) there is less curvature. The limit being a line with curvature \(0\).

If a curve is imagined to have a tangent “circle” (second order Taylor series approximation), then the curvature of that circle matches the curvature of the curve.

The torsion, \(\tau\), of a space curve (\(n=3\)), is a measure of how sharply the curve is twisting out of the plane of curvature.

The torsion is defined for smooth curves by

The torsion:

\(\tau = \frac{(\vec{r}' \times \vec{r}'') \cdot \vec{r}'''}{\|\vec{r}' \times \vec{r}''\|^2}.\)

For the torsion to be defined, the cross product \(\vec{r}' \times \vec{r}''\) must be non zero, that is the two must not be parallel or zero.

Example: Tubular surface

This last example comes from a collection of several examples provided by Discourse user @empet to illustrate plotlyjs. We adopt it to Plots with some minor changes below.

The task is to illustrate a space curve, \(c(t)\), using a tubular surface. At each time point \(t\), assume the curve has tangent, \(e_1\); normal, \(e_2\); and binormal, \(e_3\). (This assumes the defining derivatives exist and are non-zero and the cross product in the torsion is non zero.) The tubular surface is a circle of radius \(\epsilon\) in the plane determined by the normal and binormal. This curve would be parameterized by \(r(t,u) = c(t) + \epsilon (e_2(t) \cdot \cos(u) + e_3(t) \cdot \sin(u))\) for varying \(u\).

The Frenet-Serret equations setup a system of differential equations driven by the curvature and torsion. We use the DifferentialEquations package to solve this equation for two specific functions and a given initial condition. The equations when expanded into coordinates become \(12\) different equations:

# e₁, e₂, e₃, (x,y,z)
function Frenet_eq!(du, u, p, s)  #system of ODEs
    κ, τ = p
    du[1] =  κ(s) * u[4]               # e₁′ = κ ⋅ e₂
    du[2] =  κ(s) * u[5]
    du[3] =  κ(s) * u[6]
    du[4] = -κ(s) * u[1] + τ(s) * u[7] # e₂′ = - κ ⋅ e₁ + τ ⋅ e₃
    du[5] = -κ(s) * u[2] + τ(s) * u[8]
    du[6] = -κ(s) * u[3] + τ(s) * u[9]
    du[7] = -τ(s) * u[4]               # e₃′ = - τ ⋅ e₂
    du[8] = -τ(s) * u[5]
    du[9] = -τ(s) * u[6]
    du[10] = u[1]                      # c′ = e₁
    du[11] = u[2]
    du[12] = u[3]
end
Frenet_eq! (generic function with 1 method)

The last set of equations describe the motion of the spine. It follows from specifying the tangent to the curve is \(e_1\), as desired; it is parameterized by arc length, as \(\mid c'(t) \mid = 1\).

Following the example of @empet, we define a curvature function and torsion function, the latter a constant:

κ(s) = 3 * sin(s/10) * sin(s/10)
τ(s) = 0.35
τ (generic function with 1 method)

The initial condition and time span are set with:

e₁₀, e₂₀, e₃₀ = [1,0,0], [0,1,0], [0,0,1]
u₀ = [0, 0, 0]
u0 = vcat(e₁₀, e₂₀, e₃₀, u₀) # initial condition for the system of ODE
t_span = (0.0,  150.0)          # time interval for solution
(0.0, 150.0)

With this set up, the problem can be solved:

prob = ODEProblem(Frenet_eq!, u0, t_span, (κ, τ))
sol = solve(prob, Tsit5());

The “spine” is the center axis of the tube and is the \(10\)th, \(11\)th, and \(12\)th coordinates:

spine(t) = sol(t)[10:12]
spine (generic function with 1 method)

The tangent, normal, and binormal can be similarly defined using the other \(9\) indices:

e₁(t) = sol(t)[1:3]
e₂(t) = sol(t)[4:6]
e₃(t) = sol(t)[7:9]
e₃ (generic function with 1 method)

We fix a small time range and show the trace of the spine and the frame at a single point in time:

a_0, b_0 = 50, 60
ts_0 = range(a_0, b_0, length=251)

t_0 = (a_0 + b_0) / 2
ϵ = 1/5

plot_parametric(a_0..b_0, spine)

arrow!(spine(t_0), e₁(t_0))
arrow!(spine(t_0), e₂(t_0))
arrow!(spine(t_0), e₃(t_0))

r_0(t, θ) = spine(t) + ϵ * (e₂(t)*cos(θ) + e₃(t)*sin(θ))
plot_parametric!(0..2pi, θ -> r_0(t_0, θ))

The ϵ value determines the radius of the tube; we see it above as the radius of the drawn circle. The function r for a fixed t traces out such a circle centered at a point on the spine. For a fixed θ, the function r describes a line on the surface of the tube paralleling the spine.

The tubular surface is now ready to be rendered along the entire time span using a pattern for parametrically defined surfaces:

ts = range(t_span..., length=1001)
θs = range(0, 2pi, length=100)
surface(unzip(r_0.(ts, θs'))...)