daviddoji/CalculusWithJuliaNotes.jl
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

sequences and series

Browse Source
This commit is contained in:
jverzani
2025-08-14 19:08:41 -04:00
parent a192358e96
commit 042a332f71
8 changed files with 1578 additions and 97 deletions
Download Patch File Download Diff File Expand all files Collapse all files

197
quarto/differentiable_vector_calculus/vector_valued_functions.qmd
View File

@@ -38,9 +38,11 @@ A function $\vec{f}: R \rightarrow R^n$, $n > 1$ is called a vector-valued funct
$$
\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.
\begin{align*}
\vec{f}(t) &= \langle \sin(t), 2\cos(t) \rangle, \\
\vec{g}(t) &= \langle \sin(t), \cos(t), t \rangle, \\
\vec{h}(t) &= \langle 2, 3 \rangle + t \cdot \langle 1, 2 \rangle.\\
\end{align*}
$$
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.
@@ -107,21 +109,22 @@ In `Plots`, the command `plot(xs, ys)`, where, say, `xs=[x1, x2, ..., xn]` and `
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.)
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 vectors. In `Plots`, a vector of *tuples* will be read as a vector of points and plotted accordingly. On the other hand, a vector of vectors is read in as a number of series, with each element being plotted separately. (That is, `[x1,y1]` maps to `plot!([1,2], [x1,y1])`.) To get the desired graph, *either* our function can return a tuple---which makes it clumsier to work with when manipulating the output---or we can turn a vector of points into two vectors---one with the `x` values, one with the `y` values.
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]`).
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:
Visually, we have `unzip` performing this re-association:
```{verbatim}
[[x1, y1, z1], (⌈x1⌉, ⌈y1⌉, ⌈z1⌉,
[x2, y2, z2], |x2|, |y2|, |z2|,
[x3, y3, z3], --> |x3|, |y3|, |z3|,
[[x₁, y₁, z₁], (⌈x₁⌉, ⌈y₁⌉, ⌈z₁⌉,
[x₂, y₂, z₂], |x₂|, |y₂|, |z₂|,
[x₃, y₃, z₃], --> |x₃|, |y₃|, |z₃|,
⋮ ⋮
[xn, yn, zn]] ⌊xn⌋, ⌊yn⌋, ⌊zn⌋ )
[xₙ, yₙ, zₙ]] ⌊xₙ⌋, ⌊yₙ⌋, ⌊zₙ⌋ )
```
To turn a collection of vectors into separate arguments for a function, splatting (the `...`) is used.
@@ -175,6 +178,18 @@ ts = range(-2, 2, length=200)
plot(unzip(h.(ts))...)
```
### Using points, not vectors
As mentioned, there is an alternate manner to plot a vector-valued function that has some conveniences. This is to use a tuple to store the component values. For example:
```{julia}
g(t) = (cos(t) + 1/5 * cos(5t), sin(t) + 2/3*sin(3t))
ts = range(0, 2pi, 251)
plot(g.(ts); legend=false, aspect_ratio=:equal)
```
Broadcasting `g` creates a vector of tuples, which `Plots` treats as points. The drawback to this approach, as mentioned, is that manipulating the output is generally easily when the function output is a vector.
### The `plot_parametric` function
@@ -229,8 +244,8 @@ For example:
#| 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)
plot_parametric(0..2pi, f; legend=false)
scatter!([(0,0)]; markersize=4)
```
@@ -272,14 +287,14 @@ 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)
p = plot(; legend=false, aspect_ratio=:equal)
circle!([0,0], R; linecolor=:blue)
circle!(cent(t), r; linecolor=:black)
tp(t) = -R/r * t
tp(t) = -R / r * t
s(t) = cent(t) + rho * [cos(tp(t)), sin(tp(t))]
plot_parametric!(0..t, s, color=:red)
plot_parametric!(0..t, s; linecolor=:red)
p
end
@@ -479,10 +494,15 @@ 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)
plot(; legend=false, aspect_ratio=:equal)
plot_parametric!(pi/5..3pi/8, f; line=(:red, 1))
arrow!([0,0], f(t); line=(:blue,))
arrow!([0,0], f(t + Δt); line=(:blue,))
arrow!(f(t), df; line=(:black, 3,0.5))
annotate!([(f(t)..., text("f(t)", :bottom, :left)),
(f(t+Δt)..., text("f(t + Δt)", :bottom, :left)),
((f(t) + df/2)..., text("df", :top, :right)),
])
```
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:
@@ -514,9 +534,12 @@ 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)
plot(; legend=false, aspect_ratio=:equal)
p = plot_parametric!(0..2pi, f; line=(:black, 1))
for t in [1,2,3]
arrow!(f(t), f'(t)) # add arrow with tail on curve, in direction of derivative
arrow!([0,0], f(t); line=(:gray, 1, 0.5))
annotate!((2f(t)/3)..., text("f($t)", :top, :left))
arrow!(f(t), f'(t); line=(:blue, 2)) # add arrow with tail on curve, in direction of derivative
end
p
```
@@ -528,8 +551,8 @@ Were symbolic expressions used in place of functions, the vector-valued function
```{julia}
@syms 𝒕
𝒗vf = [cos(𝒕), sin(𝒕), 𝒕]
@syms t
vvf = [cos(t), sin(t), t]
```
We will see working with these expressions is not identical to working with a vector-valued function.
@@ -539,7 +562,7 @@ To plot, we can avail ourselves of the the parametric plot syntax. The following
```{julia}
plot(𝒗vf..., 0, 2pi)
plot(vvf..., 0, 2pi)
```
The `unzip` usage, as was done above, could be used, but it would be more trouble in this case.
@@ -549,7 +572,7 @@ To evaluate the function at a given value, say $t=2$, we can use `subs` with bro
```{julia}
subs.(𝒗vf, 𝒕=>2)
subs.(vvf, t=>2)
```
Limits are performed component by component, and can also be defined by broadcasting, again with the need to adjust the values:
@@ -557,21 +580,21 @@ Limits are performed component by component, and can also be defined by broadcas
```{julia}
@syms Δ
limit.((subs.(𝒗vf, 𝒕 => 𝒕 + Δ) - 𝒗vf) / Δ, Δ => 0)
limit.((subs.(vvf, t => t + Δ) - vvf) / Δ, Δ => 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, 𝒕)
diff.(vvf, t)
```
The second derivative, can be found through:
```{julia}
diff.(𝒗vf, 𝒕, 𝒕)
diff.(vvf, t, t) # or diff.(vvf, t, 2)
```
### Applications of the derivative
@@ -586,13 +609,13 @@ Here are some sample applications of the derivative.
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.
With this, the equation is simply $\vec{tl}(t) = \vec{f}(t_0) + \cdot (t - t_0) \vec{f}'(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$.
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(t) = 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:
@@ -606,7 +629,8 @@ 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)
plot(; legend=false)
plot_parametric!(0..t_0, xpos)
```
```{julia}
@@ -797,7 +821,7 @@ The dot being scalar multiplication by the derivative of the univariate function
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:
For the dot product, the combination $\vec{f}(t) \cdot \vec{g}(t)$ creates 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:
$$
@@ -842,8 +866,8 @@ In summary, these two derivative formulas hold for vector-valued functions $R \r
$$
\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}'.
(\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}'\qaud (n=3).
\end{align*}
$$
@@ -892,7 +916,9 @@ $$
\vec{F} = m \vec{a} = m \ddot{\vec{x}}.
$$
Combining, Newton states $\vec{a} = -(GM/r^2) \hat{x}$.
(The double dot is notation for two derivatives in a $t$ variable.)
Combining, Newton's law states $\vec{a} = -(GM/r^2) \hat{x}$.
Now to show the first law. Consider $\vec{x} \times \vec{v}$. It is constant, as:
@@ -901,7 +927,8 @@ 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}.
&= \vec{v} \times \vec{v} + \vec{x} \times \vec{a}\\
&= 0.
\end{align*}
$$
@@ -1000,7 +1027,7 @@ c^2 &= \|\vec{c}\|^2 \\
$$
Solving, this gives the first law. That is, the radial distance is in the form of an ellipse:
Solving for $r$, this gives the first law. That is, the radial distance is in the form of an ellipse:
$$
@@ -1231,23 +1258,25 @@ p
---
::: {.callout-note appearance="minimal"}
## Curvature of a space curve
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}.$
$$
\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
For $2$-dimensional space curves, the same formula applies after embedding a $0$ third component. It simplifies to:
$$
\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 $\sin(\theta)$ in the definition). For "curvy" curves, $\vec{r}''$ will be in a direction orthogonal of $\vec{r}'$ to the $\sin(\theta)$ term in the cross product will be closer to $1$.
:::
Curvature can also be defined by the 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 $\sin(\theta)$ in the definition). For "curvy" curves, $\vec{r}''$ will be in a direction orthogonal of $\vec{r}'$ to the $\sin(\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:
@@ -1269,17 +1298,16 @@ If a curve is imagined to have a tangent "circle" (second order Taylor series ap
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.
::: {.callout-note appearance="minimal}
## Torsion of a space curve
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}.$
$$
\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
@@ -1294,7 +1322,7 @@ This last example comes from a collection of several [examples](https://github.c
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:
The Frenet-Serret equations describe the relationship between tangent, normal, and binormal vectors through a system of differential equations driven by the curvature and torsion. Their derivation will be discussed later, here we give an example usage by using the `DifferentialEquations` package to solve these equations for a specific curvature and torsion function and given initial conditions. The vector equations along with the relationship between the space curve and its tangent vector---when expanded into coordinates---become $12$ different equations:
```{julia}
@@ -1316,7 +1344,7 @@ function Frenet_eq!(du, u, p, s) #system of ODEs
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$.
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 $\lVert c'(t) \rVert = 1$).
Following the example of `@empet`, we define a curvature function and torsion function, the latter a constant:
@@ -1345,7 +1373,7 @@ 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:
The "spine" is the center axis of the tube and is described by the $10$th, $11$th, and $12$th coordinates:
```{julia}
@@ -1372,14 +1400,15 @@ 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)
plot(; legend=false)
plot_parametric!(a_0..b_0, spine; line=(:black, 2))
arrow!(spine(t_0), e₁(t_0))
arrow!(spine(t_0), e₂(t_0))
arrow!(spine(t_0), e₃(t_0))
arrow!(spine(t_0), e₁(t_0); line=(:blue,))
arrow!(spine(t_0), e₂(t_0); line=(:red,))
arrow!(spine(t_0), e₃(t_0); line=(:green,))
r_0(t, θ) = spine(t) + ϵ * (e₂(t)*cos(θ) + e₃(t)*sin(θ))
plot_parametric!(0..2pi, θ -> r_0(t_0, θ))
plot_parametric!(0..2pi, θ -> r_0(t_0, θ); line=(:black, 1))
```
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.
@@ -1475,7 +1504,7 @@ 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$:
We see that the answer depends linearly on $a$, but otherwise is a constant involving an integral. We use `QuadGk` to provide a numeric answer for the case $a=1$:
```{julia}
@@ -1573,7 +1602,7 @@ 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.
@@ -1778,13 +1807,13 @@ 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)
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ₑ)
plot_parametric!(-4..4, x -> Pₑ + x*Vₑ; line=(:black, 1))
end
plot!(Xₑ, Yₑ, 0, 2pi, linewidth=5)
plot!(Xₑ, Yₑ, 0, 2pi, line=(:red, 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.
@@ -1858,8 +1887,10 @@ 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])
plot(; legend=false, aspect_ratio=:equal, xlims=(-6,6), ylims=(-5,5))
plot_parametric!(0..2pi, t -> rₑ₃(t)[1:2]; line=(:red,5))
plot_parametric!(0..2pi, t -> (rₑ₃(t) + Normal(rₑ₃, t)/curvature(rₑ₃, t))[1:2];
line=(:black, 1))
```
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)'$:
@@ -1877,6 +1908,7 @@ $$
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:
@@ -1891,10 +1923,10 @@ end
#| 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])
plot(; legend=false, aspect_ratio=:equal)
plot_parametric!(0..2pi, t -> r(t)[1:2]; line=(:black, 1))
plot_parametric!(0..2pi, t -> evolute(r)(t)[1:2]; line=(:red, 1))
plot_parametric!(0..2pi, t -> ((evolute∘evolute)(r)(t))[1:2]; line=(:blue,1))
```
---
@@ -1911,9 +1943,10 @@ 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)
p = plot(; legend=false)
plot_parametric!(-2..2, r; line=(:black, 1))
plot_parametric!(t0..t1, t -> beta(r, t); line=(:red,1))
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
@@ -1924,14 +1957,15 @@ This lends itself to this mathematical description, if $\vec\gamma(t)$ parameter
$$
\vec\beta(t) = \vec\gamma(t) + \left((a - \int_{t_0}^t \| \vec\gamma'(t)\| dt) \hat{T}(t)\right),
\vec\beta(t) = \vec\gamma(t) + \left(a - \int_{t_0}^t \| \vec\gamma'(t)\| dt\right) \hat{T}(t),
$$
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.
\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)$:
@@ -1989,7 +2023,7 @@ $$
$$
That is the evolute of an involute of $\vec\gamma(s)$ is $\vec\gamma(s)$.
That is, the evolute of an involute of $\vec\gamma(s)$ is $\vec\gamma(s)$.
We have:
@@ -2041,8 +2075,9 @@ 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)))
plot(; legend=false)
plot_parametric!(0..4pi, r; line=(:black, 1))
plot_parametric!(0..4pi, u -> float(subs.(ex, t => u)); line=(:blue, 1))
```
The expression `ex` is secretly `[t + sin(t), 3 + cos(t)]`, another cycloid.