using CalculusWithJulia
using Plots
plotly()
import Contour: contours, levels, level, lines, coordinates
using LinearAlgebra
using ForwardDiff59 2D and 3D plots in Julia with Plots
This section uses these add-on packages:
This covers plotting the typical 2D and 3D plots in Julia with the Plots package.
We will make use of some helper functions that will simplify plotting provided by the CalculusWithJulia package. As well, we will need to manipulate contours directly, so pull in the Contours package, using import to avoid name collisions and explicitly listing the methods we will use.
59.1 Parametrically described curves in space
Let \(r(t)\) be a vector-valued function with values in \(R^d\), \(d\) being \(2\) or \(3\). A familiar example is the equation for a line that travels in the direction of \(\vec{v}\) and goes through the point \(P\): \(r(t) = P + t \cdot \vec{v}\). A parametric plot over \([a,b]\) is the collection of all points \(r(t)\) for \(a \leq t \leq b\).
In Plots, parameterized curves can be plotted through two interfaces, here illustrated for \(d=2\): plot(f1, f2, a, b) or plot(xs, ys). The former is convenient for some cases, but typically we will have a function r(t) which is vector-valued, as opposed to a vector of functions. As such, we only discuss the latter.
An example helps illustrate. Suppose \(r(t) = \langle \sin(t), 2\cos(t) \rangle\) and the goal is to plot the full ellipse by plotting over \(0 \leq t \leq 2\pi\). As with plotting of curves, the goal would be to take many points between a and b and from there generate the \(x\) values and \(y\) values.
Let’s see this with 5 points, the first and last being identical due to the curve:
r₂(t) = [sin(t), 2cos(t)]
ts = range(0, stop=2pi, length=5)0.0:1.5707963267948966:6.283185307179586Then we can create the \(5\) points easily through broadcasting:
vs = r₂.(ts)5-element Vector{Vector{Float64}}:
[0.0, 2.0]
[1.0, 1.2246467991473532e-16]
[1.2246467991473532e-16, -2.0]
[-1.0, -3.6739403974420594e-16]
[-2.4492935982947064e-16, 2.0]This returns a vector of points (stored as vectors). The plotting function wants two collections: the set of \(x\) values for the points and the set of \(y\) values. The data needs to be generated differently or reshaped. The function unzip above takes data in this style and returns the desired format, returning a tuple with the \(x\) values and \(y\) values pulled out:
unzip(vs)([0.0, 1.0, 1.2246467991473532e-16, -1.0, -2.4492935982947064e-16], [2.0, 1.2246467991473532e-16, -2.0, -3.6739403974420594e-16, 2.0])To plot this, we “splat” the tuple so that plot gets the arguments separately:
plot(unzip(vs)...)This basic plot is lacking, of course, as there are not enough points. Using more initially is a remedy.
ts = range(0, 2pi, length=100)
plot(unzip(r₂.(ts))...)As a convenience, CalculusWithJulia provides plot_parametric to produce this plot. The interval is specified with the a..b notation of IntervalSets (which is available when the CalculusWithJulia package is loaded), the points to plot are adaptively chosen:
plot_parametric(0..2pi, r₂) # interval first59.1.1 Plotting a space curve in 3 dimensions
A parametrically described curve in 3D is similarly created. For example, a helix is described mathematically by \(r(t) = \langle \sin(t), \cos(t), t \rangle\). Here we graph two turns:
r₃(t) = [sin(t), cos(t), t]
plot_parametric(0..4pi, r₃)59.1.2 Adding a vector
The tangent vector indicates the instantaneous direction one would travel were they walking along the space curve. We can add a tangent vector to the graph. The quiver! function would be used to add a 2D vector, but Plots does not currently have a 3D analog. In addition, quiver! has a somewhat cumbersome calling pattern when adding just one vector. The CalculusWithJulia package defines an arrow! function that uses quiver for 2D arrows and a simple line for 3D arrows. As a vector incorporates magnitude and direction, but not a position, arrow! needs both a point for the position and a vector.
Here is how we can visualize the tangent vector at a few points on the helix:
plot_parametric(0..4pi, r₃, legend=false)
ts = range(0, 4pi, length=5)
for t in ts
arrow!(r₃(t), r₃'(t))
end
Note
Adding many arrows this way would be inefficient.
59.1.3 Setting a viewing angle for 3D plots
For 3D plots, the viewing angle can make the difference in visualizing the key features. In Plots, some backends allow the viewing angle to be set with the mouse by clicking and dragging. Not all do. For such, the camera argument is used, as in camera(azimuthal, elevation) where the angles are given in degrees. If the \(x\)-\(y\)-\(z\) coordinates are given, then elevation or inclination, is the angle between the \(z\) axis and the \(x-y\) plane (so 90 is a top view) and azimuthal is the angle in the \(x-y\) plane from the \(x\) axes.
59.2 Visualizing functions from \(R^2 \rightarrow R\)
If a function \(f: R^2 \rightarrow R\) then a graph of \((x,y,f(x,y))\) can be represented in 3D. It will form a surface. Such graphs can be most simply made by specifying a set of \(x\) values, a set of \(y\) values and a function \(f\), as with:
xs = range(-2, stop=2, length=100)
ys = range(-pi, stop=pi, length=100)
f(x,y) = x*sin(y)
surface(xs, ys, f)Rather than pass in a function, values can be passed in. Here they are generated with a list comprehension. The y values are innermost to match the graphic when passing in a function object:
zs = [f(x,y) for y in ys, x in xs]
surface(xs, ys, zs)Remembering if the ys or xs go first in the above can be hard. Alternatively, broadcasting can be used. The command f.(xs,ys) would return a vector, as the xs and ys match in shape–they are both column vectors. But the transpose of xs looks like a row vector and ys looks like a column vector, so broadcasting will create a matrix of values, as desired here:
surface(xs, ys, f.(xs', ys))This graph shows the tessalation algorithm. Here only the grid in the \(x\)-\(y\) plane is just one cell:
xs = ys = range(-1, 1, length=2)
f(x,y) = x*y
surface(xs, ys, f)A more accurate graph, can be seen here:
xs = ys = range(-1, 1, length=100)
f(x,y) = x*y
surface(xs, ys, f)59.2.1 Contour plots
Returning to the
The contour plot of \(f:R^2 \rightarrow R\) draws level curves, \(f(x,y)=c\), for different values of \(c\) in the \(x-y\) plane. They are produced in a similar manner as the surface plots:
xs = ys = range(-2,2, length=100)
f(x,y) = x*y
contour(xs, ys, f)The cross in the middle corresponds to \(c=0\), as when \(x=0\) or \(y=0\) then \(f(x,y)=0\).
Similarly, computed values for \(f(x,y)\) can be passed in. Here we change the function:
f(x,y) = 2 - (x^2 + y^2)
xs = ys = range(-2,2, length=100)
zs = [f(x,y) for y in ys, x in xs]
contour(xs, ys, zs)The chosen levels can be specified by the user through the levels argument, as in:
f(x,y) = 2 - (x^2 + y^2)
xs = ys = range(-2,2, length=100)
zs = [f(x,y) for y in ys, x in xs]
contour(xs, ys, zs, levels = [-1.0, 0.0, 1.0])If only a single level is desired, as scalar value can be specified. Though not with all backends for Plots. For example, this next graphic shows the \(0\)-level of the devil’s curve.
a, b = -1, 2
f(x,y) = y^4 - x^4 + a*y^2 + b*x^2
xs = ys = range(-5, stop=5, length=100)
contour(xs, ys, f, levels=[0.0])Contour plots are well known from the presence of contour lines on many maps. Contour lines indicate constant elevations. A peak is characterized by a series of nested closed paths. The following graph shows this for the peak at \((x,y)=(0,0)\).
xs = ys = range(-pi/2, stop=pi/2, length=100)
f(x,y) = sinc(sqrt(x^2 + y^2)) # sinc(x) is sin(x)/x
contour(xs, ys, f)Contour plots can be filled with colors through the contourf function:
xs = ys = range(-pi/2, stop=pi/2, length=100)
f(x,y) = sinc(sqrt(x^2 + y^2))
contourf(xs, ys, f)59.2.2 Combining surface plots and contour plots
In PyPlot it is possible to add a contour lines to the surface, or projected onto an axis. To replicate something similar, though not as satisfying, in Plots we use the Contour package.
f(x,y) = 2 + x^2 + y^2
xs = ys = range(-2, stop=2, length=100)
zs = [f(x,y) for y in ys, x in xs]
p = surface(xs, ys, zs, legend=false, fillalpha=0.5)
## we add to the graphic p, then plot
for cl in levels(contours(xs, ys, zs))
lvl = level(cl) # the z-value of this contour level
for line in lines(cl)
_xs, _ys = coordinates(line) # coordinates of this line segment
_zs = 0 * _xs
plot!(p, _xs, _ys, lvl .+ _zs, alpha=0.5) # add on surface
plot!(p, _xs, _ys, _zs, alpha=0.5) # add on x-y plane
end
end
pThere is no hidden line calculuation, in place we give the contour lines a transparency through the argument alpha=0.5.
59.2.3 Gradient and surface plots
The surface plot of \(f: R^2 \rightarrow R\) plots \((x, y, f(x,y))\) as a surface. The gradient of \(f\) is \(\langle \partial f/\partial x, \partial f/\partial y\rangle\). It is a two-dimensional object indicating the direction at a point \((x,y)\) where the surface has the greatest ascent. Illurating the gradient and the surface on the same plot requires embedding the 2D gradient into the 3D surface. This can be done by adding a constant \(z\) value to the gradient, such as \(0\).
f(x,y) = 2 - (x^2 + y^2)
xs = ys = range(-2, stop=2, length=100)
zs = [f(x,y) for y in ys, x in xs]
surface(xs, ys, zs, camera=(40, 25), legend=false)
p = [-1, 1] # in the region graphed, [-2,2] × [-2, 2]
f(x) = f(x...)
v = ForwardDiff.gradient(f, p)
# add 0 to p and v (two styles)
push!(p, -15)
scatter!(unzip([p])..., markersize=3)
v = vcat(v, 0)
arrow!(p, v)59.2.4 The tangent plane
Let \(z = f(x,y)\) describe a surface, and \(F(x,y,z) = f(x,y) - z\). The the gradient of \(F\) at a point \(p\) on the surface, \(\nabla F(p)\), will be normal to the surface and for a function, \(f(p) + \nabla f \cdot (x-p)\) describes the tangent plane. We can visualize each, as follows:
f(x,y) = 2 - x^2 - y^2
f(v) = f(v...)
F(x,y,z) = z - f(x,y)
F(v) = F(v...)
p = [1/10, -1/10]
global p1 = vcat(p, f(p...)) # note F(p1) == 0
global n⃗ = ForwardDiff.gradient(F, p1)
global tl(x) = f(p) + ForwardDiff.gradient(f, p) ⋅ (x - p)
tl(x,y) = tl([x,y])
xs = ys = range(-2, stop=2, length=100)
surface(xs, ys, f)
surface!(xs, ys, tl)
arrow!(p1, 5n⃗)From some viewing angles, the normal does not look perpendicular to the tangent plane. This is a quick verification for a randomly chosen point in the \(x-y\) plane:
a, b = randn(2)
dot(n⃗, (p1 - [a,b, tl(a,b)]))2.220446049250313e-1659.2.5 Parameterized surface plots
As illustrated, we can plot surfaces of the form \((x,y,f(x,y))\). However, not all surfaces are so readily described. For example, if \(F(x,y,z)\) is a function from \(R^3 \rightarrow R\), then \(F(x,y,z)=c\) is a surface of interest. For example, the sphere of radius one is a solution to \(F(x,y,z)=1\) where \(F(x,y,z) = x^2 + y^2 + z^2\).
Plotting such generally described surfaces is not so easy, but parameterized surfaces can be represented. For example, the sphere as a surface is not represented as a surface of a function, but can be represented in spherical coordinates as parameterized by two angles, essentially an “azimuth” and an “elevation”, as used with the camera argument.
Here we define functions that represent \((x,y,z)\) coordinates in terms of the corresponding spherical coordinates \((r, \theta, \phi)\).
# spherical: (radius r, inclination θ, azimuth φ)
X(r,theta,phi) = r * sin(theta) * sin(phi)
Y(r,theta,phi) = r * sin(theta) * cos(phi)
Z(r,theta,phi) = r * cos(theta)Z (generic function with 1 method)We can parameterize the sphere by plotting values for \(x\), \(y\), and \(z\) produced by a sequence of values for \(\theta\) and \(\phi\), holding \(r=1\):
thetas = range(0, stop=pi, length=50)
phis = range(0, stop=pi/2, length=50)
xs = [X(1, theta, phi) for theta in thetas, phi in phis]
ys = [Y(1, theta, phi) for theta in thetas, phi in phis]
zs = [Z(1, theta, phi) for theta in thetas, phi in phis]
surface(xs, ys, zs)
Note
The above may not work with all backends for Plots, even if those that support 3D graphics.
For convenience, the plot_parametric function from CalculusWithJulia can produce these plots using interval notation, a..b, and a function:
F(theta, phi) = [X(1, theta, phi), Y(1, theta, phi), Z(1, theta, phi)]
plot_parametric(0..pi, 0..pi/2, F)Example, the general cone
The general equation of cone consists of a vertex, \(V\), and a base curve. The cone consists of all line segments connecting \(V\) to the base curve, here parameterized and in the \(x-y\) plane: \(r(u) = \langle x(u), y(u), 0 \rangle\). The equations for the cone are:
\[ \frac{x - V_x}{x(u)-V_x} = \frac{y - V_y}{y(u)-V_y} = \frac{z - V_z}{z(u)-V_z} = t, \]
where \(t \in [0,1]\) This gives a vector-valued function
\[ F(u, t) = t \cdot (r(u) - V) + V, \quad \alpha \leq u \leq \beta, 0 \leq t \leq 1. \]
To illustrate, we have:
# cf. https://discourse.julialang.org/t/general-plotting-code-for-cone-in-3d-with-glmakie-or-plots/92104/3
basecurve(u) = [cos(u), sin(u) + sin(u/2), 0]
Vertex = [0, 3/4, 3]
Cone(u, t) = t * (basecurve(u) - Vertex) + Vertex
plot_parametric(0..2pi, 0..1, Cone)59.2.6 Plotting F(x,y, z) = c
There is no built in functionality in Plots to create surface described by \(F(x,y,z) = c\). An example of how to provide some such functionality for PyPlot appears here. The non-exported plot_implicit_surface function can be used to approximate this.
To use it, we see what happens when a sphere if rendered:
f(x,y,z) = x^2 + y^2 + z^2 - 25
CalculusWithJulia.plot_implicit_surface(f)This figure comes from a February 14, 2019 article in the New York Times. It shows an equation for a “heart,” as the graphic will illustrate:
a,b = 1,3
f(x,y,z) = (x^2+((1+b)*y)^2+z^2-1)^3-x^2*z^3-a*y^2*z^3
CalculusWithJulia.plot_implicit_surface(f, xlim=-2..2, ylim=-1..1, zlim=-1..2)