57  Applications with scalar functions

This section uses these add-on packages:

using CalculusWithJulia
using Plots
plotly()
using SymPy
using Roots

And the following from the Contour package:

import Contour: contours, levels, level, lines, coordinates

This section presents different applications of scalar functions.

57.1 Tangent planes, linearization

Consider the case \(f:R^2 \rightarrow R\). We visualize \(z=f(x,y)\) through a surface. At a point \((a, b)\), this surface, if \(f\) is sufficiently smooth, can be approximated by a flat area, or a plane. For example, the Northern hemisphere of the earth, might be modeled simplistically by \(z = \sqrt{R^2 - (x^2 + y^2)}\) for some \(R\) and with the origin at the earth’s core. The ancient view of a “flat earth,” can be more generously seen as identifying this tangent plane with the sphere. More apt for current times, is the use of GPS coordinates to describe location. The difference between any two coordinates is technically a distance on a curved, nearly spherical, surface. But if the two points are reasonably closes (miles, not tens of miles) and accuracy isn’t of utmost importance (i.e., not used for self-driving cars), then the distance can be found from the Euclidean distance formula, \(\sqrt{(\Delta\text{latitude})^2 + (\Delta\text{longitude})^2}\). That is, as if the points were on a plane, not a curved surface.

For the univariate case, the tangent line has many different uses. Here we see the tangent plane also does.

57.1.1 Equation of the tangent plane

The partial derivatives have the geometric view of being the derivative of the univariate functions \(f(\vec\gamma_x(t))\) and \(f(\vec\gamma_y(t))\), where \(\vec\gamma_x\) moves just parallel to the \(x\) axis (e.g. \(\langle t + a, b\rangle\)). and \(\vec\gamma_y\) moves just parallel to the \(y\) axis. The partial derivatives then are slopes of tangent lines to each curve. The tangent plane, should it exist, should match both slopes at a given point. With this observation, we can identify it.

Consider \(f(\vec\gamma_x)\) at a point \((a,b)\). The path has a tangent vector, which has “slope” \(\frac{\partial f}{\partial x}\). and in the direction of the \(x\) axis, but not the \(y\) axis, as does this vector: \(\langle 1, 0, \frac{\partial f}{\partial x} \rangle\). Similarly, this vector \(\langle 0, 1, \frac{\partial f}{\partial y} \rangle\) describes the tangent line to \(f(\vec\gamma_y)\) at the point.

These two vectors will lie in the plane. The normal vector is found by their cross product:

@syms f_x f_y
n = [1, 0, f_x] × [0, 1, f_y]

\(\left[\begin{smallmatrix}- f_{x}\\- f_{y}\\1\end{smallmatrix}\right]\)

Let \(\vec{x} = \langle a, b, f(a,b) \rangle\). The tangent plane at \(\vec{x}\) then is described by all vectors \(\vec{v}\) with \(\vec{n}\cdot(\vec{v} - \vec{x}) = 0\). Using \(\vec{v} = \langle x,y,z\rangle\), we have:

\[ [-\frac{\partial f}{\partial x}, -\frac{\partial f}{\partial y}, 1] \cdot [x-a, y-b, z - f(a,b)] = 0, \]

or,

\[ z = f(a,b) + \frac{\partial f}{\partial x} (x-a) + \frac{\partial f}{\partial y} (y-b), \]

which is more compactly expressed as

\[ z = f(a,b) + \nabla(f) \cdot \langle x-a, y-b \rangle. \]

This form would then generalize to scalar functions from \(R^n \rightarrow R\). This is consistent with the definition of \(f\) being differentiable, where \(\nabla{f}\) plays the role of the slope in the formulas.

The following figure illustrates the above for the function \(f(x,y) = 6 - x^2 - y^2\):

f(x,y) = 6 - x^2 -y^2
f(x)= f(x...)

a,b = 1, -1/2


# draw surface
xr = 7/4
xs = ys = range(-xr, xr, length=100)
surface(xs, ys, f, legend=false)

# visualize tangent plane as 3d polygon
pt = [a,b]
tplane(x) = f(pt) + gradient(f)(pt) ⋅ (x - [a,b])

pts = [[a-1,b-1], [a+1, b-1], [a+1, b+1], [a-1, b+1], [a-1, b-1]]
plot!(unzip([[pt..., tplane(pt)] for pt in pts])...)

