algorithm can be used for polynomials $a(x)$ and $b(x)$ to produce $q(x)$ and
$r(x)$ with $a = b\cdot q + r$ *and* the degree of $r(x)$ is less than the
degree of $b(x)$. This is in direct analogy to the division algorithm of
integers, only there the value of the remainder, $r(x)$, satisfies $0
\leq r < b$. Given $q(x)$ and $r(x)$ as above, we can reexpress the rational function
```math
\frac{a(x)}{b(x)} = q(x) + \frac{r(x)}{b(x)}.
```
The rational expression on the right-hand side has larger degree in the denominator.
The division algorithm is implemented in `Julia` generically through the `divrem` method:
```julia;
q, r = divrem(num, den)
```
This yields the decomposition of `num/den`:
```julia;
q + r/den
```
A similar result can be found using the `apart` function, which can be easier to use if the expression is not given in terms of a separate numerator and denominator.
```julia;
g(x) = (x-1)^2 * (x-2) / ((x+3)*(x-3)) # as a function
h = g(x) # a symbolic expression
apart(h)
```
This decomposition breaks the rational expression into two pieces: $x-4$
and $40/(3x+9) + 2/(3x-9)$. The first piece would have a graph
that is the line with slope $1$ and $y$-intercept $4$. As $x$ goes to
$\infty$, the second piece will clearly go towards ``0,`` as this simple
graph shows:
```julia;
plot(apart(h) - (x - 4), 10, 100)
```
Similarly, a plot over $[-100, -10]$ would show decay towards $0$,
though in that case from below. Combining these two facts then, it is
now no surprise that the graph of the rational function $f(x)$ should
approach a straight line, in this case $y=x-4$ as $x \rightarrow \pm
\infty$.
We can easily do most of this analysis without needing a
computer or algebra. First, we should know the four eventual shapes of a polynomial, that the graph of $y=mx$ is a line with
slope $m$, the graph of $y = c$ is a constant line at height $c$, and
the graph of $y=c/x^m$, $m > 0$ will decay towards $0$ as $x
\rightarrow \pm\infty$. The latter should be clear, as $x^m$ gets big,
so its reciprocal goes towards $0$.
The factored form, as $p$ is presented, is a bit hard to work with,
rather we use the expanded form, which we get through the `cancel`
function
```julia;
cancel(h)
```
We can see that the numerator is of degree ``3`` and the denominator
of degree $2$. The leading terms are $x^3$ and $x^2$, respectively. If
The leading term of the numerator is $17x^5$ and the leading term of the denominator is $x^5$. The ratio is $17$ (or $17x^0 = 17x^{5-5}$). As such, we would have a horizontal asymptote $y=17$.
Then we can see that the ratio of the leading terms is $x^5 / (5x^4) = (1/5)x$. We expect a slant asymptote with slope $1/5$, though we would need to divide to see the exact intercept. This is found with, say:
where both $r(c)$ and $s(c)$ are non zero. Knowing $m$ and $n$ (the multiplicities of the root $c$) allows the following to be said:
* If $m < n$ then $x=c$ will be a vertical asymptote.
* If $m \geq n$ then $x=c$ will not be vertical asymptote. (The value
$c$ will be known as a removable singularity). In this case, the
graph of $p(x)/q(x)$ and the graph of $(x-c)^{m-n}r(x)/s(x)$ will
differ, though very slightly, as the latter will include a value for
$x=c$, whereas $x=c$ is not in the domain of $p(x)/q(x)$.
Finding the multiplicity may or may not be hard, but there is a very
kludgy quick check that is often correct. With `Julia`, if you have a
rational function that has `f(c)` evaluate to `Inf` or `-Inf` then
there will be a vertical asymptote. If the expression evaluates to
`NaN`, more analysis is needed. (The value of `0/0` is `NaN`, where as
`1/0` is `Inf`.)
For example, the function
``f(x) = ((x-1)^2 \cdot (x-2)) / ((x+3) \cdot(x-3))`` has vertical asymptotes at ``-3`` and ``3``, as its graph
illustrated. Without the graph we could see this as well:
```julia; hold=true;
f(x) = (x-1)^2 * (x-2) / ((x+3)*(x-3) )
f(3), f(-3)
```
#### Graphing with vertical asymptotes
As seen in several graphs, the basic plotting algorithm does a poor
job with vertical asymptotes. For example, it may erroneously connect
their values with a steep vertical line, or the $y$-axis scale can get
so large as to make reading the rest of the graph impossible. There are some tricks to work around this.
Consider again the function $f(x) = ((x-1)^2 \cdot (x-2)) / ((x+3)
\cdot(x-3))$. Without much work, we can see that $x=3$ and $x=-3$ will
be vertical asymptotes and there will be a slant asymptote with
slope ``1``. How to graph this?
We can avoid the vertical asymptotes in our viewing window. For
example we could look at the area between the vertical asymptotes, by
plotting over $(-2.9, 2.9)$, say:
```julia;
𝒇(x) = (x-1)^2 * (x-2) / ((x+3)*(x-3) )
plot(𝒇, -2.9, 2.9)
```
This backs off by $\delta = 0.1$. As we have that $3 - 2.9$ is
$\delta$ and $1/\delta$ is 10, the $y$ axis won't get too large, and
indeed it doesn't.
This graph doesn't show well the two zeros at $x=1$ and $x=2$, for
that a narrower viewing window is needed. By successively panning
throughout the interesting part of the graph, we can get a view of the
function.
We can also clip the `y` axis. The `plot` function can be passed an
argument `ylims=(lo, hi)` to limit which values are plotted. With this,
we can have:
```julia; hold=true
plot(𝒇, -5, 5, ylims=(-20, 20))
```
This isn't ideal, as the large values are still computed, just the
viewing window is clipped. This leaves the vertical asymptotes still
effecting the graph.
There is another way, we could ask `Julia` to not plot $y$ values that
get too large. This is not a big request. If instead of the value of
`f(x)` - when it is large - -we use `NaN` instead, then the
connect-the-dots algorithm will skip those values.
This was discussed in an earlier section where the `rangeclamp` function was introduced to replace large values of `f(x)` (in absolute values) with `NaN`.
```julia;
plot(rangeclamp(𝒇, 30), -25, 25) # rangeclamp is in the CalculusWithJulia package
```
We can see the general shape of ``3`` curves broken up by the vertical
asymptotes. The two on the side heading off towards the line $x-4$ and
the one in the middle. We still can't see the precise location of the
zeros, but that wouldn't be the case with most graphs that show
asymptotic behaviors. However, we can clearly tell where to "zoom in"
were those of interest.
### Sign charts
When sketching graphs of rational functions by hand, it is useful to use sign charts.
A sign chart of a function indicates when the function is positive,
negative, $0$, or undefined. It typically is represented along the
lines of this one for $f(x) = x^3 - x$:
```verbatim
- 0 + 0 - 0 +
< ----- -1 ----- 0 ----- 1 ----- >
```
The usual recipe for construction follows these steps:
- Identify when the function is $0$ or undefined. Place those values
on a number line.
- Identify "test points" within each implied interval (these are $(-\infty, -1)$, $(-1,0)$, $(0,1)$, and $(1, \infty)$ in the example) and check for the sign of $f(x)$ at these test points. Write in `-`, `+`, `0`, or `*`, as appropriate. The value comes from the fact that "continuous" functions may only change sign when they cross $0$ or are undefined.
With the computer, where it is convenient to draw a graph, it might be better to emphasize
the sign on the graph of the function. The `sign_chart` function from `CalculusWithJulia` does this by numerically identifying points where the function is ``0`` or ``\infty`` and indicating the sign as ``x`` crosses over these points.
## The `Polynomials` package for rational functions
In the following, we import some functions from the `Polynomials` package. We avoided loading the entire namespace, as there are conflicts with `SymPy`. Here we import some useful functions and the `Polynomial` constructor:
The `Polynomials` package has support for rational functions. The `//` operator can be used to create rational expressions:
```julia;
𝒙 = variable()
𝒑 = (𝒙-1)*(𝒙-2)^2
𝒒 = (𝒙-2)*(𝒙-3)
𝒑𝒒 = 𝒑 // 𝒒
```
A rational expression is a formal object; a rational function the viewpoint that this object will be evaluated by substituting values for the indeterminate. Rational expressions made within `Polynomials` are evaluated just like functions:
```julia;
𝒑𝒒(4) # p(4)/q(4)
```
The rational expressions are not in lowest terms unless requested through the `lowest_terms` method:
```julia;
lowest_terms(𝒑𝒒)
```
For polynomials as simple as these, this computation is not a problem,
but there is the very real possibility that the lowest term
computation may be incorrect. Unlike `SymPy` which factors
symbolically, `lowest_terms` uses a numeric algorithm and does not, as
would be done by hand or with `SymPy`, factor the polynomial and then cancel common
factors.
The distinction between the two expressions is sometimes made; the
initial expression is not defined at ``x=2``; the reduced one is, so
the two are not identical when viewed as functions of the variable
``x``.
Rational expressions include polynomial expressions, just as the
rational numbers include the integers. The identification there is to
divide by ``1``, thinking of ``3`` as ``3/1``. In `Julia`, we would
just use
```julia;
3//1
```
The integer can be recovered from the rational number using `numerator`:
```julia;
numerator(3//1)
```
Similarly, we can divide a polynomial by the polynomial ``1``, which in `Julia` is returned by `one(p)`, to produce a rational expression:
```julia;
pp = 𝒑 // one(𝒑)
```
And as with rational numbers, `p` is recovered by `numerator`:
```julia;
numerator(pp)
```
One difference is the rational number `3//1` also represents other
expressions, say `6/2` or `12/4`, as `Julia`'s rational numbers are
presented in lowest terms, unlike the rational expressions in
`Polynomials`.
Rational functions also have a plot recipe defined for them that
attempts to ensure the basic features are identifiable. As previously
discussed, a plot of a rational function can require some effort to
avoid the values associated to vertical asymptotes taking up too many
of the available vertical pixels in a graph.
For the polynomial `pq` above, we have from observation that ``1`` and
``2`` will be zeros and ``x=3`` a vertical asymptote. We also can
identify a slant asymptote with slope ``1``. These are hinted at in this graph:
```julia;
plot(𝒑𝒒)
```
To better see the zeros, a plot over a narrower interval, say ``[0,2.5]``,
would be encouraged; to better see the slant asymptote, a plot over
a wider interval, say ``[-10,10]``, would be encouraged.
For one more example of the default plot recipe, we redo the graphing of the rational expression we earlier plotted with `rangeclamp`:
```julia; hold=true;
p,q = fromroots([1,1,2]), fromroots([-3,3])
plot(p//q)
```
##### Example: transformations of polynomials; real roots
We have seen some basic transformations of functions such as shifts and scales. For a polynomial expression we can implement these as follows, taking advantage of polynomial evaluation:
```julia; hold=true;
x = variable()
p = 3 + 4x + 5x^2
a = 2
p(a*x), p(x+a) # scale, shift
```
A different polynomial transformation is inversion, or the mapping ``x^d \cdot p(1/x)`` where ``d`` is the degree of ``p``. This will yield a polynomial, as perhaps this example will convince you:
```julia; hold=true
p = Polynomial([1, 2, 3, 4, 5])
d = Polynomials.degree(p) # degree is in SymPy and Polynomials, indicate which
pp = p // one(p)
x = variable(pp)
q = x^d * pp(1/x)
lowest_terms(q)
```
We had to use a rational expression so that division by the variable was possible.
The above indicates that the new polynomial, ``q``, is constructed from ``p`` by **reversing** the coefficients.
Inversion is like a funhouse mirror, flipping around parts of the polynomial. For example, the interval
``[1/4,1/2]`` is related to the interval ``[2,4]``. Of interest here, is that if ``p(x)`` had a root, ``r``, in ``[1/4,1/2]`` then ``q(x) = x^d \cdot p(1/x)`` would have a root in ``[2,4]`` at ``1/r``.
So these three transformations -- scale, shift, and inversion -- can be defined for polynomials.
Combined, the three can be used to create a [Mobius transformation](https://en.wikipedia.org/wiki/M%C3%B6bius_transformation). For two values ``a`` and ``b``, consider the polynomial derived from ``p`` (again `d=degree(p)`) by:
```math
q = (x+1)^d \cdot p(\frac{ax + b}{x + 1}).
```
Here is a non-performant implementation as a `Julia` function:
```julia;
function mobius_transformation(p, a, b)
x = variable(p)
p = p(x + a) # shift
p = p((b-a)*x) # scale
p = Polynomial(reverse(coeffs(p))) # invert
p = p(x + 1) # shift
p
end
```
We can verify this does what we want through example with the previously defined `p`:
Mobius transforms are used to map regions into other regions. In this special case, the transform ``\phi(x) = (ax + b)/(x + 1)`` takes the interval ``[0,\infty]`` and sends it to ``[a,b]`` (``0`` goes to ``(a\cdot 0 + b)/(0+1) = b``, whereas ``\infty`` goes to ``ax/x \rightarrow a``). Using this, if ``p(u) = 0``, with ``q(x) = (x-1)^d p(\phi(x))``, then setting ``u = \phi(x)`` we have ``q(x) = (\phi^{-1}(u)+1)^d p(\phi(\phi^{-1}(u))) = (\phi^{-1}(u)+1)^d \cdot p(u) = (\phi^{-1}(u)+1)^d \cdot 0 = 0``. That is, a zero of ``p`` in ``[a,b]`` will appear as a zero of ``q`` in ``[0,\infty)`` at ``\phi^{-1}(u)``.
The Descartes rule of signs applied to ``q`` then will give a bound on the number of possible roots of ``p`` in the interval ``[a,b]``. In the example we did, the Mobius transform for ``a=4, b=6`` is ``15 - x - 11x^2 - 3x^3`` with ``1`` sign change, so there must be exactly ``1`` real root of ``p=(x-1)(x-3)(x-5)`` in the interval ``[4,6]``, as we can observe from the factored form of ``p``.
Similarly, we can see there are ``2`` or ``0`` roots for ``p`` in the interval ``[2,6]`` by counting the two sign changes here:
```julia;
mobius_transformation(𝐩, 2,6)
```
This observation, along with a detailed analysis provided by [Kobel, Rouillier, and Sagraloff](https://dl.acm.org/doi/10.1145/2930889.2930937) provides a means to find intervals that enclose the real roots of a polynomial.
The basic algorithm, as presented next, is fairly simple to understand, and hints at the bisection algorithm to come. It is due to Akritas and Collins. Suppose you know the only possible positive real roots are between ``0`` and ``M`` *and* no roots are repeated. Find the transformed polynomial over ``[0,M]``:
* If there are no sign changes, then there are no roots of ``p`` in ``[0,M]``.
* If there is one sign change, then there is a single root of ``p`` in ``[0,M]``. The interval ``[0,M]`` is said to isolate the root (and the actual root can then be found by other means)
* If there is more than one sign change, divide the interval in two (``[0,M/2]`` and ``[M/2,M]``, say) and apply the same consideration to each.
Eventually, **mathematically** this will find isolating intervals for each positive real root. (The negative ones can be similarly isolated.)
Applying these steps to ``p`` with an initial interval, say ``[0,9]``, we would have:
So the three roots (``1``, ``3``, ``5``) are isolated by ``[0, 9/4]``, ``[9/4, 9/2]``, and ``[9/2, 9]``.
### The `RealPolynomialRoots` package.
For square-free polynomials, the `RealPolynomialRoots` package implements a basic version of the paper of [Kobel, Rouillier, and Sagraloff](https://dl.acm.org/doi/10.1145/2930889.2930937) to identify the real roots of a polynomial using the Descartes rule of signs and the Möbius transformations just described.
The `ANewDsc` function takes a collection of coefficients representing a polynomial and returns isolating intervals for each real root. For example:
```julia
p₀ = fromroots([1,3,5])
st = ANewDsc(coeffs(p₀))
```
These intervals can be refined to give accurate approximations to the roots:
```julia
refine_roots(st)
```
More challenging problems can be readily handled by this package. The following polynomial
```julia
𝒔 = Polynomial([0,1]) # also just variable(Polynomial{Int})
𝒖 = -1 + 254*𝒔 - 16129*𝒔^2 + 𝒔^15
```
has three real roots, two of which are clustered very close to each other:
```julia
𝒔𝒕 = ANewDsc(coeffs(𝒖))
```
and
```julia
refine_roots(𝒔𝒕)
```
The SymPy package (`sympy.real_roots`) can accurately identify the
three roots but it can take a **very** long time. The
`Polynomials.roots` function from the `Polynomials` package identifies
the cluster as complex valued. Though the implementation in
`RealPolynomialRoots` doesn't handle such large polynomials, the
authors of the algorithm have implementations that can quickly solve
polynomials with degrees as high as ``10,000``.
## Questions
###### Question
The rational expression $(x^3 - 2x + 3) / (x^2 - x + 1)$ would have
(The wikipedia page indicates that the term "asymptote" was introduced
by Apollonius of Perga in his work on conic sections, but in contrast
to its modern meaning, he used it to mean any line that does not
intersect the given curve. It can sometimes take a while to change perception.)
###### Question
Consider the two graphs of $f(x) = 1/x$ over $[10,20]$ and $[100, 200]$:
```julia; hold=true; echo=false
plot(x -> 1/x, 10, 20)
```
```julia; hold=true; echo=false
plot(x -> 1/x, 100, 200)
```
The two shapes are basically identical and do not look like straight lines. How does this reconcile with the fact that $f(x)=1/x$ has a horizontal asymptote $y=0$?
```julia; hold=true; echo=false
choices = ["The horizontal asymptote is not a straight line.",
L"The $y$-axis scale shows that indeed the $y$ values are getting close to $0$.",
L"The graph is always decreasing, hence it will eventually reach $-\infty$."
The (low-order) Pade approximation for $\sin(x)$ was seen to be $(x - 7/60 \cdot x^3)/(1 + 1/20 \cdot x^2)$. The graph showed that this approximation was fairly close
over $[-\pi, \pi]$. Without graphing would you expect the behaviour of the function and its approximation to be similar for *large* values of $x$?
```julia; hold=true; echo=false
yesnoq(false)
```
Why?
```julia; hold=true; echo=false
choices = [
L"The $\sin(x)$ oscillates, but the rational function eventually follows $7/60 \cdot x^3$",
L"The $\sin(x)$ oscillates, but the rational function has a slant asymptote",
L"The $\sin(x)$ oscillates, but the rational function has a non-zero horizontal asymptote",
L"The $\sin(x)$ oscillates, but the rational function has a horizontal asymptote of $0$"]
