diff --git a/quarto/.gitignore b/quarto/.gitignore index b19a2da..ea9ed1f 100644 --- a/quarto/.gitignore +++ b/quarto/.gitignore @@ -5,4 +5,5 @@ /*/*.ipynb/ /*/bonepile.qmd /*/references.bib -weave_support.jl \ No newline at end of file +weave_support.jl +**/*.quarto_ipynb diff --git a/quarto/_quarto.yml b/quarto/_quarto.yml index 43b8ecb..2535ef9 100644 --- a/quarto/_quarto.yml +++ b/quarto/_quarto.yml @@ -1,4 +1,4 @@ -version: "0.25" +version: "0.27" engines: ['julia'] project: @@ -10,6 +10,7 @@ comments: book: title: "Calculus with Julia" author: "John Verzani" + google-analytics: "G-LXFE6MTM4M" date: now search: true repo-url: https://github.com/jverzani/CalculusWithJuliaNotes.jl @@ -117,7 +118,7 @@ book: - part: alternatives.qmd chapters: - - alternatives/giac.qmd + #- alternatives/giac.qmd - alternatives/symbolics.qmd - alternatives/SciML.qmd #- alternatives/interval_arithmetic.qmd diff --git a/quarto/alternatives.qmd b/quarto/alternatives.qmd index 862daed..d59ba08 100644 --- a/quarto/alternatives.qmd +++ b/quarto/alternatives.qmd @@ -2,8 +2,6 @@ These notes use a particular selection of packages. This selection could have been different. For example: -* The symbolic math is provided by `SymPy`. [Giac](./alternatives/giac.html) and -[Symbolics](./alternatives/symbolics.html) (along with `SymbolicUtils` and `ModelingToolkit`) provide alternatives. * The finding of zeros of scalar-valued, univariate functions is done with `Roots`. The [NonlinearSolve](./alternatives/SciML.html#nonlinearsolve) package provides an alternative for univariate and multi-variate functions. diff --git a/quarto/alternatives/symbolics.qmd b/quarto/alternatives/symbolics.qmd index aec506f..cac5003 100644 --- a/quarto/alternatives/symbolics.qmd +++ b/quarto/alternatives/symbolics.qmd @@ -1,6 +1,8 @@ # Symbolics.jl +XXX This needs updating! XXX +XXX add https://docs.sciml.ai/SymbolicIntegration/stable/ XXX There are a few options in `Julia` for symbolic math, for example, the `SymPy` package which wraps a Python library. This section describes a collection of native `Julia` packages providing many features of symbolic math. @@ -238,27 +240,6 @@ w = x^3 + y^3 - 2z^3 substitute(w, Dict(x=>2, y=>3)) ``` -The `fold` argument can be passed `false` to inhibit evaluation of values. Compare: - - -```{julia} -ex = 1 + sqrt(x) -substitute(ex, x=>2), substitute(ex, x=>2, fold=false) -``` - -Or - - -```{julia} -ex = sin(x) -substitute(ex, x=>π), substitute(ex, x=>π, fold=false) -``` - -For the latter, it is more efficient to directly use `Term`, which creates the symbolic expression representing the calling of `sin(π)`: - -```{julia} -Symbolics.Term(sin, [π]) -``` ### Simplify diff --git a/quarto/differentiable_vector_calculus/scalar_functions_applications-backup-question.qmd b/quarto/differentiable_vector_calculus/scalar_functions_applications-backup-question.qmd new file mode 100644 index 0000000..d7f0c5b --- /dev/null +++ b/quarto/differentiable_vector_calculus/scalar_functions_applications-backup-question.qmd @@ -0,0 +1,1433 @@ + + +##### Example: Steiner's problem + + +This is from [Strang](https://ocw.mit.edu/resources/res-18-001-calculus-online-textbook-spring-2005/textbook/MITRES_18_001_strang_13.pdf) p 506. + + +We have three points in the plane, $(x_1, y_1)$, $(x_2, y_2)$, and $(x_3,y_3)$. A point $p=(p_x, p_y)$ will have $3$ distances $d_1$, $d_2$, and $d_3$. Broadly speaking we want to minimize to find the point $p$ "nearest" the three fixed points within the triangle. Locating a facility so that it can service $3$ separate cities might be one application. The answer depends on the notion of what measure of distance to use. + + +If the measure is the Euclidean distance, then $d_i^2 = (p_x - x_i)^2 + (p_y - y_i)^2$. If we sought to minimize $d_1^2 + d_2^2 + d_3^2$, then we would proceed as follows: + + +```{julia} +@syms x1 y1 x2 y2 x3 y3 +d2(p,x) = (p[1] - x[1])^2 + (p[2]-x[2])^2 +d2_1, d2_2, d2_3 = d2((x,y), (x1, y1)), d2((x,y), (x2, y2)), d2((x,y), (x3, y3)) +exₛ = d2_1 + d2_2 + d2_3 +``` + +We then find the gradient, and solve for when it is $\vec{0}$: + + +```{julia} +gradfₛ = diff.(exₛ, [x,y]) +xstarₛ = solve(gradfₛ, [x,y]) +``` + +There is only one critical point, so must be a minimum. + + +We confirm this by looking at the Hessian and noting $H_{11} > 0$: + + +```{julia} +Hₛ = subs.(hessian(exₛ, [x,y]), x=>xstarₛ[x], y=>xstarₛ[y]) +``` + +As it occurs at $(\bar{x}, \bar{y})$ where $\bar{x} = (x_1 + x_2 + x_3)/3$ and $\bar{y} = (y_1+y_2+y_3)/3$---the averages of the three values---the critical point is an interior point of the triangle. + + +As mentioned by Strang, the real problem is to minimize $d_1 + d_2 + d_3$. A direct approach with `SymPy`---just replacing `d2` above with the square root fails. Consider instead the gradient of $d_1$, say. To avoid square roots, this is taken implicitly from $d_1^2$: + + +$$ +\frac{\partial}{\partial{x}}(d_1^2) = 2 d_1 \frac{\partial{d_1}}{\partial{x}}. +$$ + +But computing directly from the expression yields $2(x - x_1)$ Solving, yields: + + +$$ +\frac{\partial{d_1}}{\partial{x}} = \frac{(x-x_1)}{d_1}, \quad +\frac{\partial{d_1}}{\partial{y}} = \frac{(y-y_1)}{d_1}. +$$ + +The gradient is then $(\vec{p} - \vec{x}_1)/\|\vec{p} - \vec{x}_1\|$, a *unit* vector, call it $\hat{u}_1$. Similarly for $\hat{u}_2$ and $\hat{u}_3$. + + +Let $f = d_1 + d_2 + d_3$. Then $\nabla{f} = \hat{u}_1 + \hat{u}_2 + \hat{u}_3$. At the minimum, the gradient is $\vec{0}$, so the three unit vectors must cancel. This can only happen if the three make a "peace" sign with angles $120^\circ$ between them. To find the minimum then within the triangle, this point and the boundary must be considered, when this point falls outside the triangle. + + +Here is a triangle, where the minimum would be within the triangle: + + +```{julia} +usₛ = [[cos(t), sin(t)] for t in (0, 2pi/3, 4pi/3)] +polygon(ps) = unzip(vcat(ps, ps[1:1])) # easier way to plot a polygon + +pₛ = scatter([0],[0], markersize=2, legend=false, aspect_ratio=:equal) + +asₛ = (1,2,3) +plot!(polygon([a*u for (a,u) in zip(asₛ, usₛ)])...) +[arrow!([0,0], a*u, alpha=0.5) for (a,u) in zip(asₛ, usₛ)] +pₛ +``` + +For this triangle we find the Steiner point outside of the triangle. + + +```{julia} +asₛ₁ = (1, -1, 3) +scatter([0],[0], markersize=2, legend=false) +psₛₗ = [a*u for (a,u) in zip(asₛ₁, usₛ)] +plot!(polygon(psₛₗ)...) +``` + +Let's see where the minimum distance point is by constructing a plot. The minimum must be on the boundary, as the only point where the gradient vanishes is the origin, not in the triangle. The plot of the triangle has a contour plot of the distance function, so we see clearly that the minimum happens at the point `[0.5, -0.866025]`. On this plot, we drew the gradient at some points along the boundary. The gradient points in the direction of greatest increase---away from the minimum. That the gradient vectors have a non-zero projection onto the edges of the triangle in a direction pointing away from the point indicates that the function `d` would increase if moved along the boundary in that direction, as indeed it does. + + +```{julia} +euclid_dist(x; ps=psₛₗ) = sum(norm(x-p) for p in ps) +euclid_dist(x,y; ps=psₛₗ) = euclid_dist([x,y]; ps=ps) +``` + +```{julia} +#| hold: true +xs = range(-1.5, 1.5, length=100) +ys = range(-3, 1.0, length=100) + +p = plot(polygon(psₛₗ)..., linewidth=3, legend=false) +scatter!(p, unzip(psₛₗ)..., markersize=3) +contour!(p, xs, ys, euclid_dist) + +# add some gradients along boundary +li(t, p1, p2) = p1 + t*(p2-p1) # t in [0,1] +for t in range(1/100, 1/2, length=3) + pt = li(t, psₛₗ[2], psₛₗ[3]) + arrow!(pt, ForwardDiff.gradient(euclid_dist, pt)) + pt = li(t, psₛₗ[2], psₛₗ[1]) + arrow!(pt, ForwardDiff.gradient(euclid_dist, pt)) +end + +p +``` + +The following graph, shows distance along each edge: + + +```{julia} +#| hold : true +li(t, p1, p2) = p1 + t*(p2-p1) +p = plot(legend=false) +for i in 1:2, j in (i+1):3 + plot!(p, t -> euclid_dist(li(t, psₛₗ[i], psₛₗ[j]); ps=psₛₗ), 0, 1) +end +p +``` + +The smallest value is when $t=0$ or $t=1$, so at one of the points, as `li` is defined above. + + +##### Example: least squares + + +We know that two points determine a line. What happens when there are more than two points? This is common in statistics where a bivariate data set (pairs of points $(x,y)$) are summarized through a linear model $\mu_{y|x} = \alpha + \beta x$, That is the average value for $y$ given a particular $x$ value is given through the equation of a line. The data is used to identify what the slope and intercept are for this line. We consider a simple case---$3$ points. The case of $n \geq 3$ being similar. + + +We have a line $l(x) = \alpha + \beta(x)$ and three points $(x_1, y_1)$, $(x_2, y_2)$, and $(x_3, y_3)$. Unless these three points *happen* to be collinear, they can't possibly all lie on the same line. So to *approximate* a relationship by a line requires some inexactness. One measure of inexactness is the *vertical* distance to the line: + + +$$ +d1(\alpha, \beta) = |y_1 - l(x_1)| + |y_2 - l(x_2)| + |y_3 - l(x_3)|. +$$ + +Another might be the vertical squared distance to the line: + +$$ +\begin{align*} +d2(\alpha, \beta) &= (y_1 - l(x_1))^2 + (y_2 - l(x_2))^2 + (y_3 - l(x_3))^2 \\ +&= (y1 - (\alpha + \beta x_1))^2 + (y2 - (\alpha + \beta x_2))^2 + (y3 - (\alpha + \beta x_3))^2 +\end{align*} +$$ + +Another might be the *shortest* distance to the line: + + +$$ +d3(\alpha, \beta) = \frac{\beta x_1 - y_1 + \alpha}{\sqrt{1 + \beta^2}} + \frac{\beta x_2 - y_2 + \alpha}{\sqrt{1 + \beta^2}} + \frac{\beta x_3 - y_3 + \alpha}{\sqrt{1 + \beta^2}}. +$$ + +The method of least squares minimizes the second one of these. That is, it chooses $\alpha$ and $\beta$ that make the expression a minimum. + + +```{julia} +@syms xₗₛ[1:3] yₗₛ[1:3] α β +li(x, alpha, beta) = alpha + beta * x +d₂(alpha, beta) = sum((y - li(x, alpha, beta))^2 for (y,x) in zip(yₗₛ, xₗₛ)) +d₂(α, β) +``` + +To identify $\alpha$ and $\beta$ we find the gradient: + + +```{julia} +grad_d₂ = diff.(d₂(α, β), [α, β]) +``` + +```{julia} +outₗₛ = solve(grad_d₂, [α, β]) +``` + +As found, the formulas aren't pretty. If $x_1 + x_2 + x_3 = 0$ they simplify. For example: + + +```{julia} +subs(outₗₛ[β], sum(xₗₛ) => 0) +``` + +Let $\vec{x} = \langle x_1, x_2, x_3 \rangle$ and $\vec{y} = \langle y_1, y_2, y_3 \rangle$ this is simply $(\vec{x} \cdot \vec{y})/(\vec{x}\cdot \vec{x})$, a formula that will generalize to $n > 3$. The assumption is not a restriction---it comes about by subtracting the mean, $\bar{x} = (x_1 + x_2 + x_3)/3$, from each $x$ term (and similarly subtract $\bar{y}$ from each $y$ term). A process called "centering." + + +With this observation, the formulas can be re-expressed through: + + +$$ +\beta = \frac{\sum{(x_i - \bar{x})(y_i - \bar{y})}}{\sum(x_i-\bar{x})^2}, +\quad +\alpha = \bar{y} - \beta \bar{x}. +$$ + +Relative to the centered values, this may be viewed as a line through $(\bar{x}, \bar{y})$ with slope given by $(\vec{x}-\bar{x})\cdot(\vec{y}-\bar{y}) / \|\vec{x}-\bar{x}\|^2$. + + +As an example, if the point are $(1,1), (2,3), (5,8)$ we get: + + +```{julia} +[k => subs(v, xₗₛ[1]=>1, yₗₛ[1]=>1, xₗₛ[2]=>2, yₗₛ[2]=>3, + xₗₛ[3]=>5, yₗₛ[3]=>8) for (k,v) in outₗₛ] +``` + +### Gradient descent + + +As seen in the examples above, extrema may be identified analytically by solving for when the gradient is $0$. Here we discuss some numeric algorithms for finding extrema. + + +An algorithm to identify where a surface is at its minimum is [gradient descent](https://en.wikipedia.org/wiki/Gradient_descent). The gradient points in the direction of the steepest ascent of the surface and the negative gradient the direction of the steepest descent. To move to a minimum then, it make intuitive sense to move in the direction of the negative gradient. How far? That is a different question and one with different answers. Let's formulate the movement first, then discuss how far. + + +Let $\vec{x}_0$, $\vec{x}_1$, $\dots$, $\vec{x}_n$ be the position of the algorithm for $n$ steps starting from an initial point $\vec{x}_0$. The difference between these points is given by: + + +$$ +\vec{x}_{n+1} = \vec{x}_n - \gamma \nabla{f}(\vec{x}_n), +$$ + +where $\gamma$ is some scaling factor for the gradient. The above quantifies the idea: to go from $\vec{x}_n$ to $\vec{x}_{n+1}$, move along $-\nabla{f}$ by a certain amount. + + +Let $\Delta_x =\vec{x}_{n}- \vec{x}_{n-1}$ and $\Delta_y = \nabla{f}(\vec{x}_{n}) - \nabla{f}(\vec{x}_{n-1})$ A variant of the Barzilai-Borwein method is to take $\gamma_n = | \Delta_x \cdot \Delta_y / \Delta_y \cdot \Delta_y |$. + + +To illustrate, take $f(x,y) = - e^{-((x-1)^2 + 2(y-1/2)^2)}$ and a starting point $\langle 0, 0 \rangle$. We have, starting with $\gamma_0 = 1$ there are $5$ steps taken: + + +```{julia} +f₂(x,y) = -exp(-((x-1)^2 + 2(y-1/2)^2)) +f₂(x) = f₂(x...) + +xs₂ = [[0.0, 0.0]] # we store a vector +gammas₂ = [1.0] + +for n in 1:5 + xn = xs₂[end] + gamma₀ = gammas₂[end] + xn1 = xn - gamma₀ * gradient(f₂)(xn) + dx, dy = xn1 - xn, gradient(f₂)(xn1) - gradient(f₂)(xn) + gamman1 = abs( (dx ⋅ dy) / (dy ⋅ dy) ) + + push!(xs₂, xn1) + push!(gammas₂, gamman1) +end + +[(x, f₂(x)) for x in xs₂] +``` + +We now visualize, using the `Contour` package to draw the contour lines in the $x-y$ plane: + + +```{julia} +#| hold: true +function surface_contour(xs, ys, f; offset=0) + p = surface(xs, ys, f, legend=false, fillalpha=0.5) + + ## we add to the graphic p, then plot + zs = [f(x,y) for x in xs, y in ys] # reverse order for use with Contour package + 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 = offset * _xs + plot!(p, _xs, _ys, _zs, alpha=0.5) # add curve on x-y plane + end + end + p +end + + +offset = 0 +us = vs = range(-1, 2, length=100) +surface_contour(us, vs, f₂, offset=offset) +pts = [[pt..., offset] for pt in xs₂] +scatter3d!(unzip(pts)...) +plot!(unzip(pts)..., linewidth=3) +``` + +### Newton's method for minimization + + +A variant of Newton's method can be used to minimize a function $f:R^2 \rightarrow R$. We look for points where both partial derivatives of $f$ vanish. Let $g(x,y) = \partial f/\partial x(x,y)$ and $h(x,y) = \partial f/\partial y(x,y)$. Then applying Newton's method, as above to solve simultaneously for when $g=0$ and $h=0$, we considered this matrix: + + +$$ +M = [\nabla{g}'; \nabla{h}'], +$$ + +and had a step expressible in terms of the inverse of $M$ as $M^{-1} [g; h]$. In terms of the function $f$, this step is $H^{-1}\nabla{f}$, where $H$ is the Hessian matrix. [Newton](https://en.wikipedia.org/wiki/Newton%27s_method_in_optimization#Higher_dimensions)'s method then becomes: + + +$$ +\vec{x}_{n+1} = \vec{x}_n - [H_f(\vec{x}_n)]^{-1} \nabla(f)(\vec{x}_n). +$$ + +The Wikipedia page states where applicable, Newton's method converges much faster towards a local maximum or minimum than gradient descent. + + +We apply it to the task of characterizing the following function, which has a few different peaks over the region $[-3,3] \times [-2,2]$: + + +```{julia} +function peaks(x, y) + z = 3 * (1 - x)^2 * exp(-x^2 - (y + 1)^2) + z += -10 * (x / 5 - x^3 - y^5) * exp(-x^2 - y^2) + z += -1/3 * exp(-(x+1)^2 - y^2) + return z +end +peaks(v) = peaks(v...) +``` + +```{julia} +#| hold: true +xs = range(-3, stop=3, length=100) +ys = range(-2, stop=2, length=100) +Ps = surface(xs, ys, peaks, legend=false) +Pc = contour(xs, ys, peaks, legend=false) +plot(Ps, Pc, layout=2) # combine plots +``` + +As we will solve for the critical points numerically, we consider the contour plot as well, as it shows better where the critical points are. + + +Over this region we see clearly 5 peaks or valleys: near $(0, 1.5)$, near $(1.2, 0)$, near $(0.2, -1.8)$, near $(-0.5, -0.8)$, and near $(-1.2, 0.2)$. To classify the $5$ critical points we need to first identify them, then compute the Hessian, and then, possibly compute $f_{xx}$ at the point. Here we do so for one of them using a numeric approach. + + +For concreteness, consider the peak or valley near $(0,1.5)$. We use Newton's method to numerically compute the critical point. The Newton step, specialized here is: + + +```{julia} +function newton_stepₚ(f, x) + M = ForwardDiff.hessian(f, x) + b = ForwardDiff.gradient(f, x) + x - M \ b +end +``` + +We perform $3$ steps of Newton's method, and see that it has found a critical point. + + +```{julia} +xₚ = [0, 1.5] +xₚ = newton_stepₚ(peaks, xₚ) +xₚ = newton_stepₚ(peaks, xₚ) +xₚ = newton_stepₚ(peaks, xₚ) +xₚ, ForwardDiff.gradient(peaks, xₚ) +``` + +The Hessian at this point is given by: + + +```{julia} +Hₚ = ForwardDiff.hessian(peaks, xₚ) +``` + +From which we see: + + +```{julia} +#| hold: true +fxx = Hₚ[1,1] +d = det(Hₚ) +fxx, d +``` + +Consequently we have a local maximum at this critical point. + + +:::{.callout-note} +## Note + +::: + +The `Optim.jl` package provides efficient implementations of these two numeric methods, and others. + + +## Constrained optimization, Lagrange multipliers + + +We considered the problem of maximizing a function over a closed region. This maximum is achieved at a critical point *or* a boundary point. Investigating the critical points isn't so difficult and the second partial derivative test can help characterize the points along the way, but characterizing the boundary points usually involves parameterizing the boundary, which is not always so easy. However, if we put this problem into a more general setting a different technique becomes available. + + +The different setting is: maximize $f(x,y)$ subject to the constraint $g(x,y) = k$. The constraint can be used to describe the boundary used previously. + + +Why does this help? The key is something we have seen prior: If $g$ is differentiable, and we take $\nabla{g}$, then it will point at directions *orthogonal* to the level curve $g(x,y) = 0$. (Parameterize the curve, then $(g\circ\vec{r})(t) = 0$ and so the chain rule has $\nabla{g}(\vec{r}(t)) \cdot \vec{r}'(t) = 0$.) For example, consider the function $g(x,y) = x^2 +2y^2 - 1$. The level curve $g(x,y) = 0$ is an ellipse. Here we plot the level curve, along with a few gradient vectors at points satisfying $g(x,y) = 0$: + + +```{julia} +#| hold: true +g(x,y) = x^2 + 2y^2 -1 +g(v) = g(v...) + +xs = range(-3, 3, length=100) +ys = range(-1, 4, length=100) + +p = plot(aspect_ratio=:equal, legend=false) +contour!(xs, ys, g, levels=[0]) + +gi(x) = sqrt(1/2*(1-x^2)) # solve for y in terms of x +pts = [[x, gi(x)] for x in (-3/4, -1/4, 1/4, 3/4)] + +for pt in pts + arrow!(pt, ForwardDiff.gradient(g, pt) ) +end + +p +``` + +From the plot we see the key property that $\nabla g$ is orthogonal to the level curve. + + +Now consider $f(x,y)$, a function we wish to maximize. The gradient points in the direction of *greatest* increase, provided $f$ is smooth. We are interested in the value of this gradient along the level curve of $g$. Consider this figure representing a portion of the level curve, it's tangent, normal, the gradient of $f$, and the contours of $f$: + + +```{julia} +#| hold: true +#| echo: false +r(t) = [cos(t), sin(t)/2] +plot_parametric(pi/12..pi/3, r, legend=false, aspect_ratio=true, linewidth=3) +T(t) = -r'(t) / norm(r'(t)) +No(t) = T'(t) / norm(T'(t)) +t = pi/4 +lambda=1/10 +scatter!(unzip([r(t)])...) +arrow!(r(t), T(t)*lambda) +arrow!(r(t), No(t)* lambda) + +f(x,y)= x^2 + y^2 +f(v) = f(v...) +arrow!(r(t), lambda*ForwardDiff.gradient(f, r(t))) + +xs = range(0.5,1, length=100) +ys = range(0.1, 0.5, length=100) +contour!(xs, ys, f) +``` + +We can identify the tangent, the normal, and subsequently the gradient of $f$. Is the point drawn a maximum of $f$ subject to the constraint $g$? + + +The answer is no, but why? By adding the contours of $f$, we see that moving along the curve from this point will increase or decrease $f$, depending on which direction we move in. As the *gradient* is the direction of greatest increase, we can see that the *projection* of the gradient on the tangent will point in a direction of *increase*. + + +It isn't just because the point picked was chosen to make a pretty picture, and not be a maximum. Rather, the fact that $\nabla{f}$ has a non-trivial projection onto the tangent vector. What does it say if we move the point in the direction of this projection? + + +The gradient points in the direction of greatest increase. If we first move in one component of the gradient we will increase, just not as fast. This is because the directional derivative in the direction of the tangent will be non-zero. In the picture, if we were to move the point to the right along the curve $f(x,y)$ will increase. + + +Now consider this figure at a different point of the figure: + + +```{julia} +#| hold: true +#| echo: false +r(t) = [cos(t), sin(t)/2] +plot_parametric(-pi/6..pi/6,r, legend=false, aspect_ratio=true, linewidth=3) +T(t) = -r'(t) / norm(r'(t)) +No(t) = T'(t) / norm(T'(t)) +t = 0 +lambda=1/10 +scatter!(unzip([r(t)])...) +arrow!(r(t), T(t)*lambda) +arrow!(r(t), No(t)* lambda) + +f(x,y)= x^2 + y^2 +f(v) = f(v...) +arrow!(r(t), lambda*ForwardDiff.gradient(f, r(t))) + +xs = range(0.5,1.5, length=100) +ys = range(-0.5, 0.5, length=100) +contour!(xs, ys, f, levels = [.7, .85, 1, 1.15, 1.3]) +``` + +We can still identify the tangent and normal directions. What is different about this point is that local movement on the constraint curve is also local movement on the contour line of $f$, so $f$ doesn't increase or decrease here, as it would if this point were an extrema along the constraint. The key to seeing this is the contour lines of $f$ are *tangent* to the constraint. The respective gradients are *orthogonal* to their tangent lines, and in dimension $2$, this implies they are parallel to each other. + +::: {.callout-note icon=false} +## The method of Lagrange multipliers + +To optimize $f(x,y)$ subject to a constraint $g(x,y) = k$ we solve for all *simultaneous* solutions to + +$$ +\begin{align*} +\nabla{f}(x,y) &= \lambda \nabla{g}(x,y), \text{and}\\ +g(x,y) &= k. +\end{align*} +$$ + + +These *possible* points are evaluated to see if they are maxima or minima. + +::: + + +The method will not work if $\nabla{g} = \vec{0}$ or if $f$ and $g$ are not differentiable. + + +--- + + +##### Example + + +We consider [again]("../derivatives/optimization.html") the problem of maximizing all rectangles subject to the perimeter being $20$. We have seen this results in a square. This time we use the Lagrange multiplier technique. We have two equations: + + +$$ +A(x,y) = xy, \quad P(x,y) = 2x + 2y = 20. +$$ + +We see $\nabla{A} = \lambda \nabla{P}$, or $\langle y, x \rangle = \lambda \langle 2, 2\rangle$. We see the solution has $x = y$ and from the constraint $x=y = 5$. + + +This is clearly the maximum for this problem, though the Lagrange technique does not imply that, it only identifies possible extrema. + + +##### Example + + +We can reverse the question: what are the ranges for the perimeter when the area is a fixed value of $25$? We have: + + +$$ +P(x,y) = 2x + 2y, \quad A(x,y) = xy = 25. +$$ + +Now we look for $\nabla{P} = \lambda \nabla{A}$ and will get, as the last example, that $\langle 2, 2 \rangle = \lambda \langle y, x\rangle$. So $x=y$ and from the constraint $x=y=5$. + + +However this is *not* the maximum perimeter, but rather the minimal perimeter. The maximum is $\infty$, which comes about in the limit by considering long skinny rectangles. + + +##### Example: A rephrasing + + +An slightly different formulation of the Lagrange method is to combine the equation and the constraint into one equation: + + +$$ +L(x,y,\lambda) = f(x,y) - \lambda (g(x,y) - k). +$$ + +The we have + + +$$ +\begin{align*} +\frac{\partial L}{\partial{x}} &= \frac{\partial{f}}{\partial{x}} - \lambda \frac{\partial{g}}{\partial{x}}\\ +\frac{\partial L}{\partial{y}} &= \frac{\partial{f}}{\partial{y}} - \lambda \frac{\partial{g}}{\partial{y}}\\ +\frac{\partial L}{\partial{\lambda}} &= 0 + (g(x,y) - k). +\end{align*} +$$ + + +But if the Lagrange condition holds, each term is $0$, so Lagrange's method can be seen as solving for point $\nabla{L} = \vec{0}$. The optimization problem in two variables with a constraint becomes a problem of finding and classifying zeros of a function with *three* variables. + + +Apply this to the optimization problem: + + +Find the extrema of $f(x,y) = x^2 - y^2$ subject to the constraint $g(x,y) = x^2 + y^2 = 1$. + + +We have: + + +$$ +L(x, y, \lambda) = f(x,y) - \lambda(g(x,y) - 1) +$$ + +We can solve for $\nabla{L} = \vec{0}$ by hand, but we do so symbolically: + + +```{julia} +@syms lambda +fₗₐ(x, y) = x^2 - y^2 +gₗₐ(x, y) = x^2 + y^2 +Lₗₐ(x, y, lambda) = fₗₐ(x,y) - lambda * (gₗₐ(x,y) - 1) +dsₗₐ = solve(diff.(Lₗₐ(x, y, lambda), [x, y, lambda])) +``` + +This has $4$ easy solutions, here are the values at each point: + + +```{julia} +[fₗₐ(d[x], d[y]) for d in dsₗₐ] +``` + +So $1$ is a maximum value and $-1$ a minimum value. + + +##### Example: Dido's problem + + +Consider a slightly different problem: What shape should a rope (curve) of fixed length make to *maximize* the area between the rope and $x$ axis? + + +Let $L$ be the length of the rope and suppose $y(x)$ describes the curve. Then we wish to + + +$$ +\text{Maximize } \int y(x) dx, \quad\text{subject to } +\int \sqrt{1 + y'(x)^2} dx = L. +$$ + +The latter being the formula for arc length. This is very much like an optimization problem that Lagrange's method could help solve, but with one big difference: the answer is *not* a point but a *function*. + + +This is a variant of [Dido](http://www.ams.org/publications/journals/notices/201709/rnoti-p980.pdf)'s problem, described by Bandle as + + +> *Dido’s problem*: The Roman poet Publius Vergilius Maro (70–19 B.C.) tells in his epic Aeneid the story of queen Dido, the daughter of the Phoenician king of the 9th century B.C. After the assassination of her husband by her brother she fled to a haven near Tunis. There she asked the local leader, Yarb, for as much land as could be enclosed by the hide of a bull. Since the deal seemed very modest, he agreed. Dido cut the hide into narrow strips, tied them together and encircled a large tract of land which became the city of Carthage. Dido faced the following mathematical problem, which is also known as the isoperimetric problem: Find among all curves of given length the one which encloses maximal area. Dido found intuitively the right answer. + + + +The problem as stated above and method of solution follows notes by [Wang](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.368.1522&rep=rep1&type=pdf) though Bandle attributes the ideas back to a 19-year old Lagrange in a letter to Euler. + + +The method of solution will be to *assume* we have the function and then characterize this function in such a way that it can be identified. + + +Following Lagrange, we generalize the problem to the following: maximize $\int_{x_0}^{x_1} f(x, y(x), y'(x)) dx$ subject to a constraint $\int_{x_0}^{x_1} g(x,y(x), y'(x)) dx = K$. Suppose $y(x)$ is a solution. + + +The starting point is a *perturbation*: $\hat{y}(x) = y(x) + \epsilon_1 \eta_1(x) + \epsilon_2 \eta_2(x)$. There are two perturbation terms, were only one term added, then the perturbation may make $\hat{y}$ not satisfy the constraint, the second term is used to ensure the constraint is not violated. If $\hat{y}$ is to be a possible solution to our problem, we would want $\hat{y}(x_0) = \hat{y}(x_1) = 0$, as it does for $y(x)$, so we *assume* $\eta_1$ and $\eta_2$ satisfy this boundary condition. + + +With this notation, and fixing $y$ we can re-express the equations in terms of $\epsilon_1$ and $\epsilon_2$: + + +$$ +\begin{align*} +F(\epsilon_1, \epsilon_2) &= \int f(x, \hat{y}, \hat{y}') dx = +\int f(x, y + \epsilon_1 \eta_1 + \epsilon_2 \eta_2, y' + \epsilon_1 \eta_1' + \epsilon_2 \eta_2') dx,\\ +G(\epsilon_1, \epsilon_2) &= \int g(x, \hat{y}, \hat{y}') dx = +\int g(x, y + \epsilon_1 \eta_1 + \epsilon_2 \eta_2, y' + \epsilon_1 \eta_1' + \epsilon_2 \eta_2') dx. +\end{align*} +$$ + + +Then our problem is restated as: + + +$$ +\text{Maximize } F(\epsilon_1, \epsilon_2) \text{ subject to } +G(\epsilon_1, \epsilon_2) = L. +$$ + +Now, Lagrange's method can be employed. This will be fruitful---even though we know the answer---it being $\epsilon_1 = \epsilon_2 = 0$! + + +Forging ahead, we compute $\nabla{F}$ and $\lambda \nabla{G}$ and set $\epsilon_1 = \epsilon_2 = 0$ where the two are equal. This will lead to a description of $y$ in terms of $y'$. + + +Lagrange's method has: + + +$$ +\frac{\partial{F}}{\partial{\epsilon_1}}(0,0) - \lambda \frac{\partial{G}}{\partial{\epsilon_1}}(0,0) = 0, \text{ and } +\frac{\partial{F}}{\partial{\epsilon_2}}(0,0) - \lambda \frac{\partial{G}}{\partial{\epsilon_2}}(0,0) = 0. +$$ + +Computing just the first one, we have using the chain rule and assuming interchanging the derivative and integral is possible: + + +$$ +\begin{align*} +\frac{\partial{F}}{\partial{\epsilon_1}} +&= \int \frac{\partial}{\partial{\epsilon_1}}( +f(x, y + \epsilon_1 \eta_1 + \epsilon_2 \eta_2, y' + \epsilon_1 \eta_1' + \epsilon_2 \eta_2')) dx\\ +&= \int \left(\frac{\partial{f}}{\partial{y}} \eta_1 + \frac{\partial{f}}{\partial{y'}} \eta_1'\right) dx\quad\quad(\text{from }\nabla{f} \cdot \langle 0, \eta_1, \eta_1'\rangle)\\ +&=\int \eta_1 \left(\frac{\partial{f}}{\partial{y}} - \frac{d}{dx}\frac{\partial{f}}{\partial{y'}}\right) dx. +\end{align*} +$$ + + +The last line by integration by parts: +$\int u'(x) v(x) dx = (u \cdot v)(x)\mid_{x_0}^{x_1} - \int u(x) \frac{d}{dx} v(x) dx = - \int u(x) \frac{d}{dx} v(x) dx$. +The last lines, as $\eta_1 = 0$ at $x_0$ and $x_1$ by assumption. We get: + +$$ +0 = \int \eta_1\left(\frac{\partial{f}}{\partial{y}} - \frac{d}{dx}\frac{\partial{f}}{\partial{y'}}\right). +$$ + +Similarly were $G$ considered, we would find a similar statement. Setting $L(x, y, y') = f(x, y, y') - \lambda g(x, y, y')$, the combination of terms gives: + + +$$ +0 = \int \eta_1\left(\frac{\partial{L}}{\partial{y}} - \frac{d}{dx}\frac{\partial{L}}{\partial{y'}}\right) dx. +$$ + +Since $\eta_1$ is arbitrary save for its boundary conditions, under smoothness conditions on $L$ this will imply the rest of the integrand *must* be $0$. + + +That is, If $y(x)$ is a maximizer of $\int_{x_0}^{x_1} f(x, y, y')dx$ and sufficiently smooth over $[x_0, x_1]$ and $y(x)$ satisfies the constraint $\int_{x_0}^{x_1} g(x, y, y')dx = K$ then there exists a constant $\lambda$ such that $L = f -\lambda g$ will satisfy: + + +$$ +\frac{d}{dx}\frac{\partial{L}}{\partial{y'}} - \frac{\partial{L}}{\partial{y}} = 0. +$$ + +If $\partial{L}/\partial{x} = 0$, this simplifies to the [Beltrami](https://en.wikipedia.org/wiki/Beltrami_identity) identity: + + +$$ +L - y' \frac{\partial{L}}{\partial{y'}} = C.\quad(\text{Beltrami identity}) +$$ + +--- + + +For Dido's problem, $f(x,y,y') = y$ and $g(x, y, y') = \sqrt{1 + y'^2}$, so $L = y - \lambda\sqrt{1 + y'^2}$ will have $0$ partial derivative with respect to $x$. Using the Beltrami identify we have: + + +$$ +(y - \lambda\sqrt{1 + y'^2}) + \lambda y' \frac{2y'}{2\sqrt{1 + y'^2}} = C. +$$ + +by multiplying through by the denominator and squaring to remove the square root, a quadratic equation in $y'^2$ can be found. This can be solved to give: + + +$$ +y' = \frac{dy}{dx} = \sqrt{\frac{\lambda^2 -(y - C)^2}{(y-C)^2}}. +$$ + +Here is a snippet of `SymPy` code to verify the above: + + +```{julia} +#| hold: true +@syms y y′ λ C +ex = Eq(-λ*y′^2/sqrt(1 + y′^2) + λ*sqrt(1 + y′^2), y - C) +Δ = sqrt(1 + y′^2) / (y - C) +ex1 = Eq(simplify(ex.lhs()*Δ), simplify(ex.rhs() * Δ)) +ex2 = Eq(ex1.lhs()^2 - 1, simplify(ex1.rhs()^2) - 1) +``` + +Now $y'$ can be integrated using the substitution $y - C = \lambda \cos\theta$ to give: $-\lambda\int\cos\theta d\theta = x + D$, $D$ some constant. That is: + + +$$ +\begin{align*} +x + D &= - \lambda \sin\theta\\ +y - C &= \lambda\cos\theta. +\end{align*} +$$ + + +Squaring gives the equation of a circle: $(x +D)^2 + (y-C)^2 = \lambda^2$. + + +We center and *rescale* the problem so that $x_0 = -1, x_1 = 1$. Then $L > 2$ as otherwise the rope is too short. From here, we describe the radius and center of the circle. + + +We have $y=0$ at $x=1$ and $-1$ giving: + + +$$ +\begin{align*} +(-1 + D)^2 + (0 - C)^2 &= \lambda^2\\ +(+1 + D)^2 + (0 - C)^2 &= \lambda^2. +\end{align*} +$$ + + +Squaring out and solving gives $D=0$, $1 + C^2 = \lambda^2$. That is, an arc of circle with radius $\sqrt{1+C^2}$ and centered at $(0, C)$. + + +$$ +x^2 + (y - C)^2 = 1 + C^2. +$$ + +Now to identify $C$ in terms of $L$. $L$ is the length of arc of circle of radius $r =\sqrt{1 + C^2}$ and angle $2\theta$, so $L = 2r\theta$ But using the boundary conditions in the equations for $x$ and $y$ gives $\tan\theta = 1/C$, so $L = 2\sqrt{1 + C^2}\tan^{-1}(1/C)$ which can be solved for $C$ provided $L \geq 2$. + + +##### Example: more constraints + + +Consider now the case of maximizing $f(x,y,z)$ subject to $g(x,y,z)=c$ and $h(x,y,z) = d$. Can something similar be said to characterize potential values for this to occur? Trying to describe where $g(x,y,z) = c$ and $h(x,y,z)=d$ in general will prove difficult. The easy case would be it the two equations were linear, in which case they would describe planes. Two non-parallel planes would intersect in a line. If the general case, imagine the surfaces locally replaced by their tangent planes, then their intersection would be a line, and this line would point in along the curve given by the intersection of the surfaces formed by the constraints. This line is similar to the tangent line in the $2$-variable case. Now if $\nabla{f}$, which points in the direction of greatest increase of $f$, had a non-zero projection onto this line, then moving the point in that direction along the line would increase $f$ and still leave the point following the constraints. That is, if there is a non-zero directional derivative the point is not a maximum. + + +The tangent planes are *orthogonal* to the vectors $\nabla{g}$ and $\nabla{h}$, so in this case parallel to $\nabla{g} \times \nabla{h}$. The condition that $\nabla{f}$ be *orthogonal* to this vector, means that $\nabla{f}$ *must* sit in the plane described by $\nabla{g}$ and $\nabla{h}$ - the plane of orthogonal vectors to $\nabla{g} \times \nabla{h}$. That is, this condition is needed: + + +$$ +\nabla{f}(x,y,z) = \lambda_1 \nabla{g}(x,y,z) + \lambda_2 \nabla{h}(x,y,z). +$$ + +At a point satisfying the above, we would have the tangent "plane" of $f$ is contained in the intersection of the tangent "plane"s to $g$ and $h$. + + +--- + + +Consider a curve given through the intersection of two expressions: $g_1(x,y,z) = x^2 + y^2 - z^2 = 0$ and $g_2(x,y,z) = x - 2z = 3$. What is the minimum distance to the origin along this curve? + + +We have $f(x,y,z) = \text{distance}(\vec{x},\vec{0}) = \sqrt{x^2 + y^2 + z^2}$, subject to the two constraints. As the square root is increasing, we can actually just consider $f(x,y,z) = x^2 + y^2 + z^2$, ignoring the square root. The Lagrange multiplier technique instructs us to look for solutions to: + + +$$ +\langle 2x, 2y ,2z \rangle = \lambda_1\langle 2x, 2y, -2z\rangle + \lambda_2 \langle 1, 0, -2 \rangle. +$$ + +Here we use `SymPy`: + + +```{julia} +@syms z lambda1 lambda2 +g1(x, y, z) = x^2 + y^2 - z^2 +g2(x, y, z) = x - 2z - 3 +fₘ(x,y,z)= x^2 + y^2 + z^2 +Lₘ(x,y,z,lambda1, lambda2) = fₘ(x,y,z) - lambda1*(g1(x,y,z) - 0) - lambda2*(g2(x,y,z) - 0) + +∇Lₘ = diff.(Lₘ(x,y,z,lambda1, lambda2), [x, y, z,lambda1, lambda2]) +``` + +Before trying to solve for $\nabla{L} = \vec{0}$ we see from the second equation that *either* $\lambda_1 = 1$ or $y = 0$. First we solve with $\lambda_1 = 1$: + + +```{julia} +solve(subs.(∇Lₘ, lambda1 .=> 1)) +``` + +There are no real solutions. Next when $y = 0$ we get: + + +```{julia} +outₘ = solve(subs.(∇Lₘ, y .=> 0)) +``` + +The two solutions have values yielding the extrema: + + +```{julia} +[fₘ(d[x], 0, d[z]) for d in outₘ] +``` + +## Taylor's theorem + + +Taylor's theorem for a univariate function states that if $f$ has $k+1$ derivatives in an open interval around $a$, $f^{(k)}$ is continuous between the closed interval from $a$ to $x$ then: + + +$$ +f(x) = \sum_{j=0}^k \frac{f^{j}(a)}{j!} (x-a)^j + R_k(x), +$$ + +where $R_k(x) = f^{k+1}(\xi)/(k+1)!(x-a)^{k+1}$ for some $\xi$ between $a$ and $x$. + + +This theorem can be generalized to scalar functions, but the notation can be cumbersome. Following [Folland](https://sites.math.washington.edu/~folland/Math425/taylor2.pdf) we use *multi-index* notation. Suppose $f:R^n \rightarrow R$, and let $\alpha=(\alpha_1, \alpha_2, \dots, \alpha_n)$. Then define the following notation: + + +$$ +\begin{align*} +|\alpha| &= \alpha_1 + \cdots + \alpha_n, \\ +\alpha! &= \alpha_1!\alpha_2!\cdot\cdots\cdot\alpha_n!, \\ +\vec{x}^\alpha &= x_1^{\alpha_1}x_2^{\alpha_2}\cdots x_n^{\alpha^n}, \\ +\partial^\alpha f &= \partial_1^{\alpha_1}\partial_2^{\alpha_2}\cdots \partial_n^{\alpha_n} f \\ +& = \frac{\partial^{|\alpha|}f}{\partial x_1^{\alpha_1} \partial x_2^{\alpha_2} \cdots \partial x_n^{\alpha_n}}. +\end{align*} +$$ + + +This notation makes many formulas from one dimension carry over to higher dimensions. For example, the binomial theorem says: + + +$$ +(a+b)^n = \sum_{k=0}^n \frac{n!}{k!(n-k)!}a^kb^{n-k}, +$$ + +and this becomes: + + +$$ +(x_1 + x_2 + \cdots + x_n)^n = \sum_{|\alpha|=k} \frac{k!}{\alpha!} \vec{x}^\alpha. +$$ + +::: {.callout-note icon=false} +## Taylor's theorem using multi-index + +If $f: R^n \rightarrow R$ is sufficiently smooth ($C^{k+1}$) on an open convex set $S$ about $\vec{a}$ then if $\vec{a}$ and $\vec{a}+\vec{h}$ are in $S$, + + +$$ +f(\vec{a} + \vec{h}) = \sum_{|\alpha| \leq k}\frac{\partial^\alpha f(\vec{a})}{\alpha!}\vec{h}^\alpha + R_{\vec{a},k}(\vec{h}), +$$ + +where $R_{\vec{a},k} = \sum_{|\alpha|=k+1}\partial^\alpha \frac{f(\vec{a} + c\vec{h})}{\alpha!} \vec{h}^\alpha$ for some $c$ in $(0,1)$. + +::: + +##### Example + + +The elegant notation masks what can be complicated expressions. Consider the simple case $f:R^2 \rightarrow R$ and $k=2$. Then this says: + +$$ +\begin{align*} +f(x + dx, y+dy) &= f(x, y) + \frac{\partial f}{\partial x} dx + \frac{\partial f}{\partial y} dy \\ +&+ \frac{\partial^2 f}{\partial x^2} \frac{dx^2}{2} + 2\frac{\partial^2 f}{\partial x\partial y} \frac{dx dy}{2}\\ +&+ \frac{\partial^2 f}{\partial y^2} \frac{dy^2}{2} + R_{\langle x, y \rangle, k}(\langle dx, dy \rangle). +\end{align*} +$$ + +Using $\nabla$ and $H$ for the Hessian and $\vec{x} = \langle x, y \rangle$ and $d\vec{x} = \langle dx, dy \rangle$, this can be expressed as: + + +$$ +f(\vec{x} + d\vec{x}) = f(\vec{x}) + \nabla{f} \cdot d\vec{x} + d\vec{x} \cdot (H d\vec{x}) +R_{\vec{x}, k}d\vec{x}. +$$ + +As for $R$, the full term involves terms for $\alpha = (3,0), (2,1), (1,2)$, and $(0,3)$. Using $\vec{a} = \langle x, y\rangle$ and $\vec{h}=\langle dx, dy\rangle$: + + +$$ +\frac{\partial^3 f(\vec{a}+c\vec{h})}{\partial x^3} \frac{dx^3}{3!}+ +\frac{\partial^3 f(\vec{a}+c\vec{h})}{\partial x^2\partial y} \frac{dx^2 dy}{2!1!} + +\frac{\partial^3 f(\vec{a}+c\vec{h})}{\partial x\partial y^2} \frac{dxdy^2}{1!2!} + +\frac{\partial^3 f(\vec{a}+c\vec{h})}{\partial y^3} \frac{dy^3}{3!}. +$$ + +The exact answer is usually not as useful as the bound: $|R| \leq M/(k+1)! \|\vec{h}\|^{k+1}$, for some finite constant $M$. + + +##### Example + + +We can encode multiindices using `SymPy`. The basic definitions are fairly straightforward using `zip` to pair variables with components of $\alpha$. We define a new type so that we can overload the familiar notation: + + +```{julia} +struct MultiIndex + alpha::Vector{Int} + end +Base.show(io::IO, α::MultiIndex) = println(io, "α = ($(join(α.alpha, ", ")))") + +## |α| = α_1 + ... + α_m +Base.length(α::MultiIndex) = sum(α.alpha) + +## factorial(α) computes α! +Base.factorial(α::MultiIndex) = prod(factorial(Sym(a)) for a in α.alpha) + +## x^α = x_1^α_1 * x_2^α^2 * ... * x_n^α_n +import Base: ^ +^(x, α::MultiIndex) = prod(u^a for (u,a) in zip(x, α.alpha)) + +## ∂^α(ex) = ∂_1^α_1 ∘ ∂_2^α_2 ∘ ... ∘ ∂_n^α_n (ex) +partial(ex::SymPy.SymbolicObject, α::MultiIndex, vars=free_symbols(ex)) = diff(ex, zip(vars, α.alpha)...) +``` + +```{julia} +@syms w +alpha = MultiIndex([1,2,1,3]) +length(alpha) # 1 + 2 + 1 + 3=7 +[1,2,3,4]^alpha +exₜ = x^3 * cos(w*y*z) +partial(exₜ, alpha, [w,x,y,z]) +``` + +The remainder term needs to know information about sets like $|\alpha| =k$. This is a combinatoric problem, even to identify the length. Here we define an iterator to iterate over all possible MultiIndexes. This is low level, and likely could be done in a much better style, so shouldn't be parsed unless there is curiosity. It manually chains together iterators. + + +```{julia} +struct MultiIndices + n::Int + k::Int +end + +function Base.length(as::MultiIndices) + n,k = as.n, as.k + n == 1 && return 1 + sum(length(MultiIndices(n-1, j)) for j in 0:k) # recursively identify length +end + +function Base.iterate(alphas::MultiIndices) + k, n = alphas.k, alphas.n + n == 1 && return ([k],(0, MultiIndices(0,0), nothing)) + + m = zeros(Int, n) + m[1] = k + betas = MultiIndices(n-1, 0) + stb = iterate(betas) + st = (k, MultiIndices(n-1, 0), stb) + return (m, st) +end + +function Base.iterate(alphas::MultiIndices, st) + + st == nothing && return nothing + k,n = alphas.k, alphas.n + k == 0 && return nothing + n == 1 && return nothing + + # can we iterate the next on + bk, bs, stb = st + + if stb==nothing + bk = bk-1 + bk < 0 && return nothing + bs = MultiIndices(bs.n, bs.k+1) + val, stb = iterate(bs) + return (vcat(bk,val), (bk, bs, stb)) + end + + resp = iterate(bs, stb) + if resp == nothing + bk = bk-1 + bk < 0 && return nothing + bs = MultiIndices(bs.n, bs.k+1) + val, stb = iterate(bs) + return (vcat(bk, val), (bk, bs, stb)) + end + + val, stb = resp + return (vcat(bk, val), (bk, bs, stb)) + +end +``` + +This returns a vector, not a `MultiIndex`. Here we get all multiindices in two variables of size $3$ + + +```{julia} +collect(MultiIndices(2, 3)) +``` + +To get all of size $3$ or less, we could do something like this: + + +```{julia} +union((collect(MultiIndices(2, i)) for i in 0:3)...) +``` + +To see the computational complexity. Suppose we had $3$ variables and were interested in the error for order $4$: + + +```{julia} +k = 4 +length(MultiIndices(3, k+1)) +``` + +Finally, to see how compact the notation issue, suppose $f:R^3 \rightarrow R$, we have the third-order Taylor series expands to $20$ terms as follows: + + +```{julia} +#| hold: true +@syms 𝐅() a[1:3] dx[1:3] + +sum(partial(𝐅(a...), α, a) / factorial(α) * dx^α for k in 0:3 for α in MultiIndex.(MultiIndices(3, k))) # 3rd order +``` + +## Questions + + +###### Question + + +Let $f(x,y) = \sqrt{x + y}$. Find the tangent plane approximation for $f(2.1, 2.2)$? + + +```{julia} +#| hold: true +#| echo: false +f(x,y) = sqrt(x + y) +f(v) = f(v...) +pt = [2,2] +dxdy = [.1, .2] +val = f(pt) + dot(ForwardDiff.gradient(f, pt), dxdy) +numericq(val) +``` + +###### Question + + +Let $f(x,y,z) = xy + yz + zx$. Using a *linear approximation* estimate $f(1.1, 1.0, 0.9)$. + + +```{julia} +#| hold: true +#| echo: false +f(x,y,z) = x*y + y*z + z*x +f(v) = f(v...) +pt = [1,1,1] +dx = [0.1, 0.0, -0.1] +val = f(pt) + ∇(f)(pt) ⋅ dx +numericq(val) +``` + +###### Question + + +Let $f(x,y,z) = xy + yz + zx - 3$. What equation describes the tangent approximation at $(1,1,1)$? + + +```{julia} +#| hold: true +#| echo: false +f(x,y,z) = x*y + y*z + z*x - 8 +f(v) = f(v...) +pt = [1,1,1] +n = ∇(f)(pt) +d = dot(n, pt) +choices = [ + raw"`` x + y + z = 3``", + raw"`` 2x + y - 2z = 1``", + raw"`` x + 2y + 3z = 6``" +] +answ = 1 +radioq(choices, answ) +``` + +###### Question + + +([Knill](http://www.math.harvard.edu/~knill/teaching/summer2018/handouts/week4.pdf)) Let $f(x,y) = xy + x^2y + xy^2$. + + +Find the gradient of $f$: + + +```{julia} +#| hold: true +#| echo: false +choices = [ + raw"`` \langle 2xy + y^2 + y, 2xy + x^2 + x\rangle``", + raw"`` y^2 + y, x^2 + x``", + raw"`` \langle 2y + y^2, 2x + x^2``" +] +answ = 1 +radioq(choices, answ) +``` + +Is this the Hessian of $f$? + + +$$ +\left[\begin{matrix}2 y & 2 x + 2 y + 1\\2 x + 2 y + 1 & 2 x\end{matrix}\right] +$$ + +```{julia} +#| hold: true +#| echo: false +yesnoq(true) +``` + +The point $(-1/3, -1/3)$ is a solution to the $\nabla{f} = 0$. What is the *determinant*, $d$, of the Hessian at this point? + + +```{julia} +#| hold: true +#| echo: false +f(x,y) = x*y + x*y^2 + x^2 * y +f(v) = f(v...) +val = det(ForwardDiff.hessian(f, [-1/3, -1/3])) +numericq(val) +``` + +Which is true of $f$ at $(-1/3, -1/3)$: + + +```{julia} +#| hold: true +#| echo: false +choices = [ + L"The function $f$ has a local minimum, as $f_{xx} > 0$ and $d >0$", + L"The function $f$ has a local maximum, as $f_{xx} < 0$ and $d >0$", + L"The function $f$ has a saddle point, as $d < 0$", + L"Nothing can be said, as $d=0$" +] +answ = 2 +radioq(choices, answ, keep_order=true) +``` + +###### Question + + +([Knill](http://www.math.harvard.edu/~knill/teaching/summer2018/handouts/week4.pdf)) Let the Tutte polynomial be $f(x,y) = x + 2x^2 + x^3 + y + 2xy + y^2$. + + +Does this accurately find the gradient of $f$? + + +```{julia} +#| hold: true +#| results: "hidden" +f(x,y) = x + 2x^2 + x^3 + y + 2x*y + y^2 +@syms x::real y::real +gradf = gradient(f(x,y), [x,y]) +``` + +```{julia} +#| hold: true +#| echo: false +yesnoq(true) +``` + +How many answers does this find to $\nabla{f} = \vec{0}$? + + +```{julia} +#| hold: true +#| results: "hidden" +f(x,y) = x + 2x^2 + x^3 + y + 2x*y + y^2 +@syms x::real y::real +gradf = gradient(f(x,y), [x,y]) + +solve(gradf, [x,y]) +``` + +```{julia} +#| hold: true +#| echo: false +numericq(2) +``` + +The Hessian is found by + + +```{julia} +#| hold: true +f(x,y) = x + 2x^2 + x^3 + y + 2x*y + y^2 +@syms x::real y::real +gradf = gradient(f(x,y), [x,y]) + +sympy.hessian(f(x,y), [x,y]) +``` + +Which is true of $f$ at $(-2/3, 1/6)$: + + +```{julia} +#| hold: true +#| echo: false +choices = [ + L"The function $f$ has a local minimum, as $f_{xx} > 0$ and $d >0$", + L"The function $f$ has a local maximum, as $f_{xx} < 0$ and $d >0$", + L"The function $f$ has a saddle point, as $d < 0$", + L"Nothing can be said, as $d=0$", + L"The test does not apply, as $\nabla{f}$ is not $0$ at this point." +] +answ = 3 +radioq(choices, answ, keep_order=true) +``` + +Which is true of $f$ at $(0, -1/2)$: + + +```{julia} +#| hold: true +#| echo: false +choices = [ + L"The function $f$ has a local minimum, as $f_{xx} > 0$ and $d >0$", + L"The function $f$ has a local maximum, as $f_{xx} < 0$ and $d >0$", + L"The function $f$ has a saddle point, as $d < 0$", + L"Nothing can be said, as $d=0$", + L"The test does not apply, as $\nabla{f}$ is not $0$ at this point." +] +answ = 1 +radioq(choices, answ, keep_order=true) +``` + +Which is true of $f$ at $(1/2, 0)$: + + +```{julia} +#| hold: true +#| echo: false +choices = [ + L"The function $f$ has a local minimum, as $f_{xx} > 0$ and $d >0$", + L"The function $f$ has a local maximum, as $f_{xx} < 0$ and $d >0$", + L"The function $f$ has a saddle point, as $d < 0$", + L"Nothing can be said, as $d=0$", + L"The test does not apply, as $\nabla{f}$ is not $0$ at this point." +] +answ = 5 +radioq(choices, answ, keep_order=true) +``` + +###### Question + + +(Strang p509) Consider the quadratic function $f(x,y) = ax^2 + bxy +cy^2$. Since the second partial derivative test is essentially done by replacing the function at a critical point by a quadratic function, understanding this $f$ is of some interest. + + +Is this the Hessian of $f$? + + +$$ +\begin{bmatrix} +2a & b\\ +b & 2c +\end{bmatrix} +$$ + +```{julia} +#| hold: true +#| echo: false +yesnoq(true) +``` + +Or is this the Hessian of $f$? + + +$$ +\begin{bmatrix} +2ax & by\\ +bx & 2cy +\end{bmatrix} +$$ + +```{julia} +#| hold: true +#| echo: false +yesnoq(false) +``` + +Explain why $4ac - b^2$ is of any interest here: + + +```{julia} +#| hold: true +#| echo: false +choices =[ + "It is the determinant of the Hessian", + L"It isn't, $b^2-4ac$ is from the quadratic formula" +] +answ = 1 +radioq(choices, answ) +``` + +Which condition on $a$, $b$, and $c$ will ensure a *local maximum*: + + +```{julia} +#| hold: true +#| echo: false +choices = [ + L"That $a>0$ and $4ac-b^2 > 0$", + L"That $a<0$ and $4ac-b^2 > 0$", + L"That $4ac-b^2 < 0$" +] +answ = 2 +radioq(choices, answ, keep_order=true) +``` + +Which condition on $a$, $b$, and $c$ will ensure a saddle point? + + +```{julia} +#| hold: true +#| echo: false +choices = [ + L"That $a>0$ and $4ac-b^2 > 0$", + L"That $a<0$ and $4ac-b^2 > 0$", + L"That $4ac-b^2 < 0$" +] +answ = 3 +radioq(choices, answ, keep_order=true) +``` + +###### Question + + +Let $f(x,y) = e^{-x^2 - y^2} (2x^2 + y^2)$. Use Lagrange's method to find the absolute maximum and absolute minimum over $x^2 + y^2 = 3$. + + +Is $\nabla{f}$ given by the following? + + +$$ +\nabla{f} =2 e^{-x^2 - y^2} \langle x(2 - 2x^2 - y^2), y(1 - 2x^2 - y^2)\rangle. +$$ + +```{julia} +#| hold: true +#| echo: false +yesnoq(true) +``` + +Which vector is orthogonal to the contour line $x^2 + y^2 = 3$? + + +```{julia} +#| echo: false +choices = [ + raw"`` \langle 2x, 2y\rangle``", + raw"`` \langle 2x, y^2\rangle``", + raw"`` \langle x^2, 2y \rangle``" +] +answ = 1 +radioq(choices, answ) +``` + +Due to the form of the gradient of the constraint, finding when $\nabla{f} = \lambda \nabla{g}$ is the same as identifying when this ratio $|f_x/f_y|$ is $1$. The following solves for this by checking each point on the constraint: + + +```{julia} +f(x,y) = exp(-x^2-y^2) * (2x^2 + y^2) +f(v) = f(v...) +r(t) = sqrt(3)*[cos(t), sin(t)] +rat(x) = abs(x[1]/x[2]) - 1 +fn = rat ∘ ∇(f) ∘ r +ts = fzeros(fn, 0, 2pi) +``` + +Using these points, what is the largest value on the boundary? + + +```{julia} +#| eval: false +#| echo: false +f(x,y) = exp(-x^2-y^2) * (2x^2 + y^2) +r(t) = sqrt(3)*[cos(t), sin(t)] +rat(x) = abs(x[1]/x[2]) - 1 +fn = rat ∘ ∇(splat(f)) ∘ r +ts = fzeros(fn, 0, 2pi) + +val = maximum((splat(u)∘r).(ts)) +numericq(val) +``` diff --git a/quarto/differentiable_vector_calculus/scalar_functions_applications.qmd b/quarto/differentiable_vector_calculus/scalar_functions_applications.qmd index 2b74579..dec255f 100644 --- a/quarto/differentiable_vector_calculus/scalar_functions_applications.qmd +++ b/quarto/differentiable_vector_calculus/scalar_functions_applications.qmd @@ -9,7 +9,7 @@ This section uses these add-on packages: ```{julia} using CalculusWithJulia using Plots -plotly() +#plotly() using SymPy using Roots ``` @@ -33,7 +33,6 @@ Consider the case $f:R^2 \rightarrow R$. We visualize $z=f(x,y)$ through a surfa For the univariate case, the tangent line has many different uses. Here we see the tangent plane also does. - ### Equation of the tangent plane @@ -119,6 +118,7 @@ arrow!(pt, [1,0,0], linestyle=:dash) arrow!(pt, [0,1,0], linestyle=:dash) ``` + #### Alternate forms @@ -170,7 +170,6 @@ 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 @@ -231,6 +230,7 @@ t = 1/2 (n1(gamma(t)) × n2(gamma(t))) × gamma'(t) ``` + #### Plotting level curves of $F(x,y,z) = c$ @@ -250,7 +250,6 @@ 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) ``` - ## Linearization @@ -329,7 +328,7 @@ Here, $\nabla{f}$ describes a *normal* to the tangent plane. The description of 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])" +println("$a x + $b y + $c z = $([a,b,c] ⋅ [2,2,2])"); ``` ### Newton's method to solve $f(x,y) = 0$ and $g(x,y)=0$. @@ -504,6 +503,9 @@ function nm(f, g, x, n=5) end ``` + + + ##### Example @@ -655,6 +657,7 @@ plot(F, 0.01, 5) # p > 0 This problem does not have a readily expressed value for $x^*$, but when $p \approx 0$ we should get similar behavior to the intersection of $y=px$ and $y=\pi/2 - x$ for $x^*$, or $x^* \approx \pi/(2(1+p))$ which has derivative of $-\pi/2$ at $p=0$, matching the above graph. For *large* $p$, the problem looks like the intersection of the line $y=1$ with $y=px$ or $x^* \approx 1/p$ which has derivative that goes to $0$ as $p$ goes to infinity, again matching this graph. + ## Optimization @@ -807,7 +810,6 @@ det(H_a) (The test is inconclusive, as it needs the function to "fall away" from the tangent plane in all directions, in this case, along a circular curve, the function touches the tangent plane, so it doesn't fall away.) - ##### Example @@ -853,6 +855,7 @@ end p ``` + ##### Example @@ -875,9 +878,8 @@ At $\vec{a}$ this has positive determinant and $f_{xx} > 0$, so $\vec{a}$ corres ```{julia} fₗ(x,y) = x^2 + 2y^2 - x -fₗ(v) = fₗ(v...) gammaₗ(t) = [cos(t), sin(t)] # traces out x^2 + y^2 = 1 over [0, 2pi] -gₗ = fₗ ∘ gammaₗ +gₗ = splat(fₗ) ∘ gammaₗ cpsₗ = find_zeros(gₗ', 0, 2pi) # critical points of g append!(cpsₗ, [0, 2pi]) @@ -907,6 +909,12 @@ So we have the maximum occurs at the angles $2\pi/3$ and $4\pi/3$. Here we visua hₗ(x,y) = fₗ(x,y) * (x^2 + y^2 <= 1 ? 1 : NaN) ``` + +```{julia} +#| echo: false +gr(); # scatter3d! and plotly() are failing!! +``` + ```{julia} #| hold: true xs = ys = range(-1,1, length=100) @@ -918,6 +926,11 @@ zs = fₗ.(xs, ys) scatter3d!(xs, ys, zs) ``` +```{julia} +#| echo: false +#plotly(); +``` + A contour plot also shows that some---and only one---extrema happens on the interior: @@ -930,6 +943,7 @@ contour(xs, ys, hₗ) The extrema are identified by the enclosing regions, in this case the one around the point $(1/2, 0)$. + ##### Example: Steiner's problem @@ -1189,6 +1203,10 @@ end We now visualize, using the `Contour` package to draw the contour lines in the $x-y$ plane: +```{julia} +#| echo: false +gr() +``` ```{julia} #| hold: true @@ -1217,6 +1235,11 @@ scatter3d!(unzip(pts)...) plot!(unzip(pts)..., linewidth=3) ``` +```{julia} +#| echo: false +#plotly() +``` + ### Newton's method for minimization @@ -2324,7 +2347,6 @@ Which vector is orthogonal to the contour line $x^2 + y^2 = 3$? ```{julia} -#| hold: true #| echo: false choices = [ raw"`` \langle 2x, 2y\rangle``", @@ -2339,7 +2361,6 @@ Due to the form of the gradient of the constraint, finding when $\nabla{f} = \la ```{julia} -#| hold: true f(x,y) = exp(-x^2-y^2) * (2x^2 + y^2) f(v) = f(v...) r(t) = sqrt(3)*[cos(t), sin(t)] @@ -2352,15 +2373,14 @@ Using these points, what is the largest value on the boundary? ```{julia} -#| hold: true +#| eval: false #| echo: false f(x,y) = exp(-x^2-y^2) * (2x^2 + y^2) -f(v) = f(v...) r(t) = sqrt(3)*[cos(t), sin(t)] rat(x) = abs(x[1]/x[2]) - 1 -fn = rat ∘ ∇(f) ∘ r +fn = rat ∘ ∇(splat(f)) ∘ r ts = fzeros(fn, 0, 2pi) -val = maximum((f∘r).(ts)) +val = maximum((splat(u)∘r).(ts)) numericq(val) ``` diff --git a/quarto/differentiable_vector_calculus/test.html b/quarto/differentiable_vector_calculus/test.html new file mode 100644 index 0000000..26d9d85 --- /dev/null +++ b/quarto/differentiable_vector_calculus/test.html @@ -0,0 +1,641 @@ + + + + + + + + + +test + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ +
+ +
+
+