# plot paths in x and y direction through (a,b)
γ_x(t) = pt + t*[1,0]
γ_y(t) = pt + t*[0,1]

plot_parametric!((-xr-a)..(xr-a), t -> [γ_x(t)..., (f∘γ_x)(t)],  linewidth=3)
plot_parametric!((-xr-b)..(xr-b), t -> [γ_y(t)..., (f∘γ_y)(t)],  linewidth=3)

# draw directional derivatives in 3d and normal
pt = [a, b, f(a,b)]
fx, fy = gradient(f)(a,b)
arrow!(pt, [1, 0, fx], linewidth=3)
arrow!(pt, [0, 1, fy], linewidth=3)
arrow!(pt, [-fx, -fy, 1], linewidth=3) # normal

# draw point in base, x-y, plane
pt = [a, b, 0]
scatter!(unzip([pt])...)
arrow!(pt, [1,0,0], linestyle=:dash)
arrow!(pt, [0,1,0], linestyle=:dash)

Alternate forms

The equation for the tangent plane is often expressed in a more explicit form. For \(n=2\), if we set \(dx = x-a\) and \(dy=y-b\), then the equation for the plane becomes:

\[ f(a,b) + \frac{\partial f}{\partial x} dx + \frac{\partial f}{\partial y} dy, \]

which is a common form for the equation, though possibly confusing, as \(\partial x\) and \(dx\) need to be distinguished. For \(n > 2\), additional terms follow this pattern. This explicit form is helpful when doing calculations by hand, but much less so when working on the computer, say with Julia, as the representations using vectors (or matrices) can be readily implemented and their representation much closer to the formulas. For example, consider these two possible functions to find the tangent plane (returned as a function) at a point in \(2\) dimensions

function tangent_plane_1st_crack(f, pt)
  fx, fy = ForwardDiff.gradient(f, pt)
  x -> f(x...) + fx * (x[1]-pt[1]) + fy * (x[2]-pt[2])
end

tangent_plane_1st_crack (generic function with 1 method)

It isn’t so bad, but as written, we specialized to the number of dimensions, used indexing, and with additional dimensions, it clearly would get tedious to generalize. Using vectors, we might have:

function tangent_plane(f, pt)
  ∇f = ForwardDiff.gradient(f, pt) # using a variable ∇f
  x -> f(pt) + ∇f ⋅ (x - pt)
end

tangent_plane (generic function with 1 method)

This is much more like the compact formula and able to handle higher dimensions without rewriting.

57.1.2 Tangent plane for level curves

Consider the surface described by \(f(x,y,z) = c\), a constant. This is more general than surfaces described by \(z = f(x,y)\). The concept of a tangent plane should still be applicable though. Suppose, \(\vec{\gamma}(t)\) is a curve in the \(x-y-z\) plane, then we have \((f\circ\vec\gamma)(t)\) is a curve on the surface and its derivative is given by the chain rule through: \(\nabla{f}(\vec\gamma(t))\cdot \vec\gamma'(t)\). But this composition is constantly the same value, so the derivative is \(0\). This says that \(\nabla{f}(\vec\gamma(t))\) is orthogonal to \(\vec\gamma'(t)\) for any curve. As these tangential vectors to \(\vec\gamma\) lie in the tangent plane, the tangent plane can be characterized by having \(\nabla{f}\) as the normal.

This computation was previously done in two dimensions, and showed the gradient is orthogonal to the contour lines (and points in the direction of greatest ascent). It can be generalized to higher dimensions.

The surface \(F(x,y,z) = z - f(x,y) = 0\) has gradient given by \(\langle -\partial{f}/\partial{x}, -\partial{f}/\partial{y}, 1\rangle\), and as seen above, this vector is normal to the tangent plane, so this generalization agrees on the easier case.

For clarity:

• The scalar function \(z = f(x,y)\) describes a surface, \((x,y,f(x,y))\); the gradient, \(\nabla{f}\), is \(2\) dimensional and points in the direction of greatest ascent for the surface.
• The scalar function \(f(x,y,z)\) also describes a surface, through level curves \(f(x,y,z) = c\), for some constant \(c\). The gradient \(\nabla{f}\) is \(3\) dimensional and orthogonal to the surface.

Example

Let \(z = f(x,y) = \sin(x)\cos(x-y)\). Find an equation for the tangent plane at \((\pi/4, \pi/3)\).

We have many possible forms to express this in, but we will use the functional description:

