# Vector-valued functions, $f:R \rightarrow R^n$


{{< include ../_common_code.qmd >}}

This section uses these add-on packages:


```{julia}
using CalculusWithJulia
using Plots
plotly()
using SymPy
using Roots
using LinearAlgebra
using QuadGK
```

and


```{julia}
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.


## 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.)


## 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:


```{julia}
f(t) = [sin(t), 2*cos(t)]
g(t) = [sin(t), cos(t), t]
h(t) = [2, 3] + t * [1, 2]
```

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


```{julia}
h(2)
```

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`.


```{julia}
F = [cos, sin, identity]
[f(1) for f in F]
```

or


```{julia}
1 .|> F
```

## 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:


```{verbatim}
[[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)`.


```{julia}
#| hold: true
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:


```{julia}
#| hold: true
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.


```{julia}
#| hold: true
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:


```{julia}
#| hold: true
ts = range(-2, 2, length=200)
plot(unzip(h.(ts))...)
```

### 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`.


```{julia}
plot_parametric(-2..2, h)
```

:::{.callout-note}
## 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:


```{julia}
#| hold: true
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](https://en.wikipedia.org/wiki/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$:


```{julia}
circle!(P, R; kwargs...) = plot_parametric!(0..2pi, t -> P + R * [cos(t), sin(t)]; kwargs...)
```

Then we have this function to visualize the spirograph for different $t$ values:


```{julia}
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
```

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


```{julia}
spiro(pi)
```

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


```{julia}
spiro(2pi)
```

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


```{julia}
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](https://arthur.io/art/marcel-duchamp/rotorelief-no-10-cage-modele-depose-verso) 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.


```{julia}
#| echo: false
let
# https://exploratorium.tumblr.com/post/33140874462/marcel-duchamp-rotoreliefs-duchamp-recognized

# coordinates and colors selected by gimp from
    # https://arthur.io/art/marcel-duchamp/rotorelief-no-10-cage-modele-depose-verso
    gr()
    circs = [466 548 513 505 556 554 # x₁,y₁,x₂,y₂,x₂,y₃
             414 549 511 455 595 549
             365 545 507 408 635 548
             319 541 506 361 673 546
             277 543 509 317 711 546
             236 539 507 272 747 551
             201 541 504 230 781 550
             166 541 503 189 816 544
             140 542 499 153 848 538
             116 537 496 119 879 538
             96 539 501 90 905 534
             81 530 500 67 930 530
             72 525 498 51 949 529
             66 520 500 36 966 527
             60 515 499 25 982 526
             35 509 499 11 1004 525 # outer edge, c₀
             ]

    greenblue= RGB(8/100, 58/100, 53/100)
    grey = RGB(76/100, 74/100, 72/100)
    white = RGB(88/100, 85/100, 81/100)

    # solve for center of circle, radius for each
    @syms h::positive k::positive r::positive
    function solve_i(i)
        eqs = [(p[1] - h)^2 + (p[2]-k)^2 ~ r^2 for
		       p ∈ (circs[i,1:2], circs[i,3:4], circs[i,5:6])]
        d = solve(tuple(eqs...))[1]
        (x=float(d[h]), y=float(d[k]), r=float(d[r]))
    end
    c₀, cs... = solve_i.(16:-1:1) # c₀ is centered

    function duchamp_rotorelief_10(θ)
        p = plot(legend=false,
		         axis=nothing, xaxis=false, yaxis=false,
		         aspect_ratio=:equal)

        O = [c₀.x, c₀.y]
        θ̂ = [cos(θ), sin(θ)]

	    circle!(O, c₀.r, # outer ring is c₀
		        linewidth=2,
		        color=grey, fill=white,
		        seriestype=:shape)

        for (i,c) ∈ enumerate(cs)   # add nested rings
	        rᵢ = sqrt((c₀.x - c.x)^2+(c₀.y - c.y)^2)
		    P = O + rᵢ * θ̂  # rotate about origin by θ
		    circle!(P, c.r,
			        linewidth = i == 1 ? 1 : i <= 3 ? 2 : 3,
			        color=greenblue)
        end

        p

    end

    # animate using Plots.@animate macro
    anim = @animate for θ ∈ range(0, -2π, length=60)
	    duchamp_rotorelief_10(θ)
    end

    fname = tempname() * ".gif"
    gif(anim, fname, fps = 40)

    plotly()
    ImageFile(fname, "Duchamp rotorelief")
end
```

```{julia}
#| echo: false
#import PlutoUI
#PlutoUI.LocalResource(fname) # to display w/in Pluto
nothing
```

##### Example


[Ivars Peterson](http://www.phschool.com/science/science_news/articles/tilt_a_whirl.html) described the carnival ride "tilt-a-whirl" as a chaotic system, whose equations of motion are presented in [American Journal of Physics](https://doi.org/10.1119/1.17742) 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:


```{julia}
#| hold: true
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)
```

## 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.


### 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:


```{julia}
#| hold: true
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:


```{julia}
#| eval: false
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:


```{julia}
#| hold: true
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
```

### Symbolic representation


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


```{julia}
@syms 𝒕
𝒗vf = [cos(𝒕), sin(𝒕), 𝒕]
```

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)`:


```{julia}
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:


```{julia}
subs.(𝒗vf, 𝒕=>2)
```

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


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

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:


```{julia}
diff.(𝒗vf, 𝒕)
```

The second derivative, can be found through:


```{julia}
diff.(𝒗vf, 𝒕, 𝒕)
```

### 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:


```{julia}
#| hold: true
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)
```

```{julia}
#| echo: false
# https://en.wikipedia.org/wiki/Tractrix
# https://sinews.siam.org/Details-Page/a-bike-and-a-catenary
# https://www.math.psu.edu/tabachni/talks/BicycleDouble.pdf
# https://www.tandfonline.com/doi/abs/10.4169/amer.math.monthly.120.03.199
# https://projecteuclid.org/download/pdf_1/euclid.em/1259158427
nothing
```

##### Example: a tractrix


A [tractrix](https://en.wikipedia.org/wiki/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](https://doi.org/10.2307/2691097) from *The Track of a Bicycle Back Tire* below, though [Levi and Tabachnikov](https://projecteuclid.org/download/pdf_1/euclid.em/1259158427) and [Foote, Levi, and Tabachnikov](https://www.tandfonline.com/doi/abs/10.4169/amer.math.monthly.120.03.199) 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](../ODEs/odes.html). 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:


```{julia}

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
```

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:


```{julia}
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())
```

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.


```{julia}
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:


```{julia}
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:


```{julia}
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:


```{julia}
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.


## 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.


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

Then the cross product has a derivative:


```{julia}
diff.(uₛ × vₛ, tₛ)
```

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}'$):


```{julia}
diff.(uₛ × vₛ, tₛ) - (diff.(uₛ, tₛ) × vₛ + uₛ × diff.(vₛ, tₛ))
```

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](https://tinyurl.com/y38wragh)'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](https://mathcs.clarku.edu/~djoyce/ma131/kepler.pdf).


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](https://ocw.mit.edu/courses/mathematics/18-02sc-multivariable-calculus-fall-2010/1.-vectors-and-matrices/part-c-parametric-equations-for-curves/session-21-keplers-second-law/MIT18_02SC_MNotes_k.pdf) and [OpenCourseWare](https://ocw.mit.edu/courses/mathematics/18-02sc-multivariable-calculus-fall-2010/index.htm).


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:


```{julia}
#| hold: true
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](https://en.wikipedia.org/wiki/Kepler%27s_laws_of_planetary_motion#Planetary_acceleration).


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.


## 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](https://tinyurl.com/y4lo29mv) curve, as is done in a [Wikipedia](https://en.wikipedia.org/wiki/Frenet%E2%80%93Serret_formulas) animation:


```{julia}
#| hold: true
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:


```{julia}
#| hold: true
@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
```

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](https://en.wikipedia.org/wiki/Torsion_of_a_curve), $\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

```{julia}
#| echo: false
gr();
```

This last example comes from a collection of several [examples](https://github.com/empet/3D-Viz-with-PlotlyJS.jl/blob/main/5-Tubular-surface.ipynb) 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:


```{julia}
# 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
```

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:


```{julia}
κ(s) = 3 * sin(s/10) * sin(s/10)
τ(s) = 0.35
```

The initial condition and time span are set with:


```{julia}
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
```

With this set up, the problem can be solved:


```{julia}
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:


```{julia}
spine(t) = sol(t)[10:12]
```

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


```{julia}
e₁(t) = sol(t)[1:3]
e₂(t) = sol(t)[4:6]
e₃(t) = sol(t)[7:9]
```

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



```{julia}
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:


```{julia}
#| hold: true
ts = range(t_span..., length=1001)
θs = range(0, 2pi, length=100)
surface(unzip(r_0.(ts, θs'))...)
```

```{julia}
#| echo: false
plotly();
```

## Arc length


In [Arc length](../integrals/arc_length.html) there is a discussion of how to find the arc length of a parameterized curve in $2$ dimensions. The general case is discussed by [Destafano](https://randomproofs.files.wordpress.com/2010/11/arc_length.pdf) who shows:


> *Arc-length*: if a curve $C$ is parameterized by a smooth function $\vec{r}(t)$ over an interval $I$, then the arc length of $C$ is:
>
> $$
> \int_I \| \vec{r}'(t)  \| dt.
> $$



If we associate $\vec{r}'(t)$ with the velocity, then this is the integral of the speed (the magnitude of the velocity).


Let $I=[a,b]$ and $s(t): [v,w] \rightarrow [a,b]$ such that $s$ is increasing and differentiable. Then $\vec{\phi} = \vec{r} \circ s$ will have


$$
\text{arc length} =
\int_v^w \| \vec{\phi}'(t)\| dt =
\int_v^w \| \vec{r}'(s(t))\| s'(t) dt =
\int_a^b \| \vec{r}'(u) \| du,
$$

by a change of variable $u=s(t)$. As such the arc length is a property of the curve and not the parameterization of the curve.


For some parameterization, we can define


$$
s(t) = \int_0^t \| \vec{r}'(u) \| du
$$

Then by the fundamental theorem of calculus, $s(t)$ is non-decreasing. If $\vec{r}'$ is assumed to be non-zero and continuous (regular), then $s(t)$ has a derivative and an inverse which is monotonic. Using the inverse function $s^{-1}$ to change variables ($\vec{\phi} = \vec{r} \circ s^{-1}$) has


$$
\int_0^c \| \phi'(t) \| dt =
\int_{s^{-1}(0)}^{s^{-1}(c)} \| \vec{r}'(u) \| du =
s(s^{-1}(c)) - s(s^{-1}(0)) =
c
$$

That is, the arc length from $[0,c]$ for $\phi$ is just $c$; the curve $C$ is parameterized by arc length.


##### Example


Viviani's curve is the intersection of sphere of radius $a$ with a cylinder of radius $a$. A parameterization was given previously by:


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

The curve is traced out over the interval $[0, 4\pi]$. We try to find the arc-length:


```{julia}
#| hold: true
@syms t::positive a::positive
speed = simplify(norm(diff.(viviani(t, a), t)))
integrate(speed, (t, 0, 4*PI))
```

We see that the answer depends linearly on $a$, but otherwise is a constant expressed as an integral. We use `QuadGk` to provide a numeric answer for the case $a=1$:


```{julia}
quadgk(t -> norm(viviani'(t)), 0, 4pi)
```

##### Example


Very few parameterized curves admit a closed-form expression for parameterization by arc-length. Let's consider the helix expressed by $\langle a\cos(t), a\sin(t), bt\rangle$, as this does allow such a parameterization.


```{julia}
@syms aₕ::positive bₕ::positive tₕ::positive alₕ::positive
helix = [aₕ * cos(tₕ), aₕ * sin(tₕ), bₕ * tₕ]
speed = simplify( norm(diff.(helix, tₕ)) )
s = integrate(speed, (tₕ, 0, alₕ))
```

So `s` is a linear function. We can re-parameterize by:


```{julia}
eqnₕ = subs.(helix, tₕ => alₕ/sqrt(aₕ^2 + bₕ^2))
```

To see that the speed, $\| \vec{\phi}' \|$, is constantly $1$:


```{julia}
simplify(norm(diff.(eqnₕ, alₕ)))
```

From this, we have the arc length is:


$$
\int_0^t \| \vec{\phi}'(u) \| du = \int_0^t 1 du = t
$$

---


Parameterizing by arc-length is only explicitly possible for a few examples, however knowing it can be done in theory, is important. Some formulas are simplified, such as the tangent, normal, and binormal. Let $\vec{r}(s)$ be parameterized by arc length, then:


$$
\hat{T}(s)= \vec{r}'(s) / \| \vec{r}'(s) \| = \hat{r}'(s),\quad
\hat{N}(s) = \hat{T}'(s) / \| \hat{T}'(s)\| = \hat{T}'(s)/\kappa,\quad
\hat{B} = \hat{T} \times \hat{N},
$$

As before, but further, we have if $\kappa$ is the curvature and $\tau$ the torsion, these relationships expressing the derivatives with respect to $s$ in terms of the components in the frame:



\begin{align*}
\hat{T}'(s) &=                    &\kappa \hat{N}(s)  &\\
\hat{N}'(s) &= -\kappa \hat{T}(s) &                   &+ \tau \hat{B}(s)\\
\hat{B}'(s) &=                    &-\tau \hat{N}(s)   &
\end{align*}


These are the [Frenet-Serret](https://en.wikipedia.org/wiki/Frenet%E2%80%93Serret_formulas) formulas.


##### Example


Continuing with our parameterization of a helix by arc length, we can compute the curvature and torsion by differentiation:


```{julia}
gammaₕ = subs.(helix, tₕ => alₕ/sqrt(aₕ^2 + bₕ^2))   # gamma parameterized by arc length
@syms uₕ::positive
gammaₕ₁ = subs.(gammaₕ, alₕ .=> uₕ)                  # u is arc-length parameterization
```

```{julia}
Tₕ = diff.(gammaₕ₁, uₕ)
norm(Tₕ)  |> simplify
```

The length is one, as the speed of a curve parameterized by arc-length is 1.


```{julia}
outₕ = diff.(Tₕ, uₕ)
```

This should be $\kappa \hat{N}$, so we do:


```{julia}
κₕ = norm(outₕ) |> simplify
Normₕ = outₕ / κₕ
κₕ, Normₕ
```

Interpreting, $a$ is the radius of the circle and $b$ how tight the coils are. If $a$ gets much larger than $b$, then the curvature is like $1/a$, just as with a circle. If $b$ gets very big, then the trajectory looks more stretched out and the curvature gets smaller.


To find the torsion, we find, $\hat{B}$ then differentiate:


```{julia}
Bₕ = Tₕ × Normₕ
outₕ₁ = diff.(Bₕ, uₕ)
τₕ = norm(outₕ₁)
```

This looks complicated, as does `Norm`:


```{julia}
Normₕ
```

However, the torsion, up to a sign, simplifies nicely:


```{julia}
τₕ |> simplify
```

Here, when $b$ gets large, the curve looks more and more "straight" and the torsion decreases. Similarly, if $a$ gets big, the torsion decreases.


##### Example


[Levi and Tabachnikov](https://projecteuclid.org/download/pdf_1/euclid.em/1259158427) consider the trajectories of the front and rear bicycle wheels. Recall the notation previously used: $\vec{F}(t)$ for the front wheel, and $\vec{B}(t)$ for the rear wheel trajectories. Consider now their parameterization by arc length, using $u$ for the arc-length parameter for $\vec{F}$ and $v$ for $\vec{B}$. We define $\alpha(u)$ to be the steering angle of the bicycle. This can be found as the angle between the tangent vector of the path of $\vec{F}$ with the vector $\vec{F} - \vec{B}$. Let $\kappa$ be the curvature of the front wheel and $k$ the curvature of the back wheel.


```{julia}
#| echo: false
let
    gr()
a = 1
F(t) = [cos(pi/2 - t), 2sin(pi/2-t)]
p = (a, F)

t0, t1 = -pi/6, pi/2.75
tspan = (t0, t1)

t = 7pi/6
B0 = F(t0) + a*[cos(t), sin(t)]
prob = ODEProblem(bicycle, B0, tspan, p)

out = solve(prob, reltol=1e-6, Tsit5())
plt = plot_parametric(t0..t1, F, linewidth=3,
                      xticks=nothing, yticks=nothing, border=:none,
                      legend=false, aspect_ratio=:equal)
plot_parametric!(t0..t1, t -> out(t),  linewidth=3)

t = pi/4
arrow!(out(t), 2*(F(t) - out(t)))
plot!(unzip([out(t), F(t)])..., linewidth=2)
arrow!(F(t), F'(t)/norm(F'(t)))
Fphat(t) = F'(t)/norm(F'(t))
arrow!( F(t), -Fphat'(t)/norm(Fphat'(t)))
using LaTeXStrings
annotate!([(-.5,1.5,L"k"),
(.775,1.55,L"\kappa"),
(.85, 1.3, L"\alpha")])
plotly()
plt
end
```

Levi and Tabachnikov prove in their Proposition 2.4:



\begin{align*}
\kappa(u) &= \frac{d\alpha(u)}{du} + \frac{\sin(\alpha(u))}{a},\\
|\frac{dv}{du}| &= |\cos(\alpha)|, \quad \text{and}\\
k &= \frac{\tan(\alpha)}{a}.
\end{align*}


The first equation relates the steering angle with the curvature. If the steering angle is not changed ($d\alpha/du=0$) then the curvature is constant and the motion is circular. It will be greater for larger angles (up to $\pi/2$). As the curvature is the reciprocal of the radius, this means the radius of the circular trajectory will be smaller. For the same constant steering angle, the curvature will be smaller for longer wheelbases, meaning the circular trajectory will have a larger radius. For cars, which have similar dynamics, this means longer wheelbase cars will take more room to make a U-turn.


The second equation may be interpreted in ratio of arc lengths. The infinitesimal arc length of the rear wheel is proportional to that of the front wheel only scaled down by $\cos(\alpha)$. When $\alpha=0$ - the bike is moving in a straight line - and the two are the same. At the other extreme - when $\alpha=\pi/2$ - the bike must be pivoting on its rear wheel and the rear wheel has no arc length. This cosine, is related to the speed of the back wheel relative to the speed of the front wheel, which was used in the initial differential equation.


The last equation, relates the curvature of the back wheel track to the steering angle of the front wheel. When $\alpha=\pm\pi/2$, the rear-wheel curvature, $k$, is infinite, resulting in a cusp (no circle with non-zero radius will approximate the trajectory). This occurs when the front wheel is steered orthogonal to the direction of motion. As was seen in previous graphs of the trajectories, a cusp can happen for quite regular front wheel trajectories.


To derive the first one, we have previously noted that when a curve is parameterized by arc length, the curvature is more directly computed: it is the magnitude of the derivative of the tangent vector. The tangent vector is of unit length, when parametrized by arc length. This implies its derivative will be orthogonal. If $\vec{r}(t)$ is a parameterization by arc length, then the curvature formula simplifies as:



\begin{align*}
\kappa(s) &= \frac{\| \vec{r}'(s) \times \vec{r}''(s) \|}{\|\vec{r}'(s)\|^3} \\
&= \frac{\| \vec{r}'(s) \times \vec{r}''(s) \|}{1}                           \\
&= \| \vec{r}'(s) \| \| \vec{r}''(s) \| \sin(\theta)                         \\
&= 1 \| \vec{r}''(s) \| 1 = \| \vec{r}''(s) \|.
\end{align*}


So in the above, the curvature is $\kappa = \| \vec{F}''(u) \|$ and $k = \|\vec{B}''(v)\|$.


On the figure, the tangent vector $\vec{F}'(u)$ is drawn, along with this unit vector rotated by $\pi/2$. We call these, for convenience, $\vec{U}$ and $\vec{V}$. We have $\vec{U} = \vec{F}'(u)$ and $\vec{V} = -(1/\kappa) \vec{F}''(u)$.


The key decomposition, is to express a unit vector in the direction of the line segment, as the vector $\vec{U}$ rotated by $\alpha$ degrees. Mathematically, this is usually expressed in matrix notation, but more explicitly by


$$
\langle \cos(\alpha) \vec{U}_1 + \sin(\alpha) \vec{U}_2,
-\sin(\alpha) \vec{U}_1 + \cos(\alpha) \vec{U}_2 \rangle =
\vec{U} \cos(\alpha) + \vec{V} \sin(\alpha).
$$

With this, the mathematical relationship between $F$ and $B$ is just a multiple of this unit vector:


$$
\vec{B}(u) = \vec{F}(u) - a \vec{U} \cos(\alpha) - a \vec{V} \sin(\alpha).
$$

It must be that the tangent line of $\vec{B}$ is parallel to $\vec{U} \cos(\alpha) + \vec{V} \sin(\alpha)$. To utilize this, we differentiate $\vec{B}$ using the facts that $\vec{U}' = -\kappa \vec{V}$ and $\vec{V}' = \kappa \vec{U}$. These coming from $\vec{U} = \vec{F}'$ and so it's derivative in $u$ has magnitude yielding the curvature, $\kappa$, and direction orthogonal to $\vec{U}$.



\begin{align*}
\vec{B}'(u) &= \vec{F}'(u)
-a \vec{U}' \cos(\alpha) -a \vec{U} (-\sin(\alpha)) \alpha'
-a \vec{V}' \sin(\alpha) - a \vec{V} \cos(\alpha) \alpha'\\
& =  \vec{U}
+a (\kappa) \vec{V} \cos(\alpha) + a \vec{U} \sin(\alpha) \alpha' -
a (\kappa) \vec{U} \sin(\alpha) - a \vec{V} \cos(\alpha) \alpha' \\
&= \vec{U}
+ a(\alpha' - \kappa) \sin(\alpha) \vec{U}
- a(\alpha' - \kappa) \cos(\alpha)\vec{V}.
\end{align*}


Extend the $2$-dimensional vectors to $3$ dimensions, by adding a zero $z$ component, then:



\begin{align*}
\vec{0} &= (\vec{U}
+ a(\alpha' - \kappa) \sin(\alpha) \vec{U}
+ a(\alpha' - \kappa) \cos(\alpha)\vec{V}) \times
(\vec{U} \cos(\alpha) +  \vec{V} \sin(\alpha)) \\
&= (\vec{U} \times \vec{V}) \sin(\alpha) +
a(\alpha' - \kappa) \sin(\alpha) \vec{U} \times  \vec{V} \sin(\alpha) -
a(\alpha' - \kappa) \cos(\alpha)\vec{V} \times \vec{U} \cos(\alpha) \\
&= (\sin(\alpha) + a(\alpha'-\kappa) \sin^2(\alpha) +
a(\alpha'-\kappa) \cos^2(\alpha)) \vec{U} \times \vec{V} \\
&= (\sin(\alpha) + a (\alpha' - \kappa)) \vec{U} \times \vec{V}.
\end{align*}


The terms $\vec{U} \times\vec{U}$ and $\vec{V}\times\vec{V}$ being $\vec{0}$, due to properties of the cross product. This says the scalar part must be $0$, or


$$
\frac{\sin(\alpha)}{a} + \alpha' = \kappa.
$$

As for the second equation, from the expression for $\vec{B}'(u)$, after setting $a(\alpha'-\kappa) = -\sin(\alpha)$:



\begin{align*}
\|\vec{B}'(u)\|^2
&= \| (1 -\sin(\alpha)\sin(\alpha)) \vec{U} +\sin(\alpha)\cos(\alpha) \vec{V} \|^2\\
&= \| \cos^2(\alpha) \vec{U} +\sin(\alpha)\cos(\alpha) \vec{V} \|^2\\
&= (\cos^2(\alpha))^2 + (\sin(\alpha)\cos(\alpha))^2\quad\text{using } \vec{U}\cdot\vec{V}=0\\
&= \cos^2(\alpha)(\cos^2(\alpha) + \sin^2(\alpha))\\
&= \cos^2(\alpha).
\end{align*}


From this $\|\vec{B}(u)\| = |\cos(\alpha)\|$. But $1 = \|d\vec{B}/dv\| = \|d\vec{B}/du \| \cdot |du/dv|$ and $|dv/du|=|\cos(\alpha)|$ follows.



## Evolutes and involutes


Following [Fuchs](https://doi.org/10.4169/amer.math.monthly.120.03.217) we discuss a geometric phenomenon known and explored by Huygens, and likely earlier. We stick to the two-dimensional case, Fuchs extends this to three dimensions. The following figure


```{julia}
#| echo: false
Xₑ(t)= 2 * cos(t)
Yₑ(t) = sin(t)
rₑ(t) = [Xₑ(t), Yₑ(t)]
unit_vec(x) = x / norm(x)
plot(legend=false, aspect_ratio=:equal)
ts = range(0, 2pi, length=50)
for t in ts
    Pₑ, Vₑ = rₑ(t), unit_vec([-Yₑ'(t), Xₑ'(t)])
    plot_parametric!(-4..4, x -> Pₑ + x*Vₑ)
end
plot!(Xₑ, Yₑ, 0, 2pi, linewidth=5)
```

is that of an ellipse with many *normal* lines drawn to it. The normal lines appear to intersect in a somewhat diamond-shaped curve. This curve is the evolute of the ellipse. We can characterize this using the language of planar curves.


Consider a parameterization of a curve by arc-length, $\vec\gamma(s) = \langle u(s), v(s) \rangle$. The unit *tangent* to this curve is $\vec\gamma'(s) = \hat{T}(s) = \langle u'(s), v'(s) \rangle$ and by simple geometry the unit *normal* will be $\hat{N}(s) = \langle -v'(s), u'(s) \rangle$. At a time $t$, a line through the curve parameterized by $\vec\gamma$ is given by $l_t(a) = \vec\gamma(t) + a \hat{N}(t)$.


Consider two nearby points $t$ and $t+\epsilon$ and the intersection of $l_t$ and $l_{t+\epsilon}$. That is, we need points $a$ and $b$ with: $l_t(a) = l_{t+\epsilon}(b)$. Setting the components equal, this is:



\begin{align*}
u(t) - av'(t) &= u(t+\epsilon) - bv'(t+\epsilon) \\
v(t) + au'(t) &= v(t+\epsilon) + bu'(t+\epsilon).
\end{align*}


This is a linear equation in two unknowns ($a$ and $b$) which can be solved. Here is the value for `a`:


```{julia}
@syms u() v() t epsilon w
@syms a b
γ(t) = [u(t),v(t)]
n(t) = subs.(diff.([-v(w), u(w)], w), w.=>t)
l(a, t) = γ(t) + a * n(t)
out = solve(l(a, t) - l(b, t+epsilon), [a,b])
out[a]
```

Letting $\epsilon \rightarrow 0$ we get an expression for $a$ that will describe the evolute at time $t$ in terms of the function $\gamma$. Looking at the expression above, we can see that dividing the *numerator* by $\epsilon$ and taking a limit will yield $u'(t)^2 + v'(t)^2$. If the *denominator* has a limit after dividing by $\epsilon$, then we can find the description sought. Pursuing this leads to:



\begin{align*}
\frac{u'(t) v'(t+\epsilon) - v'(t) u'(t+\epsilon)}{\epsilon}
&= \frac{u'(t) v'(t+\epsilon) -u'(t)v'(t) + u'(t)v'(t)- v'(t) u'(t+\epsilon)}{\epsilon} \\
&= \frac{u'(t)(v'(t+\epsilon) -v'(t))}{\epsilon} + \frac{(u'(t)- u'(t+\epsilon))v'(t)}{\epsilon},
\end{align*}


which in the limit will give $u'(t)v''(t) - u''(t) v'(t)$. All told, in the limit as $\epsilon \rightarrow 0$ we get



\begin{align*}
a &= \frac{u'(t)^2 + v'(t)^2}{u'(t)v''(t) - v'(t) u''(t)} \\
&= 1/(\|\vec\gamma'\|\kappa) \\
&= 1/(\|\hat{T}\|\kappa) \\
&= 1/\kappa,
\end{align*}


with $\kappa$ being the curvature of the planar curve. That is, the evolute of $\vec\gamma$ is described by:


$$
\vec\beta(s) = \vec\gamma(s) + \frac{1}{\kappa(s)}\hat{N}(s).
$$

Revisualizing:


```{julia}
rₑ₃(t) = [2cos(t), sin(t), 0]
Tangent(r, t) = unit_vec(r'(t))
Normal(r, t) = unit_vec((𝒕 -> Tangent(r, 𝒕))'(t))
curvature(r, t) = norm(r'(t) × r''(t) ) / norm(r'(t))^3

plot_parametric(0..2pi, t -> rₑ₃(t)[1:2], legend=false, aspect_ratio=:equal)
plot_parametric!(0..2pi, t -> (rₑ₃(t) + Normal(rₑ₃, t)/curvature(rₑ₃, t))[1:2])
```

We computed the above illustration using $3$ dimensions (hence the use of `[1:2]...`) as the curvature formula is easier to express. Recall, the curvature also appears in the [Frenet-Serret](https://en.wikipedia.org/wiki/Frenet%E2%80%93Serret_formulas) formulas: $d\hat{T}/ds = \kappa \hat{N}$ and $d\hat{N}/ds = -\kappa \hat{T}+ \tau \hat{B}$. In a planar curve, as under consideration, the binormal is $\vec{0}$. This allows the computation of $\vec\beta(s)'$:



\begin{align*}
\vec{\beta}' &= \frac{d(\vec\gamma + (1/ \kappa) \hat{N})}{ds}\\
&= \hat{T} + (-\frac{\kappa '}{\kappa ^2}\hat{N} + \frac{1}{\kappa} \hat{N}')\\
&= \hat{T} - \frac{\kappa '}{\kappa ^2}\hat{N} + \frac{1}{\kappa} (-\kappa \hat{T})\\
&= - \frac{\kappa '}{\kappa ^2}\hat{N}.
\end{align*}


We see $\vec\beta'$ is zero (the curve is non-regular) when $\kappa'(s) = 0$. The curvature changes from increasing to decreasing, or vice versa at each of the $4$ crossings of the major and minor axes - there are $4$ non-regular points, and we see $4$ cusps in the evolute.


The curve parameterized by $\vec{r}(t) = 2(1 - \cos(t)) \langle \cos(t), \sin(t)\rangle$ over $[0,2\pi]$ is cardiod. It is formed by rolling a circle of radius $r$ around another similar sized circle. The following graphically shows the evolute is a smaller cardiod (one-third the size). For fun, the evolute of the evolute is drawn:


```{julia}
function evolute(r)
    t -> r(t) + 1/curvature(r, t) * Normal(r, t)
end
```

```{julia}
#| hold: trie
r(t) = 2*(1 - cos(t)) * [cos(t), sin(t), 0]

plot(legend=false, aspect_ratio=:equal)
plot_parametric!(0..2pi, t -> r(t)[1:2])
plot_parametric!(0..2pi, t -> evolute(r)(t)[1:2])
plot_parametric!(0..2pi, t -> ((evolute∘evolute)(r)(t))[1:2])
```

---


If $\vec\beta$ is *the* **evolute** of $\vec\gamma$, then $\vec\gamma$ is *an* **involute** of $\beta$. For a given curve, there is a parameterized family of involutes. While this definition has a pleasing self-referentialness, it doesn't have an immediately clear geometric interpretation. For that, consider the image of a string of fixed length $a$ attached to the curve $\vec\gamma$ at some point $t_0$. As this curve wraps around the curve traced by $\vec\gamma$ it is held taut so that it makes a tangent at the point of contact. The end of the string will trace out a curve and this is the trace of an *involute*.


```{julia}
#| hold: true
r(t) = [t, cosh(t)]
t0, t1 = -2, 0
a = t1

beta(r, t) = r(t) - Tangent(r, t) * quadgk(t -> norm(r'(t)), a, t)[1]

p = plot_parametric(-2..2, r, legend=false)
plot_parametric!(t0..t1, t -> beta(r, t))
for t in range(t0,-0.2, length=4)
    arrow!(r(t), -Tangent(r, t) * quadgk(t -> norm(r'(t)), a, t)[1])
    scatter!(unzip([r(t)])...)
end
p
```

This lends itself to this mathematical description, if $\vec\gamma(t)$ parameterizes the planar curve, then an involute for $\vec\gamma(t)$ is described by:


$$
\vec\beta(t) = \vec\gamma(t) + \left((a - \int_{t_0}^t \| \vec\gamma'(t)\| dt) \hat{T}(t)\right),
$$

where $\hat{T}(t) = \vec\gamma'(t)/\|\vec\gamma'(t)\|$ is the unit tangent vector. The above uses two parameters ($a$ and $t_0$), but only one is needed, as there is an obvious redundancy (a point can *also* be expressed by $t$ and the shortened length of string). [Wikipedia](https://en.wikipedia.org/wiki/Involute) uses this definition for $a$ and $t$ values in an interval $[t_0, t_1]$:


$$
\vec\beta_a(t) = \vec\gamma(t) - \frac{\vec\gamma'(t)}{\|\vec\gamma'(t)\|}\int_a^t \|\vec\gamma'(t)\| dt.
$$

If $\vec\gamma(s)$ is parameterized by arc length, then this simplifies quite a bit, as the unit tangent is just $\vec\gamma'(s)$ and the remaining arc length just $(s-a)$:



\begin{align*}
\vec\beta_a(s) &= \vec\gamma(s) - \vec\gamma'(s) (s-a) \\
&=\vec\gamma(s) - \hat{T}_{\vec\gamma}(s)(s-a).\quad (a \text{ is the arc-length parameter})
\end{align*}


With this characterization, we see several properties:


  * From $\vec\beta_a'(s) = \hat{T}(s) - (\kappa(s) \hat{N}(s) (s-a) + \hat{T}(s)) = -\kappa_{\vec\gamma}(s) \cdot (s-a) \cdot \hat{N}_{\vec\gamma}(s)$,  the involute is *not* regular at $s=a$, as its derivative is zero.
  * As $\vec\beta_a(s) = \vec\beta_0(s) + a\hat{T}(s)$, the family of curves is parallel.
  * The evolute of $\vec\beta_a(s)$, $s$ the arc-length parameter of $\vec\gamma$, can be shown to be $\vec\gamma$. This requires more work:


The evolute for $\vec\beta_a(s)$ is:


$$
\vec\beta_a(s) + \frac{1}{\kappa_{\vec\beta_a}(s)}\hat{N}_{\vec\beta_a}(s).
$$

In the following we show that:



\begin{align*}
\kappa_{\vec\beta_a}(s) &= 1/(s-a),\\
\hat{N}_{\vec\beta_a}(s) &= \hat{T}_{\vec\beta_a}'(s)/\|\hat{T}_{\vec\beta_a}'(s)\| = -\hat{T}_{\vec\gamma}(s).
\end{align*}


The first shows in a different way that when $s=a$ the curve is not regular, as the curvature fails to exists. In the above figure, when the involute touches $\vec\gamma$, there will be a cusp.


With these two identifications and using $\vec\gamma'(s) = \hat{T}_{\vec\gamma(s)}$, we have the evolute simplifies to



\begin{align*}
\vec\beta_a(s) + \frac{1}{\kappa_{\vec\beta_a}(s)}\hat{N}_{\vec\beta_a}(s)
&=
\vec\gamma(s) + \vec\gamma'(s)(s-a)  + \frac{1}{\kappa_{\vec\beta_a}(s)}\hat{N}_{\vec\beta_a}(s) \\
&=
\vec\gamma(s) + \hat{T}_{\vec\gamma}(s)(s-a) + \frac{1}{1/(s-a)} (-\hat{T}_{\vec\gamma}(s)) \\
&= \vec\gamma(s).
\end{align*}


That is the evolute of an involute of $\vec\gamma(s)$ is $\vec\gamma(s)$.


We have:



\begin{align*}
\beta_a(s) &= \vec\gamma - \vec\gamma'(s)(s-a)\\
\beta_a'(s) &= -\kappa_{\vec\gamma}(s)(s-a)\hat{N}_{\vec\gamma}(s)\\
\beta_a''(s) &= (-\kappa_{\vec\gamma}(s)(s-a))' \hat{N}_{\vec\gamma}(s) + (-\kappa_{\vec\gamma}(s)(s-a))(-\kappa_{\vec\gamma}\hat{T}_{\vec\gamma}(s)),
\end{align*}


the last line by the Frenet-Serret  formula for *planar* curves which show $\hat{T}'(s) = \kappa(s) \hat{N}$ and $\hat{N}'(s) = -\kappa(s)\hat{T}(s)$.


To compute the curvature of $\vec\beta_a$, we need to compute both:



\begin{align*}
\| \vec\beta' \|^3 &= |\kappa^3 (s-a)^3|\\
\| \vec\beta' \times \vec\beta'' \| &= |\kappa(s)^3 (s-a)^2|,
\end{align*}


the last line using both $\hat{N}\times\hat{N} = \vec{0}$ and $\|\hat{N}\times\hat{T}\| = 1$. The curvature then is $\kappa_{\vec\beta_a}(s) = 1/(s-a)$.


Using the formula for $\vec\beta'$ above, we get $\hat{T}_\beta(s)=\hat{N}_{\vec\gamma}(s)$ so $\hat{T}_\beta(s)' = -\kappa_{\vec\gamma}(s) \hat{T}_{\vec\gamma}(s)$ with unit vector just $\hat{N}_{\vec\beta_a} = -\hat{T}_{\vec\gamma}(s)$.


---


Show that *an* involute of the cycloid $\vec{r}(t) = \langle t - \sin(t), 1 - \cos(t) \rangle$ is also a cycloid. We do so graphically:


```{julia}
#| hold: true
r(t) = [t - sin(t), 1 - cos(t)]
## find *involute*: r - r'/|r'| * int(|r'|, a, t)
t0, t1, a = 2PI, PI, PI
@syms t::real
rp = diff.(r(t), t)
speed = 2sin(t/2)

ex = r(t) - rp/speed * integrate(speed, t)

plot_parametric(0..4pi, r, legend=false)
plot_parametric!(0..4pi, u -> float.(subs.(ex, t .=> u)))
```

The expression `ex` is secretly `[t + sin(t), 3 + cos(t)]`, another cycloid.


##### Example: goats


An old problem of calculus is called the [goat problem](https://en.wikipedia.org/wiki/Goat_problem). This formulation – with horses – is from $1748$:


> Observing a horse tied to feed in a gentlemen’s park, with one end of a rope to his fore foot, and the other end to one of the circular iron rails, inclosing a pond, the circumference of which rails being $160$ yards, equal to the length of the rope, what quantity of ground at most, could the horse feed?



Let $r$ be the radius of a circle and for concreteness we position it at $(-r, 0)$. Let $R$ be the length of a rope, and suppose $R \ge 2\pi r$. (It is equal in the problem). Then the question can be rephrased as what is *twice* the area suggested by this graphic which is drawn in pieces:


  * Between angles $0$ and $\pi/2$ the horse has unconstrained access, so they can graze a wedge of radius $R$.
  * Between angles $\pi/2$ and until the horse's $y$ position is $0$ when the tether is taut the boundary of what can be eaten is described by the involute.
  * The horse can't eat from within the circle or radius $r$.


```{julia}
#| echo: false
let
    r,R = 1, 10
    R = max(R, pi*r) # R ≥ 1/2 circumference

    γ(θ) = -2r*cos(θ) * [cos(θ), sin(θ)] # parameterize the circle of radius r
    involute(t) = γ(t) + γ'(t)/norm(γ'(t))* (R - quadgk(u -> norm(γ'(u)), pi/2, t)[1])
    t₀ = find_zero(t -> round(involute(t)[2], digits=4), (3pi/4, pi))

    p = plot(; legend=false)
    plot_polar!(0..(pi/2), t -> R) # unobstructed -> quarter circle
    plot_parametric!((pi/2)..t₀, involute)
    plot_parametric!((pi/2)..pi, γ)
    plot!([0,R],[0,0])
end
```

To solve for the area we parameterize the circle of radius $r$ between $\pi/2$ and when the involute would cross the $x$ axis. We use `find_zero` to identify the value.


```{julia}
let
r,R = 160/(2π), 160
R = max(R, pi*r) # R ≥ 1/2 circumference
γ(θ) = -2r*cos(θ) * [cos(θ), sin(θ)]
## find *involute*: r - r'/|r'| * int(|r'|, a, t)
involute(t) = γ(t) + γ'(t)/norm(γ'(t))* (R - quadgk(u -> norm(γ'(u)), pi/2, t)[1])

t₀ = find_zero(t -> round(involute(t)[2], digits=4), (3pi/4, pi))

A₁ = π * R^2 / 4
y(t) = involute(t)[2]
x′(t) = (h=1e-4; (involute(t+h)[1]-involute(t-h)[1])/(2h))
A₂ = quadgk(t -> -y(t)*x′(t), pi/2, t₀)[1]  # A₂ = -∫ y dx, as counterclockwise parameterization
A₃ = (1/2) * π * r^2
2 * (A₁ + A₂ - A₃)
end
```

The calculation for $A_1$ and $A_3$ are from the familiar formula for the area of a circle. However, $A_2$ requires the formula for area above the $x$ axis when the curve is parameterized: $A = -\int_a^b y(t) x'(t) dt$, given how the curve is parameterized. As written, the automatic derivative of the numeric integral gives an error, so a central-difference approximation is used for $x'(t)$.


## Questions


###### Question


A cycloid is formed by pushing a wheel on a surface without slipping. The position of a fixed point on the outer rim of the wheel traces out the cycloid. Suppose the wheel has radius $R$ and the initial position of the point is at the bottom, $(0,0)$. Let $t$ measure angle measurement, in radians. Then the point of contact of the wheel will be at $Rt$, as that is the distance the wheel will have rotated. That is, the hub of the wheel will move according to $\langle Rt,~ R\rangle$. Relative to the hub, the point on the rim will have coordinates $\langle -R\sin(t), -R\cos(t) \rangle$, so the superposition gives:


$$
\vec{r}(t) = \langle Rt - R\sin(t), R - R\cos(t) \rangle.
$$

What is the position at $t=\pi/4$?


```{julia}
#| hold: true
#| echo: false
choices = [
q"[0.0782914, 0.292893 ]",
q"[0.181172, 0.5]",
q"[0.570796, 1.0]"]
answ = 1
radioq(choices, answ)
```

And the position at $\pi/2$?


```{julia}
#| hold: true
#| echo: false
choices = [
q"[0.0782914, 0.292893 ]",
q"[0.181172, 0.5]",
q"[0.570796, 1.0]"]
answ = 3
radioq(choices, answ)
```

###### Question


Suppose instead of keeping track of a point on the outer rim of the wheel, a point a distance $r < R$ from the hub is chosen in the above description of a cycloid (a [Curtate](http://mathworld.wolfram.com/CurtateCycloid.html) cycloid). If we start at $\langle 0,~ R-r \rangle$, what will be the position at $t$?


```{julia}
#| hold: true
#| echo: false
choices = [
" ``\\langle Rt - r\\sin(t),~ R - r\\cos(t) \\rangle``",
" ``\\langle Rt - R\\sin(t),~ R - R\\cos(t) \\rangle``",
" ``\\langle -r\\sin(t),~ -r\\cos(t) \\rangle``"
]
answ = 1
radioq(choices, answ)
```

###### Question


For the  cycloid $\vec{r}(t) = \langle t - \sin(t),~ 1 - \cos(t) \rangle$, find a simplified expression for $\| \vec{r}'(t)\|$.


```{julia}
#| hold: true
#| echo: false
choices = [
    " ``\\sqrt{2 - 2\\cos(t)}``",
    " ``1``",
    " ``1 - \\cos(t)``",
    " ``1 + \\cos(t) + \\cos(2t)``"
]
answ = 1
radioq(choices, answ)
```

###### Question


The cycloid $\vec{r}(t) = \langle t - \sin(t),~ 1 - \cos(t) \rangle$ has a formula for the arc length from $0$ to $t$ given by: $l(t) = 4 - 4\cos(t/2)$.


Plot the following two equations over $[0,8]$ which are a reparameterization of the cycloid by $l^{-1}(t)$.


```{julia}
#| hold: true
γ(s) = 2 * acos(1-s/4)
x1(s) = γ(s) - sin(γ(s))
y1(s) = 1 - cos(γ(s))
```

How many arches of the cycloid are traced out?


```{julia}
#| hold: true
#| echo: false
radioq(1:3, 1, keep_order=true)
```

###### Question


Consider the cycloid  $\vec{r}(t) = \langle t - \sin(t),~ 1 - \cos(t) \rangle$


What is the derivative at $t=\pi/2$?


```{julia}
#| hold: true
#| echo: false
choices = [
q"[1,1]",
q"[2,0]",
q"[0,0]"
]
answ = 1
radioq(choices, answ)
```

What is the derivative at $t=\pi$?


```{julia}
#| hold: true
#| echo: false
choices = [
q"[1,1]",
q"[2,0]",
q"[0,0]"
]
answ = 2
radioq(choices, answ)
```

###### Question


Consider the circle $\vec{r}(t) = R \langle \cos(t),~ \sin(t) \rangle$, $R > 0$. Find the norm of $\vec{r}'(t)$:


```{julia}
#| hold: true
#| echo: false
choices = [
    " ``1``",
    " ``1/R``",
    " ``R``",
    " ``R^2``"
]
answ = 3
radioq(choices, answ, keep_order=true)
```

###### Question


The curve described by $\vec{r}(t) = \langle 10t,~ 10t - 16t^2\rangle$ models the flight of an arrow. Compute the length traveled from when it is launched to when it returns to the ground.


```{julia}
#| hold: true
#| echo: false
x(t) = 10t
y(t) = 10t - 16t^2
a,b = sort(find_zeros(y, -10, 10))
f(x,y) = 1
val, _ = quadgk(t -> f(x(t), y(t)) * sqrt(D(x)(t)^2 + D(y)(t)^2), a, b)
numericq(val)
```

###### Question


Let $\vec{r}(t) = \langle t, t^2 \rangle$ describe a parabola. What is the arc length between $0 \leq t \leq 1$? First, what is a formula for the speed ($\| \vec{r}'(t)\|$)?


```{julia}
#| hold: true
#| echo: false
choices = [
    " ``\\sqrt{1 + 4t^2}``",
    " ``1 + 4t^2``",
    " ``1``",
    " ``t + t^2``"
]
answ = 1
radioq(choices, answ)
```

Numerically find the arc length.


```{julia}
#| hold: true
#| echo: false
val,err = quadgk(t -> (1 + 4t^2)^(1/2), 0, 1)
numericq(val)
```

###### Question


Let $\vec{r}(t) = \langle t, t^2 \rangle$ describe a parabola. What is the curvature of $\vec{r}(t)$ at $t=0$?


```{julia}
#| hold: true
#| echo: false
@syms t::positive
rt = [t, t^2, 0]
rp = diff.(rt, t)
rpp = diff.(rt, t, t)
kappa =  norm(rp × rpp) / norm(rp)^3
#val = N(kappa(t=>0)) #2
val = 2
numericq(val)
```

The curvature at $1$ will be


```{julia}
#| hold: true
#| echo: false
choices = [
"greater than the curvature at ``t=0``",
"less than the curvature at ``t=0``",
"the same as the curvature at ``t=0``"]
answ = 2
radioq(choices, answ)
```

The curvature as $t\rightarrow \infty$ will be


```{julia}
#| hold: true
#| echo: false
choices = [
    " ``0``",
    " ``\\infty``",
    " ``1``"
]
answ = 1
radioq(choices, answ)
```

---


Now, if we have a more general parabola by introducing a parameter $a>0$: $\vec{r}(t) = \langle t, a\cdot t^2 \rangle$, What is the curvature of $\vec{r}(t)$ at $t=0$?


```{julia}
#| hold: true
#| echo: false
choices = [
    " ``2a``",
    " ``2/a``",
    " ``2``",
    " ``1``"
]
answ = 1
radioq(choices, answ)
```

###### Question


Projectile motion with constant acceleration is expressed parametrically by $\vec{x}(t) = \vec{x}_0 + \vec{v}_0 t + (1/2) \vec{a} t^2$, where $\vec{x}_0$ and $\vec{v}_0$ are initial positions and velocity respectively. In [Strang](https://ocw.mit.edu/resources/res-18-001-calculus-online-textbook-spring-2005/textbook/MITRES_18_001_strang_12.pdf) p451, we find an example utilizing this formula to study the curve of a baseball. Place the pitcher at the origin, the batter along the $x$ axis, then a baseball thrown with spin around its $z$ axis will have acceleration in the $y$ direction in addition to the acceration due to gravity in the $z$ direction. Suppose the ball starts $5$ feet above the ground when pitched ($\vec{x}_0 = \langle 0,0, 5\rangle$), and has initial velocity $\vec{v}_0 = \langle 120, -2, 2 \rangle$. ($120$ feet per second is about $80$ miles per hour). Suppose the pitcher can produce an acceleration in the $y$ direction of $16ft/sec^2$, then $\vec{a} = \langle 0, 16, -32\rangle$ in these units. (Gravity is $9.8m/s^2$ or $32ft/s^2$.)


The plate is $60$ feet away. How long will it take for the ball to reach the batter? (When the first component is $60$?)


```{julia}
#| hold: true
#| echo: false
x0 = [0,0,5]
v0 = [120, -2, 2]
a = [0, 16, -32]
r(t) = x0 + v0*t + 1/2*a*t^2
answ = 60/v0[1]
numericq(answ)
```

At $t=1/4$ the ball is half-way to home. If the batter reads the ball at this point, where in the $y$ direction is the ball?


```{julia}
#| hold: true
#| echo: false
x0 = [0,0,5]
v0 = [120, -2, 2]
a = [0, 16, -32]
r(t) = x0 + v0*t + 1/2*a*t^2
t = 1/4
answ = r(t)[2]
numericq(answ)
```

At $t=1/2$ has the ball moved more than $1/2$ foot in the $y$ direction?


```{julia}
#| hold: true
#| echo: false
x0 = [0,0,5]
v0 = [120, -2, 2]
a = [0, 16, -32]
r(t) = x0 + v0*t + 1/2*a*t^2
t = 1/2
answ = abs(r(t)[2]) > 1/2
yesnoq(answ)
```

###### Question


In [Strang](https://ocw.mit.edu/resources/res-18-001-calculus-online-textbook-spring-2005/textbook/MITRES_18_001_strang_12.pdf) we see this picture describing  a curve:


```{julia}
#| hold: true
#| echo: false
a = 1

plot(t -> 0, -2, 2,aspect_ratio=:equal, legend=false)
plot!(t -> 2a)
r(t) = [0, a] + a*[cos(t), sin(t)]
plot!(unzip(r, 0, 2pi)...)
theta = pi/3
plot!([0, 2a/tan(theta)], [0, 2a], linestyle=:dash)
A = [2a*cot(theta), 2a]
B = 2a*sin(theta)^2 *[ 1/tan(theta),1]
scatter!(unzip([A,B])...)
plot!([B[1],A[1],A[1]], [B[2],B[2],A[2]], linestyle=:dash)
delta = 0.2
annotate!([(B[1],B[2]-delta,"B"),(A[1]+delta,A[2]-delta,"A")])
r(theta) = [2a*cot(theta), 2a*sin(theta)^2 ]
theta0 = pi/4
plot!(unzip(r, theta0, pi-theta0)..., linewidth=3)
P = r(theta)
annotate!([(P[1],P[2]-delta, "P")])
```

Strang notes that the curve is called the "witch of Agnesi" after Maria Agnesi, the author of the first three-semester calculus book. (L'Hopital's book did not contain integration.)


We wish to identify the parameterization. Using $\theta$ an angle in standard position, we can see that the component functions $x(\theta)$ and $y(\theta)$ may be found using trigonometric analysis.


What is the $x$ coordinate of point $A$? (Also the $x$ coordinate of $P$.)


```{julia}
#| hold: true
#| echo: false
choices = [
           " ``2\\cot(\\theta)``",
           " ``\\cot(\\theta)``",
           " ``2\\tan(\\theta)``",
           " ``\\tan(\\theta)``"
       ]
answ = 1
radioq(choices, answ)
```

Using the polar form of a circle, the length between the origin and $B$ is given by $2\cos(\theta-\pi/2) = 2\sin(\theta)$. Using this, what is the $y$ coordinate of $B$?


```{julia}
#| hold: true
#| echo: false
choices = [
    " ``2\\sin^2(\\theta)``",
    " ``2\\sin(\\theta)``",
    " ``2``",
    " ``\\sin(\\theta)``"
]
answ=1
radioq(choices, answ)
```

###### Question


Let $n > 0$, $\vec{r}(t) = \langle t^{(n+1)},t^n\rangle$. Find the speed, $\|\vec{r}'(t)\|$.


```{julia}
#| hold: true
#| echo: false
choices = [
    " ``\\frac{\\sqrt{n^{2} t^{2 n} + t^{2 n + 2} \\left(n + 1\\right)^{2}}}{t}``",
    " ``t^n + t^{n+1}``",
    " ``\\sqrt{n^2 + t^2}``"
]
answ=1
radioq(choices, answ)
```

For $n=2$, the arc length of $\vec{r}$ can be found exactly. What is the arc-length between $0 \leq t \leq a$?


```{julia}
#| hold: true
#| echo: false
choices = [
    " ``\\frac{a^{2} \\sqrt{9 a^{2} + 4}}{3} + \\frac{4 \\sqrt{9 a^{2} + 4}}{27} - \\frac{8}{27}``",
    " ``\\frac{2 a^{\\frac{5}{2}}}{5}``",
    " ``\\sqrt{a^2 + 4}``"
]
answ = 1
radioq(choices, answ)
```

###### Question


The [astroid](http://www-history.mcs.st-and.ac.uk/Curves/Astroid.html) is one of the few curves with an exactly computable arc-length. The curve is parametrized by $\vec{r}(t) = \langle a\cos^3(t), a\sin^3(t)\rangle$. For $a=1$ find the arc-length between $0 \leq t \leq \pi/2$.


```{julia}
#| hold: true
#| echo: false
choices = [
    " ``\\sqrt{2}``",
    " ``3/2``",
    " ``\\pi/2``",
    " ``2``"
]
answ = 2
radioq(choices, answ, keep_order=true)
```

###### Question


```{julia}
#| echo: false
let
    t0, t1 = pi/12, pi/3
    tspan = (t0, t1)  # time span to consider

    a = 1
    r(theta) = -cos(theta) + 4*2cos(theta)*sin(theta)^2
    F(t) = r(t) * [cos(t), sin(t)]
    p = (a, F)      # combine parameters

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

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

    plt = plot(unzip(F, t0, t1)..., legend=false, color=:red)
    plot!(plt, unzip(t->out(t),  t0, t1)..., color=:blue)
end
```

Let $F$ and $B$ be pictured above. Which is the red curve?


```{julia}
#| hold: true
#| echo: false
choices = [
"The front wheel",
"The back wheel"
]
answ=1
radioq(choices, answ)
```

###### Question


```{julia}
#| echo: false
let
    t0, t1 = 0.0, pi/3
    tspan = (t0, t1)  # time span to consider

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

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

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

    plt = plot(unzip(F, t0, t1)..., legend=false, color=:blue)
    plot!(plt, unzip(t->out(t),  t0, t1)..., color=:red)
end
```

Let $F$ and $B$ be pictured above. Which is the red curve?


```{julia}
#| hold: true
#| echo: false
choices = [
"The front wheel",
"The back wheel"
]
answ=2
radioq(choices, answ)
```

###### Question


Let $\vec{\gamma}(s)$ be a parameterization of a curve by arc length and $s(t)$ some continuous increasing function of $t$. Then $\vec{\gamma} \circ s$ also parameterizes the curve. We have


$$
\text{velocity}  = \frac{d (\vec{\gamma} \circ s)}{dt} = \frac{d\vec{\gamma}}{ds} \frac{ds}{dt} = \hat{T} \frac{ds}{dt}.
$$

Continuing with a second derivative


$$
\text{acceleration} = \frac{d^2(\vec{\gamma}\circ s)}{dt^2} =
\frac{d\hat{T}}{ds} \frac{ds}{dt} \frac{ds}{dt} + \hat{T} \frac{d^2s}{dt^2} = \frac{d^2s}{dt^2}\hat{T} + \kappa (\frac{ds}{dt})^2 \hat{N},
$$

Using $d\hat{T}/{ds} = \kappa\hat{N}$ when parameterized by arc length.


This expresses the acceleration in terms of the tangential part and the normal part. [Strang](https://ocw.mit.edu/resources/res-18-001-calculus-online-textbook-spring-2005/textbook/MITRES_18_001_strang_12.pdf) views this in terms of driving where the car motion is determined by the gas pedal and the brake pedal only giving acceleration in the $\hat{T}$ direction) and the steering wheel (giving acceleration in the $\hat{N}$ direction).


If a car is on a straight road, then $\kappa=0$. Is the acceleration along the $\hat{T}$ direction or the $\hat{N}$ direction?


```{julia}
#| hold: true
#| echo: false
choices = [
    "The ``\\hat{T}`` direction",
    "The ``\\hat{N}`` direction"]
answ = 1
radioq(choices, answ)
```

Suppose no gas or brake is applied for a duration of time. The tangential acceleration will be $0$. During this time, which of these must  be $0$?


```{julia}
#| hold: true
#| echo: false
choices = [
    " ``\\vec{\\gamma} \\circ s``",
    " ``ds/dt``",
    " ``d^2s/dt^2``"
]
answ = 3
radioq(choices, answ, keep_order=true)
```

In going around a corner (with non-zero curvature), which is true?


```{julia}
#| hold: true
#| echo: false
choices = [
"The acceleration in the normal direction depends on both the curvature and the speed (``ds/dt``)",
"The acceleration in the normal direction depends only on the curvature and not the speed (``ds/dt``)",
"The acceleration in the normal direction depends only on the speed (``ds/dt``) and not the curvature"
]
answ = 1
radioq(choices, answ)
```

###### Question


The evolute comes from the formula $\vec\gamma(T) - (1/\kappa(t)) \hat{N}(t)$. For hand computation, this formula can be explicitly given by two components $\langle X(t), Y(t) \rangle$ through:



\begin{align*}
r(t) &= x'(t)^2 + y'(t)^2\\
k(t) &= x'(t)y''(t) - x''(t) y'(t)\\
X(t) &= x(t) - y'(t) r(t)/k(t)\\
Y(t) &= y(t) + x'(t) r(t)/k(t)
\end{align*}


Let $\vec\gamma(t) = \langle t, t^2 \rangle = \langle x(t), y(t)\rangle$ be a parameterization of a parabola.


  * Compute $r(t)$


```{julia}
#| hold: true
#| echo: false
choices = [
    " ``1 + 4t^2``",
    " ``1 - 4t^2``",
    " ``1 + 2t``",
    " ``1 - 2t``"
]
answ = 1
radioq(choices, answ)
```

  * Compute $k(t)$


```{julia}
#| hold: true
#| echo: false
choices = [
    " ``2``",
    " ``-2``",
    " ``8t``",
    " ``-8t``"
]
answ = 1
radioq(choices, answ)
```

  * Compute $X(t)$


```{julia}
#| hold: true
#| echo: false
choices = [
    " ``t - 2t(1 + 4t^2)/2``",
    " ``t - 4t(1+2t)/2``",
    " ``t - 2(8t)/(1-2t)``",
    " ``t - 1(1+4t^2)/2``"
]
answ = 1
radioq(choices, answ)
```

  * Compute $Y(t)$


```{julia}
#| hold: true
#| echo: false
choices = [
    " ``t^2 + 1(1 + 4t^2)/2``",
    " ``t^2 + 2t(1+4t^2)/2``",
    " ``t^2 - 1(1+4t^2)/2``",
    " ``t^2 - 2t(1+4t^2)/2``"
]
answ = 1
radioq(choices, answ)
```

###### Question


The following will compute the evolute of an ellipse:


```{julia}
#| eval: false
@syms t a b
x = a * cos(t)
y = b * sin(t)
xp, xpp, yp, ypp = diff(x, t), diff(x,t,t), diff(y,t), diff(y,t,t)
r2 = xp^2 + yp^2
k = xp * ypp - xpp * yp
X = x - yp * r2 / k |> simplify
Y = y + xp * r2 / k |> simplify
[X, Y]
```

What is the resulting curve?


```{julia}
#| hold: true
#| echo: false
choices = [
"An astroid of the form ``c \\langle \\cos^3(t), \\sin^3(t) \\rangle``",
"An cubic parabola of the form ``\\langle ct^3, dt^2\\rangle``",
"An ellipse of the form ``\\langle a\\cos(t), b\\sin(t)``",
"A cyloid of the form ``c\\langle t + \\sin(t), 1 - \\cos(t)\\rangle``"
]
answ = 1
radioq(choices, answ)
```
