There are a few options for plotting equations in `Julia`. We will use `ImplicitPlots` in this section, but note both `ImplicitEquations` and `IntervalConstraintProgramming` offer alternatives that are a bit more flexible.
To plot an implicit equation using `ImplicitPlots` requires expressing the relationship in terms of a function, and then plotting the equation `f(x,y) = 0`. In practice this simply requires all the terms be moved to one side of an equals sign.
To plot the circle of radius ``2``, or the equations ``x^2 + y^2 = 2^2`` we would move all terms to one side ``x^2 + y^2 - 2^2 = 0`` and then express the left hand side through a function:
The `f` is a function of *two* variables, used here to express one side of an equation. `Julia` makes this easy to do - just make sure two variables are in the signature of `f` when it is defined. Using functions like this, we can express our equation in the form ``f(x,y) = c`` or, more generally, as ``f(x,y) = g(x,y)``. The latter of which can be expressed as ``h(x,y) = f(x,y) - g(x,y) = 0``. That is, only the form ``f(x,y)=0`` is needed to represent an equation.
There are two different styles in `Julia` to add simple plot recipes. `ImplicitPlots` adds a new plotting function (`implicit_plot`); alternatively many packages add a new recipe for the generic `plot` method using new types. (For example, `SymPy` has a plot recipe for symbolic types.
need us to express ``y`` as a function of ``x`` and then differentiate
that function. Of course, in this example, we would need two functions:
``f(x) = \sqrt{1-x^2}`` and ``g(x) = - \sqrt{1-x^2}`` to do this
completely.
In general though, we may not be able to solve for ``y`` in terms of ``x``. What then?
The idea is to *assume* that ``y`` is representable by some function of
``x``. This makes sense, moving on the curve from ``(x,y)`` to some nearby
point, means changing ``x`` will cause some change in ``y``. This
assumption is only made *locally* - basically meaning a complicated
graph is reduced to just a small, well-behaved, section of its graph.
With this assumption, asking what ``dy/dx`` is has an obvious meaning -
what is the slope of the tangent line to the graph at ``(x,y)``. (The assumption eliminates the question of what a tangent line would mean when a graph self intersects.)
The method of implicit differentiation allows this question to be
investigated. It begins by differentiating both sides of the equation
assuming ``y`` is a function of ``x`` to derive a new equation involving ``dy/dx``.
For example, starting with ``x^2 + y^2 = 1``, differentiating both sides in ``x`` gives:
The chain rule was used to find ``(d/dx)(y^2) = [y(x)^2]' = 2y \cdot dy/dx``. From this we can solve for ``dy/dx`` (the resulting equations are linear in ``dy/dx``, so can always be solved explicitly):
Basically this includes all the same steps as if done "by hand." Some effort could have been saved in plotting, had
values for the parameters been substituted initially, but not doing so
shows their dependence in the derivative.
```julia; echo=false
alert("The use of `lambdify(H)` is needed to turn the symbolic expression, `H`, into a function.")
```
```julia; echo=false
note("""
While `SymPy` itself has the `plot_implicit` function for plotting implicit equations, this works only with `PyPlot`, not `Plots`, so we use the `ImplicitPlots` package in these examples.
""")
```
## Higher order derivatives
Implicit differentiation can be used to find ``d^2y/dx^2`` or other higher-order derivatives. At each stage, the same technique is applied. The
only "trick" is that some simplifications can be made.
For example, consider ``x^3 - y^3=3``. To find ``d^2y/dx^2``, we first find ``dy/dx``:
```math
3x^2 - (3y^2 \frac{dy}{dx}) = 0.
```
We could solve for ``dy/dx`` at this point - it always appears as a linear factor - to get:
* as we change quadrants from the third to the fourth to the first the
concavity changes from down to up to down, as the sign of the second
derivative changes from negative to positive to negative;
* and that at these inflection points, the "tangent" line is vertical
when ``y=0`` and flat when ``x=0``.
```julia;
K(x,y) = x^3 - y^3 - 3
implicit_plot(K, xlims=(-3, 3), ylims=(-3, 3))
```
The same problem can be done symbolically. The steps are similar, though the last step (replacing ``x^3 - y^3`` with ``3``) isn't done without explicitly asking.
As [mentioned](../precalc/inversefunctions.html), an [inverse](http://en.wikipedia.org/wiki/Inverse_function) function for ``f(x)`` is a function ``g(x)`` satisfying:
``y = f(x)`` if and only if ``g(y) = x`` for all ``x`` in the domain of ``f`` and ``y`` in the range of ``f``.
In short, both ``f \circ g`` and ``g \circ f`` are identify functions on their respective domains.
As inverses are unique, their notation, ``f^{-1}(x)``, reflects the name of the related function.
The chain rule can be used to give the derivative of an inverse
function when applied to ``f(f^{-1}(x)) = x``. Solving gives,
This is great - if we can remember the rules. If not, sometimes implicit
differentiation can also help.
Consider the inverse function for the tangent, which exists when the domain of the tangent function is restricted to ``(-\pi/2, \pi/2)``. The function solves ``y = \tan^{-1}(x)`` or ``\tan(y) = x``. Differentiating this yields:
But ``\sec(y)^2 = 1 + \tan(y)^2 = 1 + x^2``, as can be seen by right-triangle trigonometry. This yields the formula ``dy/dx = [\tan^{-1}(x)]' = 1 / (1 + x^2)``.
##### Example
For a more complicated example, suppose we have a moving trajectory ``(x(t),
y(t))``. The angle it makes with the origin satisfies
```math
\tan(\theta(t)) = \frac{y(t)}{x(t)}.
```
Suppose ``\theta(t)`` can be defined in terms of the inverse to some
function (``\tan^{-1}(x)``). We can differentiate implicitly to find ``\theta'(t)`` in
But ``\sec^2(\theta(t)) = (r(t)/x(t))^2 = (x(t)^2 + y(t)^2) / x(t)^2``, so moving to the other side the secant term gives an explicit, albeit complicated, expression for the derivative of ``\theta`` in terms of the functions ``x`` and ``y``:
Okay, now to find the lowest point. This will be when the derivative
is ``0``. We solve by assuming ``y`` is a function of ``x`` called `u`. We have already defined symbolic variables `a`, `b`, `x`, and `y`, here we define `L`:
```julia;
@syms L
```
Then
```julia;
eqn = F₀(x,y,a,b) - L
```
```julia;
eqn_1 = diff(eqn(y => u(x)), x)
eqn_2 = solve(eqn_1, diff(u(x), x))[1]
dydx₂ = eqn_2(u(x) => y)
```
We are looking for when the tangent line has ``0`` slope, or when
`dydx` is ``0``:
```julia;
cps = solve(dydx₂, x)
```
There are two answers, as we could guess from the graph, but we want the one for the smallest value of ``y``, which is the second.
The values of `dydx` depend on any pair ``(x,y)``, but our solution must
also satisfy the equation. That is for our value of ``x``, we need to find
the corresponding ``y``. This should be possible by substituting:
```julia;
eqn1 = eqn(x => cps[2])
```
We would try to solve `eqn1` for `y` with `solve(eqn1, y)`, but
`SymPy` can't complete this problem. Instead, we will approach this
numerically using `find_zero` from the `Roots` package. We make the above a function of `y` alone
```julia;
eqn2 = eqn1(a=>3, b=>3, L=>10)
ystar = find_zero(eqn2, -3)
```
Okay, now we need to put this value back into our expression for the `x` value and also substitute in for the parameters:
```julia;
xstar = N(cps[2](y => ystar, a =>3, b => 3, L => 3))
```
Our minimum is at `(xstar, ystar)`, as this graphic shows:
Now, were we lucky and just happened to take ``a=3``, ``b = 3`` in such a way to
make this work? Well, no. But convince yourself by doing the above for
different values of ``b``.
----
In the above, we started with ``F(x,y) = L`` and solved symbolically for ``y=f(x)`` so that ``F(x,f(x)) = L``. Then we took a derivative of ``f(x)`` and set this equal to ``0`` to solve for the minimum ``y`` values.
Here we try the same problem numerically, using a zero-finding approach to identify ``f(x))``.
Starting with ``F(x,y) = \sqrt{x^2 + y^2} + \sqrt{(x-1)^2 + (b-2)^2}`` and ``L=3``, we have:
```julia
F₁(x,y) = F₀(x,y, 1, 2) - 3 # a,b,L = 1,2,3
implicit_plot(F₁)
```
Trying to find the lowest ``y`` value we have from the graph it is near ``x=0.1``. We can do better.
First, we could try so solve for the ``f`` using `find_zero`. Here is one way:
```julia
f₀(x) = find_zero(y -> F₁(x, y), 0)
```
We use ``0`` as an initial guess, as the ``y`` value is near ``0``. More on this later. We could then just sample many ``x`` values between ``-0.5`` and ``1.5`` and find the one corresponding to the smallest ``t`` value:
```julia
findmin([f₀(x) for x ∈ range(-0.5, 1.5, length=100)])
```
This shows the smallest value is around ``-0.414`` and occurs in the ``33``rd position of the sampled ``x`` values. Pretty good, but we can do better. We just need to differentiate ``f``, solve for ``f'(x) = 0`` and then put that value back into ``f`` to find the smallest ``y``.
**However** there is one subtle point. Using automatic differentiation, as implemented in `ForwardDiff`, with `find_zero` requires the `x0` initial value to have a certain type. In this case, the same type as the "`x`" passed into ``f(x)``. So rather than use an initial value of ``0``, we must use an initial value `zero(x)`! (Otherwise, there will be an error "`no method matching Float64(::ForwardDiff.Dual{...`".)
With this slight modification, we have:
```julia
f₁(x) = find_zero(y -> F₁(x, y), zero(x))
plot(f₁', -0.5, 1.5)
```
The zero of `f'` is a bit to the right of ``0``, say ``0.2``; we use `find_zero` again to find it:
```julia
xstar₁ = find_zero(f₁', 0.2)
xstar₁, f₁(xstar₁)
```
It is important to note that the above uses of `find_zero` required *good* initial guesses, which we were fortunate enough to identify.
## Questions
###### Question
Is ``(1,1)`` on the graph of
```math
x^2 - 2xy + y^2 = 1?
```
```julia; hold=true; echo=false
x,y=1,1
yesnoq(x^2 - 2x*y + y^2 ==1)
```
###### Question
For the equation
```math
x^2y + 2y - 4 x = 0,
```
if ``x=4``, what is a value for ``y`` such that ``(x,y)`` is a point on the graph of the equation?
An ellipse is when ``n=1``. Take ``n=3``, ``a=1``, and ``b=2``.
Find a *positive* value of ``y`` when ``x=1/2``.
```julia; hold=true; echo=false
a,b,n=1,2,3
val = b*(1 - ((1/2)/a)^n)^(1/n)
numericq(val)
```
What expression gives ``dy/dx``?
```julia; hold=true; echo=false
choices = [
"`` -(y/x) \\cdot (x/a)^n \\cdot (y/b)^{-n}``",
"``b \\cdot (1 - (x/a)^n)^{1/n}``",
"``-(x/a)^n / (y/b)^n``"
]
ans = 1
radioq(choices, ans)
```
###### Question
Let ``y - x^2 = -\log(x)``. At the point ``(1/2, 0.9431...)``, the graph has a tangent line. Find this line, then find its intersection point with the ``y`` axes.
This intersection is:
```julia; hold=true; echo=false
f(x) = x^2 - log(x)
x0 = 1/2
tl(x) = f(x0) + f'(x0) * (x - x0)
numericq(tl(0))
```
###### Question
The [witch](http://www-history.mcs.st-and.ac.uk/Curves/Witch.html) of [Agnesi](http://www.maa.org/publications/periodicals/convergence/mathematical-treasures-maria-agnesis-analytical-institutions) is the curve given by the equation
```math
y(x^2 + a^2) = a^3.
```
If ``a=1``, numerically find a a value of ``y`` when ``x=2``.
```julia; hold=true; echo=false
a = 1
f(x,y) = y*(x^2 + a^2) - a^3
val = find_zero(y->f(2,y), 1)
numericq(val)
```
What expression yields ``dy/dx`` for this curve:
```julia; hold=true; echo=false
choices = [
"``-2xy/(x^2 + a^2)``",
"``2xy / (x^2 + a^2)``",
"``a^3/(x^2 + a^2)``"
]
ans = 1
radioq(choices, ans)
```
###### Question
```julia; hold=true; echo=false
### {{{lhopital_35}}}
imgfile = "figures/fcarc-may2016-fig35-350.gif"
caption = """
Image number 35 from L'Hospitals calculus book (the first). Given a description of the curve, identify the point ``E`` which maximizes the height.
"""
ImageFile(:derivatives, imgfile, caption)
```
The figure above shows a problem appearing in L'Hospital's first calculus book. Given a function defined implicitly by ``x^3 + y^3 = axy`` (with ``AP=x``, ``AM=y`` and ``AB=a``) find the point ``E`` that maximizes the height. In the [AMS feature column](http://www.ams.org/samplings/feature-column/fc-2016-05) this problem is illustrated and solved in the historical manner, with the comment that the concept of implicit differentiation wouldn't have occurred to L'Hospital.
Using Implicit differentiation, find when ``dy/dx = 0``.
choices = ["concave up", "concave down", "both concave up and down"]
ans = 3
radioq(choices, ans, keep_order=true)
```
## Appendix
There are other packages in the `Julia` ecosystem that can plot implicit equations.
### The ImplicitEquations package
The `ImplicitEquations` packages can plot equations and inequalities. The use is somewhat similar to the examples above, but the object plotted is a predicate, not a function. These predicates are created with functions like `Eq` or `Lt`.
For example, the `ImplicitPlots` manual shows this function ``f(x,y) = (x^4 + y^4 - 1) \cdot (x^2 + y^2 - 2) + x^5 \cdot y`` to plot. Using `ImplicitEquations`, this equation would be plotted with:
The rendered plots look "blocky" due to the algorithm used to plot the
equations. As there is no rule defining ``(x,y)`` pairs to plot, a
search by regions is done. A region is initially labeled
undetermined. If it can be shown that for any value in the region the
equation is true (equations can also be inequalities), the region is
colored black. If it can be shown it will never be true, the region is
dropped. If a black-and-white answer is not clear, the region is
subdivided and each subregion is similarly tested. This continues
until the remaining undecided regions are smaller than some
threshold. Such regions comprise a boundary, and here are also colored
black. Only regions are plotted - not ``(x,y)`` pairs - so the results
are blocky. Pass larger values of ``N=M`` (with defaults of ``8``) to
`plot` to lower the threshold at the cost of longer computation times,
as seen in the last example.
### The IntervalConstraintProgramming package
The `IntervalConstraintProgramming` package also can be used to graph implicit equations. For certain problem descriptions it is significantly faster and makes better graphs. The usage is slightly more involved. We show the commands, but don't run them here, as there are minor conflicts with the `CalculusWithJulia`package.
We specify a problem using the `@constraint` macro. Using a macro
allows expressions to involve free symbols, so the problem is
specified in an equation-like manner:
```julia; eval=false
S = @constraint x^2 + y^2 <= 2
```
The right hand side must be a number.
The area to plot over must be specified as an `IntervalBox`, basically a pair of intervals. The interval ``[a,b]`` is expressed through `a..b`:
```julia; eval=false
J = -3..3
X = IntervalArithmetic.IntervalBox(J, J)
```
The `pave` command does the heavy lifting:
```julia; eval=false
region = IntervalConstraintProgramming.pave(S, X)
```
A plot can be made of either the boundary, the interior, or both.
```julia; eval=false
plot(region.inner) # plot interior; use r.boundary for boundary