Applications with scalar functions

+
+ + + +
+ + + + +
+ + + +
+ + +

This section uses these add-on packages:

+
+
using CalculusWithJulia
+using Plots
+plotly()
+using SymPy
+using Roots
+
+
+
Example
+

Consider the function \(f(x,y) = x^2 + 3y^2 -x\) over the region \(x^2 + y^2 \leq 1\). This is a continuous function over a closed set, so will have both an absolute maximum and minimum. Find these from an investigation of the critical points and the boundary points.

+

The gradient is easily found: \(\nabla{f} = \langle 2x - 1, 6y \rangle\), and is \(\vec{0}\) only at \(\vec{a} = \langle 1/2, 0 \rangle\). The Hessian is:

+

\[ +H = +\begin{bmatrix} +2 & 0\\ +0 & 6 +\end{bmatrix}. +\]

+

At \(\vec{a}\) this has positive determinant and \(f_{xx} > 0\), so \(\vec{a}\) corresponds to a local minimum with values \(f(\vec{a}) = (1/2)^2 + 3(0) - 1/2 = -1/4\). The absolute maximum and minimum may occur here (well, not the maximum) or on the boundary, so that must be considered. In this case we can easily parameterize the boundary and turn this into the univariate case:

+
+
fₗ(x,y) = x^2 + 2y^2 - x
+gammaₗ(t) = [cos(t), sin(t)]  # traces out x^2 + y^2 = 1 over [0, 2pi]
+gₗ = splat(fₗ) ∘ gammaₗ
+
+cpsₗ = find_zeros(gₗ', 0, 2pi) # critical points of g
+append!(cpsₗ, [0, 2pi])
+unique!(cpsₗ)
+gₗ.(cpsₗ)
+
+
5-element Vector{Float64}:
+ 0.0
+ 2.25
+ 2.0
+ 2.25
+ 0.0
+
+
+

We see that maximum value is 2.25 and that the interior point, \(\vec{a}\), will be where the minimum value occurs. To see exactly where the maximum occurs, we look at the values of gamma:

+
+
inds = [2,4]
+cpsₗ[inds]
+
+
2-element Vector{Float64}:
+ 2.0943951023931953
+ 4.1887902047863905
+
+
+

These are multiples of \(\pi\):

+
+
cpsₗ[inds]/pi
+
+
2-element Vector{Float64}:
+ 0.6666666666666666
+ 1.3333333333333333
+
+
+

So we have the maximum occurs at the angles \(2\pi/3\) and \(4\pi/3\). Here we visualize, using a hacky trick of assigning NaN values to the function to avoid plotting outside the circle:

+
+
hₗ(x,y) = fₗ(x,y) * (x^2 + y^2 <= 1 ? 1 : NaN)
+
+
hₗ (generic function with 1 method)
+
+
+
+
gr()
+
+
Plots.GRBackend()
+
+
+
+
xs = ys = range(-1,1, length=100)
+plt = surface(xs, ys, hₗ)
+
+ts = cpsₗ  # 2pi/3 and 4pi/3 by above
+xs, ys = cos.(ts), sin.(ts)
+zs = fₗ.(xs, ys)
+scatter3d!(xs, ys, zs)
+#tuple.(xs, ys, zs)
+#plot!(plt, tuple.(xs, ys, zs); linetype=:scatter)
+#plt
+
+ +
+
+

A contour plot also shows that some—and only one—extrema happens on the interior:

+ +
+ +
+ + +
+ + + + + \ No newline at end of file diff --git a/quarto/differentiable_vector_calculus/test.jl b/quarto/differentiable_vector_calculus/test.jl new file mode 100644 index 0000000..fe58175 --- /dev/null +++ b/quarto/differentiable_vector_calculus/test.jl @@ -0,0 +1,827 @@ +@show 4 +using QuizQuestions +using LaTeXStrings +using CalculusWithJulia +using Plots +plotly() +using SymPy +using Roots +@show 6 +import Contour: contours, levels, level, lines, coordinates +@show 15 +@syms f_x f_y +n = [1, 0, f_x] × [0, 1, f_y] +@show 27 +#| hold: true +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) +@show 33 +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 +@show 35 +function tangent_plane(f, pt) + ∇f = ForwardDiff.gradient(f, pt) # using a variable ∇f + x -> f(pt) + ∇f ⋅ (x - pt) +end +@show 46 +@syms x, y +@show 47 +#| hold: true +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) +@show 55 +#| hold: true +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) +@show 60 +#| hold: true +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) +@show 71 +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) +@show 73 +V(a₁ + dx₁) - V(a₁) +@show 85 +#| hold: true +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])" +#@show 92 +#| hold: true +@syms a b c d u v +M = [a b; c d] +B = [u, v] +M \ B .|> simplify +@show 96 +#| hold: true +#| echo: false +f(x,y) = 2 - x^2 - y^2 +g(x,y) = 3 - 2x^2 - (1/3)y^2 +xs = ys = range(-3, stop=3, length=100) +zfs = [f(x,y) for x in xs, y in ys] +zgs = [g(x,y) for x in xs, y in ys] + + +ps = Any[] +pf = surface(xs, ys, f, alpha=0.5, legend=false) + +for cl in levels(contours(xs, ys, zfs, [0.0])) + for line in lines(cl) + _xs, _ys = coordinates(line) + plot!(pf, _xs, _ys, 0*_xs, linewidth=3, color=:blue) + end +end + + +pg = surface(xs, ys, g, alpha=0.5, legend=false) +for cl in levels(contours(xs, ys, zgs, [0.0])) + for line in lines(cl) + _xs, _ys = coordinates(line) + plot!(pg, _xs, _ys, 0*_xs, linewidth=3, color=:red) + end +end + +pcnt = plot(legend=false) +for cl in levels(contours(xs, ys, zfs, [0.0])) + for line in lines(cl) + _xs, _ys = coordinates(line) + plot!(pcnt, _xs, _ys, linewidth=3, color=:blue) + end +end + +for cl in levels(contours(xs, ys, zgs, [0.0])) + for line in lines(cl) + _xs, _ys = coordinates(line) + plot!(pcnt, _xs, _ys, linewidth=3, color=:red) + end +end + +l = @layout([a b c]) +plot(pf, pg, pcnt, layout=l) +@show 106 +function newton_step(f, g, xn) + M = [ForwardDiff.gradient(f, xn)'; ForwardDiff.gradient(g, xn)'] + b = -[f(xn), g(xn)] + Delta = M \ b + xn + Delta +end +@show 108 +𝒇(x,y) = 2 - x^2 - y^2 +𝒈(x,y) = 3 - 2x^2 - (1/3)y^2 +𝒇(v) = 𝒇(v...); 𝒈(v) = 𝒈(v...) +𝒙₀ = [1,1] +𝒙₁ = newton_step(𝒇, 𝒈, 𝒙₀) +@show 110 +𝒇(𝒙₁), 𝒈(𝒙₁) +@show 112 +𝒙₂ = newton_step(𝒇, 𝒈, 𝒙₁) +𝒙₃ = newton_step(𝒇, 𝒈, 𝒙₂) +𝒙₄ = newton_step(𝒇, 𝒈, 𝒙₃) +𝒙₅ = newton_step(𝒇, 𝒈, 𝒙₄) +𝒙₅, 𝒇(𝒙₅), 𝒈(𝒙₅) +@show 116 +function nm(f, g, x, n=5) + for i in 1:n + x = newton_step(f, g, x) + end + x +end +@show 123 +#| hold: true +c = 1/2 +f(x,y) = 1 - y^2 - c^2 +g(x,y) = (1 - x^2) - c^2 +f(v) = f(v...); g(v) = g(v...) +nm(f, g, [1/2, 1/3]) +@show 148 +#| hold: true +@syms x, y, Z() +∂x = solve(diff(x^4 -x^3 + y^2 + Z(x,y)^2, x), diff(Z(x,y),x)) +∂y = solve(diff(x^4 -x^3 + y^2 + Z(x,y)^2, y), diff(Z(x,y),y)) +∂x, ∂y +@show 158 +f(x, p) = cos(x) - p*x +p = 2 +xᵅ = find_zero(f, (0, pi/2), p) +@show 160 +p = 2 +xᵅ = find_zero(f, (0, pi/2), p) +fₓ = ForwardDiff.derivative(x -> f(x,p), xᵅ) +fₚ = ForwardDiff.derivative(p -> f(xᵅ, p), p) +- fₚ / fₓ +@show 163 +function find_zero_derivative(f, x₀, p) + xᵅ = find_zero(f, x₀, p) + fₓ = ForwardDiff.derivative(x -> f(x,p), xᵅ) + fₚ = ForwardDiff.derivative(p -> f(xᵅ, p), p) + - fₚ / fₓ +end +F(p) = find_zero_derivative(f, (0, pi/2), p) +plot(F, 0.01, 5) # p > 0 +@show 183 +#| hold: true +f(x,y)= exp(-(x^2 + y^2)/5) * cos(x^2 + y^2) +xs = ys = range(-4, 4, length=100) +surface(xs, ys, f, legend=false) +@show 190 +#| hold: true +f(x,y) = x*y +xs = ys = range(-3, 3, length=100) +surface(xs, ys, f, legend=false) + +plot_parametric!(-4..4, t -> [t, 0, f(t, 0)], linewidth=5) +plot_parametric!(-4..4, t -> [0, t, f(0, t)], linewidth=5) +@show 203 +fₖ(x,y) = exp(-(x^2 + y^2)/5) * cos(x^2 + y^2) +Hₖ = sympy.hessian(fₖ(x,y), (x,y)) +@show 205 +H₀₀ = subs.(Hₖ, x=>0, y=>0) +@show 207 +H₀₀[1,1] < 0 && det(H₀₀) > 0 +@show 209 +#| hold: true +gradfₖ = diff.(fₖ(x,y), [x,y]) +a = [sqrt(2PI + atan(-Sym(1)//5)), 0] +subs.(gradfₖ, x => a[1], y => a[2]) +@show 211 +#| hold: true +a = [sqrt(PI + atan(-Sym(1)//5)), 0] +H_a = subs.(Hₖ, x => a[1], y => a[2]) +det(H_a) +@show 216 +fⱼ(x,y) = 4x*y - x^4 - y^4 +gradfⱼ = diff.(fⱼ(x,y), [x,y]) +@show 217 +all_ptsⱼ = solve(gradfⱼ, [x,y]) +ptsⱼ = filter(u -> all(isreal.(u)), all_ptsⱼ) +@show 219 +Hⱼ = sympy.hessian(fⱼ(x,y), (x,y)) +function classify(H, pt) + Ha = subs.(H, x => pt[1], y => pt[2]) + (det=det(Ha), f_xx=Ha[1,1]) +end +[classify(Hⱼ, pt) for pt in ptsⱼ] +@show 221 +#| hold: true +xs = ys = range(-3/2, 3/2, length=100) +p = surface(xs, ys, fⱼ, legend=false) +for pt ∈ ptsⱼ + scatter!(p, unzip([N.([pt...,fⱼ(pt...)])])..., + markercolor=:black, markersize=5) # add each pt on surface +end +p +@show 228 +fₗ(x,y) = x^2 + 2y^2 - x +fₗ(v) = fₗ(v...) +gammaₗ(t) = [cos(t), sin(t)] # traces out x^2 + y^2 = 1 over [0, 2pi] +gₗ = fₗ ∘ gammaₗ + +cpsₗ = find_zeros(gₗ', 0, 2pi) # critical points of g +append!(cpsₗ, [0, 2pi]) +unique!(cpsₗ) +gₗ.(cpsₗ) +@show 230 +inds = [2,4] +cpsₗ[inds] +@show 232 +cpsₗ[inds]/pi +@show 234 +hₗ(x,y) = fₗ(x,y) * (x^2 + y^2 <= 1 ? 1 : NaN) +@show 235 +#| hold: true +xs = ys = range(-1,1, length=100) +surface(xs, ys, hₗ) + +ts = cpsₗ # 2pi/3 and 4pi/3 by above +xs, ys = cos.(ts), sin.(ts) +zs = fₗ.(xs, ys) +scatter3d!(xs, ys, zs) +@show 237 +#| hold: true +xs = ys = range(-1,1, length=100) +contour(xs, ys, hₗ) +@show 243 +@syms x1 y1 x2 y2 x3 y3 +d2(p,x) = (p[1] - x[1])^2 + (p[2]-x[2])^2 +d2_1, d2_2, d2_3 = d2((x,y), (x1, y1)), d2((x,y), (x2, y2)), d2((x,y), (x3, y3)) +exₛ = d2_1 + d2_2 + d2_3 +@show 245 +gradfₛ = diff.(exₛ, [x,y]) +xstarₛ = solve(gradfₛ, [x,y]) +@show 248 +Hₛ = subs.(hessian(exₛ, [x,y]), x=>xstarₛ[x], y=>xstarₛ[y]) +@show 259 +usₛ = [[cos(t), sin(t)] for t in (0, 2pi/3, 4pi/3)] +polygon(ps) = unzip(vcat(ps, ps[1:1])) # easier way to plot a polygon + +pₛ = scatter([0],[0], markersize=2, legend=false, aspect_ratio=:equal) + +asₛ = (1,2,3) +plot!(polygon([a*u for (a,u) in zip(asₛ, usₛ)])...) +[arrow!([0,0], a*u, alpha=0.5) for (a,u) in zip(asₛ, usₛ)] +pₛ +@show 261 +asₛ₁ = (1, -1, 3) +scatter([0],[0], markersize=2, legend=false) +psₛₗ = [a*u for (a,u) in zip(asₛ₁, usₛ)] +plot!(polygon(psₛₗ)...) +@show 263 +euclid_dist(x; ps=psₛₗ) = sum(norm(x-p) for p in ps) +euclid_dist(x,y; ps=psₛₗ) = euclid_dist([x,y]; ps=ps) +@show 264 +#| hold: true +xs = range(-1.5, 1.5, length=100) +ys = range(-3, 1.0, length=100) + +p = plot(polygon(psₛₗ)..., linewidth=3, legend=false) +scatter!(p, unzip(psₛₗ)..., markersize=3) +contour!(p, xs, ys, euclid_dist) + +# add some gradients along boundary +li(t, p1, p2) = p1 + t*(p2-p1) # t in [0,1] +for t in range(1/100, 1/2, length=3) + pt = li(t, psₛₗ[2], psₛₗ[3]) + arrow!(pt, ForwardDiff.gradient(euclid_dist, pt)) + pt = li(t, psₛₗ[2], psₛₗ[1]) + arrow!(pt, ForwardDiff.gradient(euclid_dist, pt)) +end + +p +@show 266 +#| hold : true +li(t, p1, p2) = p1 + t*(p2-p1) +p = plot(legend=false) +for i in 1:2, j in (i+1):3 + plot!(p, t -> euclid_dist(li(t, psₛₗ[i], psₛₗ[j]); ps=psₛₗ), 0, 1) +end +p +@show 280 +@syms xₗₛ[1:3] yₗₛ[1:3] α β +li(x, alpha, beta) = alpha + beta * x +d₂(alpha, beta) = sum((y - li(x, alpha, beta))^2 for (y,x) in zip(yₗₛ, xₗₛ)) +d₂(α, β) +@show 282 +grad_d₂ = diff.(d₂(α, β), [α, β]) +@show 283 +outₗₛ = solve(grad_d₂, [α, β]) +@show 285 +subs(outₗₛ[β], sum(xₗₛ) => 0) +@show 292 +[k => subs(v, xₗₛ[1]=>1, yₗₛ[1]=>1, xₗₛ[2]=>2, yₗₛ[2]=>3, + xₗₛ[3]=>5, yₗₛ[3]=>8) for (k,v) in outₗₛ] +@show 302 +f₂(x,y) = -exp(-((x-1)^2 + 2(y-1/2)^2)) +f₂(x) = f₂(x...) + +xs₂ = [[0.0, 0.0]] # we store a vector +gammas₂ = [1.0] + +for n in 1:5 + xn = xs₂[end] + gamma₀ = gammas₂[end] + xn1 = xn - gamma₀ * gradient(f₂)(xn) + dx, dy = xn1 - xn, gradient(f₂)(xn1) - gradient(f₂)(xn) + gamman1 = abs( (dx ⋅ dy) / (dy ⋅ dy) ) + + push!(xs₂, xn1) + push!(gammas₂, gamman1) +end + +[(x, f₂(x)) for x in xs₂] +@show 304 +#| hold: true +function surface_contour(xs, ys, f; offset=0) + p = surface(xs, ys, f, legend=false, fillalpha=0.5) + + ## we add to the graphic p, then plot + zs = [f(x,y) for x in xs, y in ys] # reverse order for use with Contour package + 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 = offset * _xs + plot!(p, _xs, _ys, _zs, alpha=0.5) # add curve on x-y plane + end + end + p +end + + +offset = 0 +us = vs = range(-1, 2, length=100) +surface_contour(us, vs, f₂, offset=offset) +pts = [[pt..., offset] for pt in xs₂] +scatter3d!(unzip(pts)...) +plot!(unzip(pts)..., linewidth=3) +@show 314 +function peaks(x, y) + z = 3 * (1 - x)^2 * exp(-x^2 - (y + 1)^2) + z += -10 * (x / 5 - x^3 - y^5) * exp(-x^2 - y^2) + z += -1/3 * exp(-(x+1)^2 - y^2) + return z +end +peaks(v) = peaks(v...) +@show 315 +#| hold: true +xs = range(-3, stop=3, length=100) +ys = range(-2, stop=2, length=100) +Ps = surface(xs, ys, peaks, legend=false) +Pc = contour(xs, ys, peaks, legend=false) +plot(Ps, Pc, layout=2) # combine plots +@show 319 +function newton_stepₚ(f, x) + M = ForwardDiff.hessian(f, x) + b = ForwardDiff.gradient(f, x) + x - M \ b +end +@show 321 +xₚ = [0, 1.5] +xₚ = newton_stepₚ(peaks, xₚ) +xₚ = newton_stepₚ(peaks, xₚ) +xₚ = newton_stepₚ(peaks, xₚ) +xₚ, ForwardDiff.gradient(peaks, xₚ) +@show 323 +Hₚ = ForwardDiff.hessian(peaks, xₚ) +@show 325 +#| hold: true +fxx = Hₚ[1,1] +d = det(Hₚ) +fxx, d +@show 335 +#| hold: true +g(x,y) = x^2 + 2y^2 -1 +g(v) = g(v...) + +xs = range(-3, 3, length=100) +ys = range(-1, 4, length=100) + +p = plot(aspect_ratio=:equal, legend=false) +contour!(xs, ys, g, levels=[0]) + +gi(x) = sqrt(1/2*(1-x^2)) # solve for y in terms of x +pts = [[x, gi(x)] for x in (-3/4, -1/4, 1/4, 3/4)] + +for pt in pts + arrow!(pt, ForwardDiff.gradient(g, pt) ) +end + +p +@show 338 +#| hold: true +#| echo: false +r(t) = [cos(t), sin(t)/2] +plot_parametric(pi/12..pi/3, r, legend=false, aspect_ratio=true, linewidth=3) +T(t) = -r'(t) / norm(r'(t)) +No(t) = T'(t) / norm(T'(t)) +t = pi/4 +lambda=1/10 +scatter!(unzip([r(t)])...) +arrow!(r(t), T(t)*lambda) +arrow!(r(t), No(t)* lambda) + +f(x,y)= x^2 + y^2 +f(v) = f(v...) +arrow!(r(t), lambda*ForwardDiff.gradient(f, r(t))) + +xs = range(0.5,1, length=100) +ys = range(0.1, 0.5, length=100) +contour!(xs, ys, f) +@show 344 +#| hold: true +#| echo: false +r(t) = [cos(t), sin(t)/2] +plot_parametric(-pi/6..pi/6,r, legend=false, aspect_ratio=true, linewidth=3) +T(t) = -r'(t) / norm(r'(t)) +No(t) = T'(t) / norm(T'(t)) +t = 0 +lambda=1/10 +scatter!(unzip([r(t)])...) +arrow!(r(t), T(t)*lambda) +arrow!(r(t), No(t)* lambda) + +f(x,y)= x^2 + y^2 +f(v) = f(v...) +arrow!(r(t), lambda*ForwardDiff.gradient(f, r(t))) + +xs = range(0.5,1.5, length=100) +ys = range(-0.5, 0.5, length=100) +contour!(xs, ys, f, levels = [.7, .85, 1, 1.15, 1.3]) +@show 381 +@syms lambda +fₗₐ(x, y) = x^2 - y^2 +gₗₐ(x, y) = x^2 + y^2 +Lₗₐ(x, y, lambda) = fₗₐ(x,y) - lambda * (gₗₐ(x,y) - 1) +dsₗₐ = solve(diff.(Lₗₐ(x, y, lambda), [x, y, lambda])) +@show 383 +[fₗₐ(d[x], d[y]) for d in dsₗₐ] +@show 432 +#| hold: true +@syms y y′ λ C +ex = Eq(-λ*y′^2/sqrt(1 + y′^2) + λ*sqrt(1 + y′^2), y - C) +Δ = sqrt(1 + y′^2) / (y - C) +ex1 = Eq(simplify(ex.lhs()*Δ), simplify(ex.rhs() * Δ)) +ex2 = Eq(ex1.lhs()^2 - 1, simplify(ex1.rhs()^2) - 1) +@show 457 +@syms z lambda1 lambda2 +g1(x, y, z) = x^2 + y^2 - z^2 +g2(x, y, z) = x - 2z - 3 +fₘ(x,y,z)= x^2 + y^2 + z^2 +Lₘ(x,y,z,lambda1, lambda2) = fₘ(x,y,z) - lambda1*(g1(x,y,z) - 0) - lambda2*(g2(x,y,z) - 0) + +∇Lₘ = diff.(Lₘ(x,y,z,lambda1, lambda2), [x, y, z,lambda1, lambda2]) +@show 459 +solve(subs.(∇Lₘ, lambda1 .=> 1)) +@show 461 +outₘ = solve(subs.(∇Lₘ, y .=> 0)) +@show 463 +[fₘ(d[x], 0, d[z]) for d in outₘ] +@show 498 +struct MultiIndex + alpha::Vector{Int} + end +Base.show(io::IO, α::MultiIndex) = println(io, "α = ($(join(α.alpha, ", ")))") + +## |α| = α_1 + ... + α_m +Base.length(α::MultiIndex) = sum(α.alpha) + +## factorial(α) computes α! +Base.factorial(α::MultiIndex) = prod(factorial(Sym(a)) for a in α.alpha) + +## x^α = x_1^α_1 * x_2^α^2 * ... * x_n^α_n +import Base: ^ +^(x, α::MultiIndex) = prod(u^a for (u,a) in zip(x, α.alpha)) + +## ∂^α(ex) = ∂_1^α_1 ∘ ∂_2^α_2 ∘ ... ∘ ∂_n^α_n (ex) +partial(ex::SymPy.SymbolicObject, α::MultiIndex, vars=free_symbols(ex)) = diff(ex, zip(vars, α.alpha)...) +@show 499 +@syms w +alpha = MultiIndex([1,2,1,3]) +length(alpha) # 1 + 2 + 1 + 3=7 +[1,2,3,4]^alpha +exₜ = x^3 * cos(w*y*z) +partial(exₜ, alpha, [w,x,y,z]) +@show 501 +struct MultiIndices + n::Int + k::Int +end + +function Base.length(as::MultiIndices) + n,k = as.n, as.k + n == 1 && return 1 + sum(length(MultiIndices(n-1, j)) for j in 0:k) # recursively identify length +end + +function Base.iterate(alphas::MultiIndices) + k, n = alphas.k, alphas.n + n == 1 && return ([k],(0, MultiIndices(0,0), nothing)) + + m = zeros(Int, n) + m[1] = k + betas = MultiIndices(n-1, 0) + stb = iterate(betas) + st = (k, MultiIndices(n-1, 0), stb) + return (m, st) +end + +function Base.iterate(alphas::MultiIndices, st) + + st == nothing && return nothing + k,n = alphas.k, alphas.n + k == 0 && return nothing + n == 1 && return nothing + + # can we iterate the next on + bk, bs, stb = st + + if stb==nothing + bk = bk-1 + bk < 0 && return nothing + bs = MultiIndices(bs.n, bs.k+1) + val, stb = iterate(bs) + return (vcat(bk,val), (bk, bs, stb)) + end + + resp = iterate(bs, stb) + if resp == nothing + bk = bk-1 + bk < 0 && return nothing + bs = MultiIndices(bs.n, bs.k+1) + val, stb = iterate(bs) + return (vcat(bk, val), (bk, bs, stb)) + end + + val, stb = resp + return (vcat(bk, val), (bk, bs, stb)) + +end +@show 503 +collect(MultiIndices(2, 3)) +@show 505 +union((collect(MultiIndices(2, i)) for i in 0:3)...) +@show 507 +k = 4 +length(MultiIndices(3, k+1)) +@show 509 +#| hold: true +@syms 𝐅() a[1:3] dx[1:3] + +sum(partial(𝐅(a...), α, a) / factorial(α) * dx^α for k in 0:3 for α in MultiIndex.(MultiIndices(3, k))) # 3rd order +@show 513 +#| hold: true +#| echo: false +f(x,y) = sqrt(x + y) +f(v) = f(v...) +pt = [2,2] +dxdy = [.1, .2] +val = f(pt) + dot(ForwardDiff.gradient(f, pt), dxdy) +numericq(val) +@show 516 +#| hold: true +#| echo: false +f(x,y,z) = x*y + y*z + z*x +f(v) = f(v...) +pt = [1,1,1] +dx = [0.1, 0.0, -0.1] +val = f(pt) + ∇(f)(pt) ⋅ dx +numericq(val) +@show 519 +#| hold: true +#| echo: false +f(x,y,z) = x*y + y*z + z*x - 8 +f(v) = f(v...) +pt = [1,1,1] +n = ∇(f)(pt) +d = dot(n, pt) +choices = [ + raw"`` x + y + z = 3``", + raw"`` 2x + y - 2z = 1``", + raw"`` x + 2y + 3z = 6``" +] +answ = 1 +radioq(choices, answ) +@show 523 +#| hold: true +#| echo: false +choices = [ + raw"`` \langle 2xy + y^2 + y, 2xy + x^2 + x\rangle``", + raw"`` y^2 + y, x^2 + x``", + raw"`` \langle 2y + y^2, 2x + x^2``" +] +answ = 1 +radioq(choices, answ) +@show 527 +#| hold: true +#| echo: false +yesnoq(true) +@show 529 +#| hold: true +#| echo: false +f(x,y) = x*y + x*y^2 + x^2 * y +f(v) = f(v...) +val = det(ForwardDiff.hessian(f, [-1/3, -1/3])) +numericq(val) +@show 531 +#| hold: true +#| echo: false +choices = [ + L"The function $f$ has a local minimum, as $f_{xx} > 0$ and $d >0$", + L"The function $f$ has a local maximum, as $f_{xx} < 0$ and $d >0$", + L"The function $f$ has a saddle point, as $d < 0$", + L"Nothing can be said, as $d=0$" +] +answ = 2 +radioq(choices, answ, keep_order=true) +@show 535 +#| hold: true +#| results: "hidden" +f(x,y) = x + 2x^2 + x^3 + y + 2x*y + y^2 +@syms x::real y::real +gradf = gradient(f(x,y), [x,y]) +@show 536 +#| hold: true +#| echo: false +yesnoq(true) +@show 538 +#| hold: true +#| results: "hidden" +f(x,y) = x + 2x^2 + x^3 + y + 2x*y + y^2 +@syms x::real y::real +gradf = gradient(f(x,y), [x,y]) + +solve(gradf, [x,y]) +@show 539 +#| hold: true +#| echo: false +numericq(2) +@show 541 +#| hold: true +f(x,y) = x + 2x^2 + x^3 + y + 2x*y + y^2 +@syms x::real y::real +gradf = gradient(f(x,y), [x,y]) + +sympy.hessian(f(x,y), [x,y]) +@show 543 +#| hold: true +#| echo: false +choices = [ + L"The function $f$ has a local minimum, as $f_{xx} > 0$ and $d >0$", + L"The function $f$ has a local maximum, as $f_{xx} < 0$ and $d >0$", + L"The function $f$ has a saddle point, as $d < 0$", + L"Nothing can be said, as $d=0$", + L"The test does not apply, as $\nabla{f}$ is not $0$ at this point." +] +answ = 3 +radioq(choices, answ, keep_order=true) +@show 545 +#| hold: true +#| echo: false +choices = [ + L"The function $f$ has a local minimum, as $f_{xx} > 0$ and $d >0$", + L"The function $f$ has a local maximum, as $f_{xx} < 0$ and $d >0$", + L"The function $f$ has a saddle point, as $d < 0$", + L"Nothing can be said, as $d=0$", + L"The test does not apply, as $\nabla{f}$ is not $0$ at this point." +] +answ = 1 +radioq(choices, answ, keep_order=true) +@show 547 +#| hold: true +#| echo: false +choices = [ + L"The function $f$ has a local minimum, as $f_{xx} > 0$ and $d >0$", + L"The function $f$ has a local maximum, as $f_{xx} < 0$ and $d >0$", + L"The function $f$ has a saddle point, as $d < 0$", + L"Nothing can be said, as $d=0$", + L"The test does not apply, as $\nabla{f}$ is not $0$ at this point." +] +answ = 5 +radioq(choices, answ, keep_order=true) +@show 553 +#| hold: true +#| echo: false +yesnoq(true) +@show 557 +#| hold: true +#| echo: false +yesnoq(false) +@show 559 +#| hold: true +#| echo: false +choices =[ + "It is the determinant of the Hessian", + L"It isn't, $b^2-4ac$ is from the quadratic formula" +] +answ = 1 +radioq(choices, answ) +@show 561 +#| hold: true +#| echo: false +choices = [ + L"That $a>0$ and $4ac-b^2 > 0$", + L"That $a<0$ and $4ac-b^2 > 0$", + L"That $4ac-b^2 < 0$" +] +answ = 2 +radioq(choices, answ, keep_order=true) +@show 563 +#| hold: true +#| echo: false +choices = [ + L"That $a>0$ and $4ac-b^2 > 0$", + L"That $a<0$ and $4ac-b^2 > 0$", + L"That $4ac-b^2 < 0$" +] +answ = 3 +radioq(choices, answ, keep_order=true) +@show 569 +#| hold: true +#| echo: false +yesnoq(true) +@show 571 +#| echo: false +choices = [ + raw"`` \langle 2x, 2y\rangle``", + raw"`` \langle 2x, y^2\rangle``", + raw"`` \langle x^2, 2y \rangle``" +] +answ = 1 +radioq(choices, answ) +@show 573 +f(x,y) = exp(-x^2-y^2) * (2x^2 + y^2) +f(v) = f(v...) +r(t) = sqrt(3)*[cos(t), sin(t)] +rat(x) = abs(x[1]/x[2]) - 1 +fn = rat ∘ ∇(f) ∘ r +ts = fzeros(fn, 0, 2pi) +@show 575 +#| eval: false +#| echo: false +f(x,y) = exp(-x^2-y^2) * (2x^2 + y^2) +r(t) = sqrt(3)*[cos(t), sin(t)] +rat(x) = abs(x[1]/x[2]) - 1 +fn = rat ∘ ∇(splat(f)) ∘ r +ts = fzeros(fn, 0, 2pi) + +val = maximum((splat(u)∘r).(ts)) +numericq(val) diff --git a/quarto/differentiable_vector_calculus/test.qmd b/quarto/differentiable_vector_calculus/test.qmd new file mode 100644 index 0000000..37ab0eb --- /dev/null +++ b/quarto/differentiable_vector_calculus/test.qmd @@ -0,0 +1,98 @@ +# Applications with scalar functions + + +{{< include ../_common_code.qmd >}} + +This section uses these add-on packages: + + +```{julia} +using CalculusWithJulia +using Plots +plotly() +using SymPy +using Roots +``` + +##### Example + + +Consider the function $f(x,y) = x^2 + 3y^2 -x$ over the region $x^2 + y^2 \leq 1$. This is a continuous function over a closed set, so will have both an absolute maximum and minimum. Find these from an investigation of the critical points and the boundary points. + + +The gradient is easily found: $\nabla{f} = \langle 2x - 1, 6y \rangle$, and is $\vec{0}$ only at $\vec{a} = \langle 1/2, 0 \rangle$. The Hessian is: + + +$$ +H = +\begin{bmatrix} +2 & 0\\ +0 & 6 +\end{bmatrix}. +$$ + +At $\vec{a}$ this has positive determinant and $f_{xx} > 0$, so $\vec{a}$ corresponds to a *local* minimum with values $f(\vec{a}) = (1/2)^2 + 3(0) - 1/2 = -1/4$. The absolute maximum and minimum may occur here (well, not the maximum) or on the boundary, so that must be considered. In this case we can easily parameterize the boundary and turn this into the univariate case: + + +```{julia} +fₗ(x,y) = x^2 + 2y^2 - x +gammaₗ(t) = [cos(t), sin(t)] # traces out x^2 + y^2 = 1 over [0, 2pi] +gₗ = splat(fₗ) ∘ gammaₗ + +cpsₗ = find_zeros(gₗ', 0, 2pi) # critical points of g +append!(cpsₗ, [0, 2pi]) +unique!(cpsₗ) +gₗ.(cpsₗ) +``` + +We see that maximum value is `2.25` and that the interior point, $\vec{a}$, will be where the minimum value occurs. To see exactly where the maximum occurs, we look at the values of gamma: + +```{julia} +inds = [2,4] +cpsₗ[inds] +``` + +These are multiples of $\pi$: + + +```{julia} +cpsₗ[inds]/pi +``` + +So we have the maximum occurs at the angles $2\pi/3$ and $4\pi/3$. Here we visualize, using a hacky trick of assigning `NaN` values to the function to avoid plotting outside the circle: + + +```{julia} +hₗ(x,y) = fₗ(x,y) * (x^2 + y^2 <= 1 ? 1 : NaN) +``` + +```{julia} +gr() +``` + +```{julia} +#| hold: true +xs = ys = range(-1,1, length=100) +plt = surface(xs, ys, hₗ) + +ts = cpsₗ # 2pi/3 and 4pi/3 by above +xs, ys = cos.(ts), sin.(ts) +zs = fₗ.(xs, ys) +scatter3d!(xs, ys, zs) +#tuple.(xs, ys, zs) +#plot!(plt, tuple.(xs, ys, zs); linetype=:scatter) +#plt +``` + +A contour plot also shows that some---and only one---extrema happens on the interior: + + diff --git a/quarto/differentiable_vector_calculus/vector_fields.qmd b/quarto/differentiable_vector_calculus/vector_fields.qmd index 70136e0..c315d8a 100644 --- a/quarto/differentiable_vector_calculus/vector_fields.qmd +++ b/quarto/differentiable_vector_calculus/vector_fields.qmd @@ -9,7 +9,7 @@ This section uses these add-on packages: ```{julia} using CalculusWithJulia using Plots -plotly() +#plotly() using SymPy using ForwardDiff using LinearAlgebra @@ -454,7 +454,7 @@ surface(unzip(surf.(ts, θs'))...; legend=false) ```{julia} #| echo: false -plotly() +#plotly() nothing ```