@syms x, y

(x, y)

f(x,y) = sin(x) * cos(x-y)
f(x) = f(x...)
vars = [x, y]

gradf = diff.(f(x,y), vars)  # or use gradient(f, vars) or ∇((f,vars))

pt = [PI/4, PI/3]
gradfa = subs.(gradf, x=>pt[1], y=>pt[2])

f(pt) + gradfa ⋅ (vars - pt)

\(\left(x - \frac{\pi}{4}\right) \left(- \frac{\sqrt{2} \left(- \frac{\sqrt{6}}{4} + \frac{\sqrt{2}}{4}\right)}{2} + \frac{\sqrt{2} \left(\frac{\sqrt{2}}{4} + \frac{\sqrt{6}}{4}\right)}{2}\right) + \frac{\sqrt{2} \left(- \frac{\sqrt{6}}{4} + \frac{\sqrt{2}}{4}\right) \left(y - \frac{\pi}{3}\right)}{2} + \frac{\sqrt{2} \left(\frac{\sqrt{2}}{4} + \frac{\sqrt{6}}{4}\right)}{2}\)

Example

A cylinder \(f(x,y,z) = (x-a)^2 + y^2 = (2a)^2\) is intersected with a sphere \(g(x,y,z) = x^2 + y^2 + z^2 = a^2\). Let \(V\) be the line of intersection. (Viviani’s curve). Let \(P\) be a point on the curve. Describe the tangent to the curve.

We have the line of intersection will have tangent line lying in the tangent plane to both surfaces. These two surfaces have normal vectors given by the gradient, or \(\vec{n}_1 = \langle 2(x-a), 2y, 0 \rangle\) and \(\vec{n}_2 = \langle 2x, 2y, 2z \rangle\). The cross product of these two vectors will lie in both tangent planes, so we have:

\[ P + t (\vec{n}_1 \times \vec{n}_2), \]

will describe the tangent.

The curve may be described parametrically by \(\vec\gamma(t) = a \langle 1 + \cos(t), \sin(t), 2\sin(t/2) \rangle\). Let’s see that the above is correct by verifying that the cross product of the tangent vector computed two ways is \(0\):

a = 1
gamma(t) = a * [1 + cos(t), sin(t), 2sin(t/2) ]
P = gamma(1/2)
n1(x,y,z)= [2*(x-a), 2y, 0]
n2(x,y,z) = [2x,2y,2z]
n1(x) = n1(x...)
n2(x) = n2(x...)

t = 1/2
(n1(gamma(t)) × n2(gamma(t))) × gamma'(t)

3-element Vector{Float64}:
 0.0
 0.0
 0.0

Plotting level curves of \(F(x,y,z) = c\)

The wireframe plot can be used to visualize a surface of the type z=f(x,y), as previously illustrated. However we have no way of plotting \(3\)-dimensional implicit surfaces (of the type \(F(x,y,z)=c\)) as we do for \(2\)-dimensional implicit surfaces with Plots. (The MDBM or IntervalConstraintProgramming packages can be used along with Makie plotting package to produce one.)

The CalculusWithJulia package provides a stop-gap function, plot_implicit_surface for this task. The basic idea is to slice an axis, by default the \(z\) axis up and for each level plot the contours of \((x,y) \rightarrow f(x,y,z)-c\), which becomes a \(2\)-dimensional problem. The function allows any of 3 different axes to be chosen to slice over, the default being just the \(z\) axis.

We demonstrate with an example 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)

57.2 Linearization

The tangent plane is the best “linear approximation” to a function at a point. “Linear” refers to mathematical properties of the tangent plane, but at a practical level it means easy to compute, as it will involve only multiplication and addition. “Approximation” is useful in that if a bit of error is an acceptable tradeoff for computational ease, the tangent plane may be used in place of the function. In the univariate case, this is known as linearization, and the tradeoff is widely used in the derivation of theoretical relationships, as well as in practice to get reasonable numeric values.

Formally, this is saying:

\[ f(\vec{x}) \approx f(\vec{a}) + ∇f(\vec{a}) ⋅ (\vec{x} - \vec{a}). \]

The explicit meaning of \(\approx\) will be made clear when the generalization of Taylor’s theorem is to be stated.

Example: Linear approximation

The volume of a cylinder is \(V=\pi r^2 h\). It is thought a cylinder has \(r=1\) and \(h=2\). If instead, the amounts are \(r=1.01, h=2.01\), what is the difference in volume?

