In these notes, we will need these packages. The first is in the standard library, the others must previously have been installed.
using LinearAlgebra # for norm, dot, and cross
using Plots # for plotting
using SymPy # for symbolic math
using Roots # to find zeros numerically
using QuadGK # to integerate numerically
(If the MTH229 package is installed, the last four lines could simply be
using MTH229`.)
In julia
, creating vectors is very straightforward: in place of the "book" notation $\langle 1,2,3 \rangle$, we use square brackets, as in
[1,2,3]
3-element Array{Int64,1}: 1 2 3
In Julia, this creates a vector. The output describes itself as Array{Int64,1}
. This translates into an array of integers with one dimension. An array is a collection of objects stored in a multi-dimensional grid. For vectors the dimension is 1, for matrices – two-dimensional, rectangular arrays – the dimension is 2. This is important, note the difference if we leave the commas out of the above:
[1 2 3]
1×3 Array{Int64,2}: 1 2 3
This creates a 1x3 Array{Int64,2}
, which has 1 row, 3 columns, 2 dimensions and stores integers. The only moral here – use commas to create vectors.
The
[]
notation is used in many different ways injulia
. In the above we see its use to combine like-type items into a collection, in this case a vector or matrix. As well, the[]
notation is used to access vector components, for examplex[2]
would be the second component of the collectionx
. It is also common with list comprehensions of the form[x^2 for x in 1:5]
, though this is related to the first use. More generally, identifying the different uses of matched braces:[]
,()
, and{}
is important in trying to readjulia
syntax.
Vectors have many arithmetic operations defined for them. These fit naturally into julia
's syntax.
u = [1,2,3]
v = [4,5,6]
w = [1,3,5]
2*u, u + v, v - w
([2, 4, 6], [5, 7, 9], [3, 2, 1])
The number of entries in a vector is determined by length
:
length(u)
3
The norm of the vector is returned by norm
, but first you need to load the LinearAlgebra
package, as was done above. The default is to use the Euclidean norm (p=2
), though others are possible:
norm(u)
3.7416573867739413
This makes finding a unit vector trivial:
u/norm(u)
3-element Array{Float64,1}: 0.2672612419124244 0.5345224838248488 0.8017837257372732
We make a function to compute the unit vector of any vector:
uvec(u) = u/norm(u)
uvec (generic function with 1 method)
If we give our variable names u\hat<tab>
then a "hat" will match the book for a unit vector:
û = uvec(u)
3-element Array{Float64,1}: 0.2672612419124244 0.5345224838248488 0.8017837257372732
(There is nothing but a naming convention. Using a hat does not by itself make something a unit vector.)
Here we verify that the triangle inequality holds for u
and v
:
norm(u+v) <= norm(u) + norm(v)
true
The dot product is defined for any two vectors of the same length and is implemented via dot
:
dot(u, v)
32
For those who like to match the text book, you can use \cdot<tab>
as an infix operation:
u ⋅ v
32
using $\hat{u} \cdot \hat{v} = \cos(\theta)$ To find the angle between the two vectors u
and v
we have:
acos(uvec(u) ⋅ uvec(v))
0.2257261285527342
The projection of u
along v
is given by the projection formula $u \cdot \hat{v} \hat{v}$:
ev = uvec(v)
(u ⋅ ev) * ev # parentheses are not needed here, as dot happens before *
3-element Array{Float64,1}: 1.6623376623376622 2.0779220779220777 2.493506493506493
Lets break u
up into a parallel and perpendicular components in terms of v
:
ev = uvec(v)
upar = (u ⋅ ev) * ev
uperp = u - upar
3-element Array{Float64,1}: -0.6623376623376622 -0.07792207792207773 0.506493506493507
And we check that uperp
is indeed perpendicular, as it should have a 0 dot product with v
:
uperp ⋅ v # 0 up to roundoff errors
4.440892098500626e-15
For 3D vectors, like u
and v
, the cross product is defined and implemented in cross
:
cross(u, v)
3-element Array{Int64,1}: -3 6 -3
Again, an infix operator is available: \times<tab>
:
u × v
3-element Array{Int64,1}: -3 6 -3
Here we verify for these vectors that the cross product is not commutative, but rather anti-commutative:
(u × v) - (v × u), (u × v) + (v × u)
([-6, 12, -6], [0, 0, 0])
Here we verify that the cross product is not-associative
left = (u × v) × w # also cross(cross(u, v), w)
right = u × (v × w) # and cross(u, cross(v, w))
left - right
3-element Array{Int64,1}: -17 -2 13
And this finds, again, the angle between u
and v
:
û, v̂ = uvec(u), uvec(v) # u\hat<tab> and v\hat<tab> to enter
asin(norm( û × v̂ ))
0.22572612855273397
Let's verify that $\| u \times v \|^2 = \| u \|^2 \| v\|^2 - (u \cdot v)^2$:
left = norm( u × v )^2
right = norm(u)^2 * norm(v)^2 - (u ⋅ v)^2
left - right # up to round off error
-2.2737367544323206e-13
A general definition of a function is a mapping from a domain to a range. A vector-valued function is a function which takes a value from the real line and returns a vector. In shorthand, $f: R \rightarrow R^n$, where $n=2,3, ...$. In this section, we see how to define such functions, how to operate on them, and visualize them. In this specific case of a function from the real line to 2 or 3 dimensions, the term parameterized curve is used to describe the function.
In julia
two reasonable approaches to implementing a vector-valued function are:
a function that returns a vector, as in
vv(t) = [sin(t), cos(t), t]
vv (generic function with 1 method)
a vector of functions, as in
vf = [sin, cos, t -> t]
3-element Array{Function,1}: sin cos getfield(WeavePynb.ZKWQZSUDT, Symbol("##1#2"))()
The latter uses an anonymous function in its last component. One can convert between the two representations using a comprehension:
as_vf(vv,n=length(vv(0))) = [t -> vv(t)[i] for i in 1:n]
as_vv(vf) = t -> [f(t) for f in vf]
as_vv (generic function with 1 method)
(The as_vf
function needs to know the $n$ in $R^n$. It can be specifed as a second argument or, by default, is computed from the value of vv(0)
.)
The function that returns a vector is generally easier to work with, except for the case of plotting, so we will functions that return a vector (vv
).
As with most cases in math, we begin with lines. A line in 2 or 3 dimensions can be thought of in terms of scalar multiples of a vector from a fixed point. Let $p$ be the fixed point, and $v$ the vector, then we can represent any point on a line going through $p$ in the direction of $v$ by:
In julia
, we have, for example:
p = [1,2,3]
v = [3,2,1]
f(t) = p + t*v
f(0), f(1), f(2)
([1, 2, 3], [4, 4, 4], [7, 6, 5])
A typical helix is parameterized by
helix(t) = [sin(t), cos(t), t]
helix (generic function with 1 method)
Another example we will use, is this function, which when graphed will make a "pringle":
pringle(t) = [cos(t), sin(t), sin(2t)]
pringle (generic function with 1 method)
Two-dimensional vector-valued are produced in a similar manner with just two components being satisfied.
We have used the Plots
package for plotting, primarily as it has many conveniences. In particular, we can plot the function $f$ over $[a,b]$ through the command plot(f, a, b)
. This is a convenience for the task of:
generating $x$ values between $a$ and $b$;
computing the $y$ values, $f(x)$;
plotting the paired points $(x,y)$ and connecting with a line, as needed.
For plots in 3 dimensions, however, at times must follow steps like the above, rather than use a more simplified interface.
To review how to do the three steps above, we have
generating $x$ values between $a$ and $b$:
We have two basic ways to do this: Either through the notation xs=a:h:b
which steps from a
towards b
by steps of size h
. Or xs=range(a, stop=b, length=n)
which finds n
evenly spaced points from a
to b
.
computing the $y$ values, $f(x)$:
If we have a function already defined, then the "dot" syntax allows one to easily broadcast over the x
values, as in ys = f.(xs)
.
If we have a mathematical expression to apply to the x
values, we can use a comprenshion, as in ys=[sin(2x) for x in xs]
.
plotting the paired points $(x,y)$ and connecting with a line:
The call plot(xs,ys)
will generate a line plot. The scatter
function is used to plot just the points.
We make a plot of the sine curve over $[0,2\pi]$ and its tangent line at $c=\pi/4$, as follows:
First, we plot the sine curve:
f(x) = sin(x)
xs = range(0, stop=2pi, length=100) # 100 points from 0 to 2pi, evenly spaced
ys = sin.(xs)
plot(xs, ys)
We add to the graphic, using the related plot!
function. Here we add a tangent line at $c=\pi/4$:
c = pi/4; m = cos(c)
tl(x) = f(c) + m * (x-c)
tls = tl.(xs)
plot(xs, ys) # redo plot
plot!(xs, tls) # add tangent line
To make a parametric plot of $(f(x), g(x))$ is similar, though both the xs
and ys
are generated in terms of $t$ values. For example, here we plot $(\sin(t), \sin(2t))$ over $[0, 2pi]$:
ts = range(0, stop=2pi, length=100)
xs = sin.(ts)
ys = sin.(2ts)
plot(xs, ys)
(There is also the interface plot(f1, f2, a, b)
, but we will not use this.)
For a 3-dimensional parametric plot, the plot
function is also used, the steps are identical. The helix, a plot of $(\sin(t), \cos(t), t)$ may be visualized through:
ts = range(0, stop=2pi, length=100)
xs = sin.(ts)
ys = cos.(ts)
zs = ts
plot(xs, ys, zs)
How to plot the vector-valued function helix(t)
?
The output of helix.(ts)
is a vector of vectors:
ts = 1:4
helix.(ts)
4-element Array{Array{Float64,1},1}: [0.841471, 0.540302, 1.0] [0.909297, -0.416147, 2.0] [0.14112, -0.989992, 3.0] [-0.756802, -0.653644, 4.0]
This is in the wrong order, as we need xs
to be the first entry of each answer, ys
the second, and zs
the third. Some data reshaping is necessary. The following function xs_ys
extracts values for xs
, ys
, and optionally zs
, from a vector of vectors:
xs_ys(vs) = (A=hcat(vs...); Tuple([A[i,:] for i in eachindex(vs[1])]))
xs_ys (generic function with 1 method)
For our example, you can see how the values are collected:
xs_ys(helix.(ts))
([0.841471, 0.909297, 0.14112, -0.756802], [0.540302, -0.416147, -0.989992, -0.653644], [1.0, 2.0, 3.0, 4.0])
We add a few different calling methods, that may prove convenient:
xs_ys(v,vs...) = xs_ys([v, vs...])
xs_ys(r::Function, a, b, n=100) = (ts = range(a, stop=b, length=n); xs_ys(r.(ts)))
xs_ys (generic function with 4 methods)
The first allows the vectors to be specified as arguments to the function, the second a convenient wrapper to the process of creating some ts
over [a,b]
and then corresponding xs
, ys
, and (optionally) zs
.
(This function is in the MTH229
convenient package.)
With such a funtion, the above 3D graph could also have been generated through:
ts = range(0, stop=2pi, length=100)
xs, ys, zs = xs_ys(helix.(ts))
plot(xs, ys, zs)
The last two steps can be combined as:
ts = range(0, stop=2pi, length=100)
plot(xs_ys(helix.(ts))...)
Using "splatting" to create 3 arguments for plot
from the container with 3 components returned by xs_ys
.
Finally, these two steps could be just one, through the alternate iterface:
plot(xs_ys(helix, 0, 2pi)...)
As another example, this will plot the pringle
function:
ts = range(0, stop=2pi, length=100)
plot(xs_ys(pringle.(ts))...)
For function $f: R^1 \rightarrow R^1$ derivatives can be taken using automatic differentiation, as implemented in ForwardDiff
or by finite differencing. The latter allows the same definition in the vector valued case, so we will use that:
D1(r::Function, h=1e-4) = t -> (r(t+h) - r(t-h)) / (2h)
D1 (generic function with 2 methods)
r(t) = [sin(t), cos(t), t]
D1(r)(0)
3-element Array{Float64,1}: 0.9999999983333334 0.0 1.0
Though in general it can be a bad idea to iterate differencing, here the approximation isn't bad. This tests the maximum difference of the second derivative's error
ts = range(0, stop=2pi, length=1000)
D2(f) = D1(D1(f))
maximum([norm(D2(r)(t) - [-sin(t), -cos(t), 0]) for t in ts])
1.7411125000106153e-8
As with single variable calculus, the derivative may be found using automatic differentiation, implemented in the ForwardDiff
package. This is more accurate than the above and nearly as applicable. Here we overload the adjoint
function to provide the mathematical notation f'
:
using ForwardDiff
D(f::Function, n=1) = n > 1 ? D(D(f), n-1) : t -> ForwardDiff.derivative(f, float(t))
Base.adjoint(f::Function) = D(f)
note("""Again, this is in the `MTH229` convenience package. The adjoint function is understood to be a matrix operation and the postfix prime notation is an alias. Redefining it, as above, is generally considered bad form, as it will conflict with usage for arrays of functions. As this is atpyical usage for this task, and the familiar calculus notation is very useful, we do so here.
""")
<div class="alert alert-info" role="alert"> <div class="markdown"><p>Again, this is in the <code>MTH229</code> convenience package. The adjoint function is understood to be a matrix operation and the postfix prime notation is an alias. Redefining it, as above, is generally considered bad form, as it will conflict with usage for arrays of functions. As this is atpyical usage for this task, and the familiar calculus notation is very useful, we do so here.</p> </div> </div>
Then we have much improved accuracy:
ts = range(0, stop=2pi, length=1000)
maximum([norm(r''(t) - [-sin(t), -cos(t), 0]) for t in ts])
0.0
Adding an arrow to the plot is not as easy as it should be. In Plots
the quiver
function is used to add an arrow, but as of writing, there is no 3D support for this function. In the following, we use a line segment for that case. As well, the quiver
function is well suited for adding many arrows, but a bit cumbersome for just one, given data coming from a vector-valued function. This arrow!
function makes adding as single arrow fairly easy, by specifing the point to anchor it and the vector representing it.
function arrow!(plt::Plots.Plot, p, v; kwargs...)
if length(p) == 2
quiver!(plt, xs_ys([p])..., quiver=Tuple(xs_ys([v])); kwargs...)
elseif length(p) == 3
# 3d quiver needs support
# https://github.com/JuliaPlots/Plots.jl/issues/319#issue-159652535
# headless arrow instead
plot!(plt, xs_ys(p, p+v)...; kwargs...)
end
end
arrow!(p,v;kwargs...) = arrow!(Plots.current(), p, v; kwargs...)
arrow! (generic function with 2 methods)
Here we see how to use arrow!
by adding vectors along the pringle space curve:
ts = range(0, stop=2pi, length=100)
plt = plot(xs_ys(pringle.(ts))..., legend=false)
ts = range(0, stop=2pi, length=5)
for t0 in ts
p = pringle(t0)
T = uvec(pringle'(t0))
arrow!(plt, p, T, color="red")
end
plt
A smooth space curve has a tangent vector, as drawn, but also a normal and binormal vector for each time $t$. These may be difficult to compute by hand, but relatively straight forward numerically, as seen below.
The (unit) tangent vector is a the unit vector or $f'$, $T=f'/\|f'\|$. The (unit) normal vector is the unit vector of $T'$, $N=T'/\|T'\|$. Finally, the binormal is a unit vector perpendicular to both, or proportional to $T × N$.
r(t) = pringle(t)
a, b = 0, 2pi
t0 = 1.0
ts = range(a, stop=b, length=100)
plt = plot(xs_ys(r.(ts))...)
T(t) = uvec(r'(t))
N(t) = uvec(T'(t))
B(t) = uvec(T(t) × N(t))
p = pringle(t0)
arrow!(p, T(t0))
arrow!(p, N(t0))
arrow!(p, B(t0))
There is a fact for fixed-length vector functions (that is, $\|r(t)\|=c$ for all $t$) that the derivative and function are orthogonal. It should so so then that not only is $B$ orthogonal to $T$ and $N$, but also $N$ is orthogonal to $T$. This might be visible in the picture, here we verify it for each point in ts
:
maximum(abs.([T(t) ⋅ N(t) for t in ts]))
4.163336342344337e-16
This is not the case, were unit vectors not involved:
T1(t) = r'(t)
N1(t) = T1'(t)
maximum(abs.([T1(t) ⋅ N1(t) for t in ts]))
3.9994965106955003
The length of a parameterized curve is given by the formula
We can compute this easily in julia
. For our pringle we have the length of one turn is found with:
using QuadGK
r(t) = [cos(t), sin(t), sin(2t)]
ds(t) = norm(r'(t))
quadgk(ds, 0, 2pi)[1] # first value returned is answer (second is error)
10.540734326381326
This can also be done with with composition (∘
typed as \circ[tab]
):
quadgk(norm ∘ r', 0, 2pi)[1]
10.540734326381326
We put this into a function for convenient usage:
arclength(r::Function, a::Real, b::Real) = quadgk(norm ∘ r', a, b)[1]
arclength (generic function with 1 method)
Let $r(t)$ be a parameterization of a curve. The same space curve can have many different parameterizations, in fact any monotonically increasing function g
gives a new parameterization via r(g(t))
. Parameterizing by arclength is a useful way to talk about a specific parameter. Basically, this is a parameterization so the that arclength between 0 and $s$ is just $s$. Mathematically, the arclength function is monotonic, so we just need to take its inverse and call that g
. Finding the inverse is always possible due to monotonicity, but may not be algebraically possible. However, we have numeric tools to approximate it. For example,
r(t) = [cos(t), sin(t), sin(2t)] # pringle
s(t) = arclength(r, 0, t)
s (generic function with 1 method)
How to invert? We use the fzero
method from Roots
and bracketing, making an assumption that t
is positive and the arclength is less than 100:
using Roots
g(u) = fzero(t -> s(t) - u, (-1, 100))
g (generic function with 1 method)
Then the parameterization is found by composition
rs(t) = r(g(t))
rs (generic function with 1 method)
Here automatic differentiation won't work, as it doesn't work over g
, so we use D1
, defined above, to find the approximate derivative
quadgk(t -> norm(D1(rs)(t)), 0, 2)[1] # should be around 2
1.9999999911498048
The curvature can be computed different ways. One is that the curvature is the norm of the derivative of the tangent vector when the curve is parameterized by arc length. But more importantly, we have this formula which does not require a special parameterization, but does require a 3D vector-valued function:
kappa(t) = norm( r'(t) × r''(t) ) / norm(r'(t))^3
kappa (generic function with 1 method)
Many of the above computations can be done symbolically to match the work done with paper and pencil. Here we see how:
using SymPy
@vars t # create a symbolic variable
(t,)
With the function
r(t) = [sin(t), cos(t), t]
r (generic function with 1 method)
We can create a symbolic expression by evaluating r
at t
:
u = r(t)
The usual operations work as expected:
2u
u + u
u ⋅ u
u × u
More complicated expressions are possible
u ⋅ u - norm(u)^2
The above should be $0$ – for real valued vectors, but as SymPy
does not make that assumption unless asked. Here we try again:
@vars t real=true
u = r(t)
u ⋅ u - norm(u)^2
We can compute quantities symbolically. For example, the arclength (where we use t twice, though differently):
integrate(norm(diff.(u)), (t, 0, t)) |> simplify # integrate(ex, (var, a, b))
And the curvature:
κ = norm(diff.(u) × diff.(u,t,2)) / norm(diff.(u))^3 |> simplify
Of course, not all functions are so tractable as this example and numeric integration may still be required to get an answer.