- It imports some of the functionality provided by `SymPy`, including the ability to create symbolic variables.
- It overloads many `Julia` functions to work seamlessly with symbolic expressions. This makes working with polynomials quite natural.
- It gives access to a wide range of SymPy's functionality through the `sympy` object.
To illustrate, using the just defined `x`, here is how we can create the polynomial $-16x^2 + 100$:
```julia;
𝒑 = -16x^2 + 100
```
That is, the expression is created just as you would create it within
a function body. But here the result is still a symbolic object. We
have assigned this expression to a variable `p`, and have not defined
it as a function `p(x)`. Mentally keeping the distinction between symbolic
expressions and functions is very important.
The `typeof` function shows that `𝒑` is of a symbolic type (`Sym`):
```julia;
typeof(𝒑)
```
We can mix and match symbolic objects. This command creates an
arbitrary quadratic polynomial:
```julia;
quad = a*x^2 + b*x + c
```
Again, this is entered in a manner nearly identical to how we see such
expressions typeset ($ax^2 + bx+c$), though we must remember to
explicitly place the multiplication operator, as the symbols are not
numeric literals.
We can apply many of `Julia`'s mathematical functions and the result will still be symbolic:
```julia;
sin(a*(x - b*pi) + c)
```
Another example, might be the following combination:
```julia;
quad + quad^2 - quad^3
```
One way to create symbolic expressions is simply to call a `Julia` function with symbolic arguments. The first line in the next example defines a function, the second evaluates it at the symbols `x`, `a`, and `b` resulting in a symbolic expression `ex`:
```julia
f(x, m, b) = m*x + b
ex = f(x, a, b)
```
## Substitution: subs, replace
Algebraically working with symbolic expressions is straightforward. A
different symbolic task is substitution. For example, replacing each
instance of `x` in a polynomial, with, say, `(x-1)^2`. Substitution
requires three things to be specified: an expression to work on, a
variable to substitute, and a value to substitute in.
SymPy provides its `subs` function for this. This function is available in `Julia`, but it is easier to use notation reminiscent of function evaluation.
To illustrate, to do
the task above for the polynomial $-16x^2 + 100$ we could have:
```julia;
𝒑(x => (x-1)^2)
```
This "call" notation takes pairs (designated by `a=>b`) where the left-hand side is the variable to substitute for, and the right-hand side the new value.
The value to substitute can depend on the variable, as illustrated; be
a different variable; or be a numeric value, such as $2$:
```julia;
𝒚 = 𝒑(x=>2)
```
The result will always be of a symbolic type, even if the answer is
just a number:
```julia;
typeof(𝒚)
```
If there is just one free variable in an expression, the pair notation can be dropped:
```julia;
𝒑(4) # substitutes x=>4
```
##### Example
Suppose we have the polynomial $p = ax^2 + bx +c$. What would it look
like if we shifted right by $E$ units and up by $F$ units?
```julia;
@syms E F
p₂ = a*x^2 + b*x + c
p₂(x => x-E) + F
```
And expanded this becomes:
```julia;
expand(p₂(x => x-E) + F)
```
### Conversion of symbolic numbers to Julia numbers
In the above, we substituted `2` in for `x` to get `y`:
```julia; hold=true
p = -16x^2 + 100
y = p(2)
```
The value, $36$ is still symbolic, but clearly an integer. If we
are just looking at the output, we can easily translate from the
symbolic value to an integer, as they print similarly. However the
conversion to an integer, or another type of number, does not happen
automatically. If a number is needed to pass along to another `Julia`
function, it may need to be converted. In general, conversions between
different types are handled through various methods of
`convert`. However, with `SymPy`, the `N` function will attempt to do
the conversion for you:
```julia;hold=true
p = -16x^2 + 100
N(p(2))
```
Where `convert(T,x)` requires a specification of the type to convert `x` to, `N` attempts to match the data type used by SymPy to store the number. As such, the output type of `N` may vary (rational, a BigFloat, a float, etc.)
For getting more digits of accuracy, a
precision can be passed to `N`. The following command will take
the symbolic value for $\pi$, `PI`, and produce about ``60`` digits worth
as a `BigFloat` value:
```julia;
N(PI, 60)
```
Conversion by `N` will fail if the value to be converted contains free
symbols, as would be expected.
### Converting symbolic expressions into Julia functions
Evaluating a symbolic expression and returning a numeric value can be done by composing the two just discussed concepts. For example:
```julia;
𝐩 = 200 - 16x^2
N(𝐩(2))
```
This approach is direct, but can be slow *if* many such evaluations were needed (such as with a plot). An alternative is to turn the symbolic expression into a `Julia` function and then evaluate that as usual.
The `lambdify` function turns a symbolic expression into a `Julia` function
```julia;hold=true
pp = lambdify(𝐩)
pp(2)
```
The `lambdify` function uses the name of the similar `SymPy` function which is named after Pythons convention of calling anoynmous function "lambdas." The use above is straightforward. Only slightly more complicated is the use when there are multiple symbolic values. For example:
```julia; hold=true
p = a*x^2 + b
pp = lambdify(p)
pp(1,2,3)
```
This evaluation matches `a` with `1`, `b` with`2`, and `x` with `3` as that is the order returned by the function call `free_symbols(p)`. To adjust that, a second `vars` argument can be given:
```julia; hold=true
pp = lambdify(p, (x,a,b))
pp(1,2,3) # computes 2*1^2 + 3
```
## Graphical properties of polynomials
Consider the graph of the polynomial `x^5 - x + 1`:
```julia;
plot(x^5 - x + 1, -3/2, 3/2)
```
(Plotting symbolic expressions is similar to plotting a function, in
that the expression is passed in as the first argument. The expression
must have only one free variable, as above, or an error will occur.)
This graph illustrates the key features of polynomial graphs:
* there may be values for `x` where the graph crosses the $x$ axis
(real roots of the polynomial);
* there may be peaks and valleys (local maxima and local minima)
* except for constant polynomials, the ultimate behaviour for large
values of $\lvert x\rvert$ is either both sides of the graph going to positive
infinity, or negative infinity, or as in this graph one to the
positive infinity and one to negative infinity. In particular, there
is no *horizontal asymptote*.
To investigate this last point, let's consider the case of the
monomial $x^n$. When $n$ is even, the following animation shows that
larger values of $n$ have greater growth once outside of $[-1,1]$:
Suppose $p = a_n x^n + \cdots + a_1 x + a_0$ with $a_n > 0$. Then by
the above, eventually for large $x > 0$ we have $p > 0$, as that is the
behaviour of $a_n x^n$. Were $a_n < 0$, then eventually for large
$x>0$, $p < 0$.
Now consider the related polynomial, $q$, where we multiply $p$ by $x^n$ and substitute in $1/x$ for $x$. This is the "reversed" polynomial, as we see in this illustration for $n=2$:
```julia;hold=true
p = a*x^2 + b*x + c
n = 2 # the degree of p
q = expand(x^n * p(x => 1/x))
```
In particular, from the reversal, the behavior of $q$ for large $x$
depends on the sign of $a_0$. As well, due to the $1/x$, the behaviour
of $q$ for large $x>0$ is the same as the behaviour of $p$ for small
*positive* $x$. In particular if $a_n > 0$ but $a_0 < 0$, then `p` is
eventually positive and `q` is eventually negative.
That is, if $p$ has $a_n > 0$ but $a_0 < 0$ then the graph of $p$ must cross the $x$ axis.
This observation is the start of Descartes' rule of
[signs](http://sepwww.stanford.edu/oldsep/stew/descartes.pdf), which
counts the change of signs of the coefficients in `p` to say something
about how many possible crossings there are of the $x$ axis by the
graph of the polynomial $p$.
## Factoring polynomials
Among numerous others, there are two common ways of representing a
non-zero polynomial:
* expanded form, as in $a_n x^n + a_{n-1}x^{n-1} + \cdots a_1 x + a_0, a_n \neq 0$; or
* factored form, as in $a\cdot(x-r_1)\cdot(x-r_2)\cdots(x-r_n), a \neq 0$.
The latter writes $p$ as a product of linear factors, though this is
only possible in general if we consider complex roots. With real roots
only, then the factors are either linear or quadratic, as will be
discussed later.
There are values to each representation. One value of the expanded
form is that polynomial addition and scalar multiplication is much easier than in factored form. For example, adding polynomials just requires matching
up the monomials of similar powers. For the factored form, polynomial multiplication is much easier than expanded form. For the factored form it is easy to read off *roots* of the polynomial (values of $x$ where $p$ is $0$), as
a product is $0$ only if a term is $0$, so any zero must be a zero of
a factor. Factored form has other technical advantages. For example,
the polynomial $(x-1)^{1000}$ can be compactly represented using the
factored form, but would require ``1001`` coefficients to store in expanded
form. (As well, due to floating point differences, the two would
evaluate quite differently as one would require over a ``1000`` operations
to compute, the other just two.)
Translating from factored form to expanded form can be done by
carefully following the distributive law of multiplication. For
The `SymPy` function `expand` will perform these algebraic
manipulations without fuss:
```julia;
expand((x-1)*(x-2)*(x-3))
```
Factoring a polynomial is several weeks worth of lessons, as there is
no one-size-fits-all algorithm to follow. There are some tricks that
are taught: for example factoring differences of perfect squares,
completing the square, the rational root theorem, $\dots$. But in
general the solution is not automated. The `SymPy` function `factor`
will find all rational factors (terms like $(qx-p)$), but will leave
terms that do not have rational factors alone. For example:
```julia;
factor(x^3 - 6x^2 + 11x -6)
```
Or
```julia;
factor(x^5 - 5x^4 + 8x^3 - 8x^2 + 7x - 3)
```
But will not factor things that are not hard to see:
```julia;
x^2 - 2
```
The factoring $(x-\sqrt{2})\cdot(x + \sqrt{2})$ is not found, as
$\sqrt{2}$ is not rational.
(For those, it may be possible to solve to get the roots, which
can then be used to produce the factored form.)
### Polynomial functions and polynomials.
Our definition of a polynomial is in terms of algebraic expressions
which are easily represented by `SymPy` objects, but not objects from
base `Julia`. (Later we discuss the `Polynomials` package for representing polynomials. There is also the `AbstractAlbegra` package for a more algebraic treatment of polynomials.)
However, *polynomial functions* are easily represented by `Julia`, for
example,
```julia;
f(x) = -16x^2 + 100
```
The distinction is subtle, the expression is turned into a function
just by adding the `f(x) =` preface. But to `Julia` there is a big
distinction. The function form never does any computation until after a value
of $x$ is passed to it. Whereas symbolic expressions can be
manipulated quite freely before any numeric values are specified.
It is easy to create a symbolic expression from a function - just
evaluate the function on a symbolic value:
```julia;
f(x)
```
This is easy - but can also be confusing. The function object is `f`,
the expression is `f(x)` - the function evaluated on a symbolic
object. Moreover, as seen, the symbolic expression can be evaluated
using the same syntax as a function call:
```julia;
p = f(x)
p(2)
```
For many uses, the distinction is unnecessary to make, as the many functions will work with any callable expression. One such is `plot` -- either
`plot(f, a, b)` or `plot(f(x),a, b)` will produce the same plot using the