That is, if \(V(r,h) = \pi r^2 h\), what is \(V(1.01, 2.01) - V(1,2)\)?

We can use linear approximation to see that this difference is approximately \(\nabla{V} \cdot \langle 0.01, 0.01 \rangle\). This is:

V(r, h) = pi * r^2 * h
V(v) = V(v...)
a₁ = [1,2]
dx₁ = [0.01, 0.01]
ForwardDiff.gradient(V, a₁) ⋅ dx₁   # or use ∇(V)(a)

0.15707963267948966

The exact difference can be computed:

V(a₁ + dx₁) - V(a₁)

0.15833941133357854

Example

Let \(f(x,y) = \sin(\pi x y^2)\). Estimate \(f(1.1, 0.9)\).

Using linear approximation with \(dx=0.1\) and \(dy=-0.1\), this is

\[ f(1,1) + \nabla{f}(1,1) \cdot \langle 0.1, -0.1\rangle, \]

where \(f(1,1) = \sin(\pi) = 0\) and \(\nabla{f} = \langle \pi y^2\cos(\pi x y^2), \cos(\pi x y^2) 2\pi x y\rangle = \pi y \cos(\pi x y^2)\langle y,2x\rangle\). So, the answer is:

\[ 0 + \pi \cos(\pi) \langle 1,2\rangle\cdot \langle 0.1, -0.1 \rangle = (-\pi)(0.1 - 2(0.1)) = 0.1\pi. \]

Example

A piriform is described by the quartic surface \(f(x,y,z) = x^4 -x^3 + y^2+z^2 = 0\). Find the tangent line at the point \(\langle 2,2,2 \rangle\).

Here, \(\nabla{f}\) describes a normal to the tangent plane. The description of a plane may be described by \(\hat{N}\cdot(\vec{x} - \vec{x}_0) = 0\), where \(\vec{x}_0\) is identified with a point on the plane (the point \((2,2,2)\) here). With this, we have \(\hat{N}\cdot\vec{x} = ax + by + cz = \hat{N}\cdot\langle 2,2,2\rangle = 2(a+b+c)\). For this problem, \(\nabla{f}(2,2,2) = \langle a, b, c\rangle\) is given by:

f(x,y,z) = x^4 -x^3 + y^2 + z^2
f(v) = f(v...)
a, b,c = ∇(f)(2,2,2)
"$a x + $b y  + $c z = $([a,b,c] ⋅ [2,2,2])"

"20 x + 4 y  + 4 z = 56"

57.2.1 Newton’s method to solve \(f(x,y) = 0\) and \(g(x,y)=0\).

The level curve \(f(x,y)=0\) and the level curve \(g(x,y)=0\) may intersect. Solving algebraically for the intersection may be difficult in most cases, though the linear case is not. (The linear case being the intersection of two lines).

To elaborate, consider two linear equations written in a general form:

\[\begin{align*} ax + by &= u\\ cx + dy &= v \end{align*}\]

A method to solve this by hand would be to solve for \(y\) from one equation, replace this expression into the second equation and then solve for \(x\). From there, \(y\) can be found. A more advanced method expresses the problem in a matrix formulation of the form \(Mx=b\) and solves that equation. This form of solving is implemented in Julia, through the “backslash” operator. Here is the general solution:

@syms a b c d u v
M = [a b; c d]
B = [u, v]
M \ B .|> simplify

\(\left[\begin{smallmatrix}\frac{- b v + d u}{a d - b c}\\\frac{a v - c u}{a d - b c}\end{smallmatrix}\right]\)

The term \(\det(M) = ad-bc\) term is important, as evidenced by its appearance in the denominator of each term. When this is zero there is not a unique solution, as in the typical case.

Using Newton’s method to solve for intersection points, uses linearization of the surfaces to replace the problem to the intersection of level curves for tangent planes. This is the linear case that can be readily solved. As with Newton’s method for the univariate case, the new answer is generally a better approximation to the answer, and the process is iterated to get a good enough approximation, as defined through some tolerance.

Consider the functions \(f(x,y) =2 - x^2 - y^2\) and \(g(x,y) = 3 - 2x^2 - (1/3)y^2\). These graphs show their surfaces with the level sets for \(c=0\) drawn and just the levels sets, showing they intersect in \(4\) places.