Symbolic math with julia

The julia language bills itself as "fresh approach to technical computing." By saying "fresh" the implication is that there exists many older approaches to technical computing. Indeed there are. For mathematical areas there are three different philosophies for computing: symbolic, numeric, and general purpose. The symbolic approach is the domain of Computer Algebra Systems (CAS), and is exemplified by very comprehensive programs like Mathematica, Maple and the open-source alternative Sage. The numeric approach is the domain of tailored programming languages such as MATLAB and R. A general purpose approach would be to leverage a widely used programming language, such as Python or Haskell, for a specific use. The julia language is an alternative approach to MATLAB or R for numerical computation.

One strength of julia is how well it plays with others. This is leveraged in the SymPy package for julia to provide a symbolic math interface through a connection to Python and its SymPy package via julia's PyCall package. We see in this project how this additional functionality affords an alternative approach to performing calculus problems.

The SymPy package for julia is an add on. If not already installed, one can install it through the command


In order for this to work, you should have a working version of Python on you computer and the SymPy package for Python installed for that.


The SymPy program extends julia by providing a type for symbolic expressions. Such an expression is encapsulated by a symbolic variable x instantiated through:

using SymPy
x = Sym("x")

(We typeset symbolic expressions differently in this project.)

The "x" on the right-hand side is a character argument to the Sym constructor which returns a symbolic object stored as x:

x |> typeof
## Sym

That was painless. While we are here, lets define symbolic objects h and y:

h, y = Sym("h", "y")        # two at once

Most of the typical math functions have been overloaded to work with these symbolic expressions: the functions accept a symbolic expression and return a newly computed one. For example, the expression x^2 -2x +2 when evaluated becomes a new symbolic expression:

x^2 - 2x + 2 |> typeof
## Sym

Similarly, when working with functions a symbolic expression is returned:

f(x) = exp(-x^2/2)             ## a julia function
f(x)                           ## takes a symbolic object and returns a new one

This shows that function object f will map a symbolic object to another symbolic object.

In SymPy, the subs function allows you to evaluate a symbolic expression for a given value. For example,

subs(f(x), x, 1)        ## set x equal to 1

The subs function has an alternate syntax:

f(x) |> replace(x, 1)

The output of subs is still a symbolic object. To pull the numeric value back into julia it can be coerced through integer (if applicable) or float:

f(x) |> replace(x, 1) |> float
## 0.6065306597126334

For the most part, one can work with symbolic expressions without pulling them back into julia expressions until needed.


SymPy can do much of the basic tasks learned during algebra: simplification, factoring and solving equations. Just a few new commands are needed.

Basic algebra

SymPy does some automatic simplification of expressions. For example, terms are combined

(x+1) + (x+2) + (x+3)
3⋅x + 6

However, not everything is as simplified as possible. The simplify function can be called to (sometimes) go further:

sin(x)^2 + cos(x)^2     # not simplified
   2         2   
sin (x) + cos (x)
simplify(sin(x)^2 + cos(x)^2)   # 1

One can factor polynomials:

factor(x^2 + 2x + 1)
(x + 1) 

SymPy will leave terms in factored form. To multiply them out, can be requested through the function expand:

expand( (x+1)*(x+2)*(x+3) )
 3      2           
x  + 6⋅x  + 11⋅x + 6

Solving equations

We can now see how to use symbolic expressions. First we discuss solving for when an expression is 0.

For polynomials, the real_roots function will find the real roots, returning them as julia numbers:

real_roots(x^2 + x - 1) 
2-element Array{Sym,1}:
  1   ╲╱ 5 
- ─ + ─────
  2     2  
  ╲╱ 5    1
- ───── - ─
    2     2

If the formatting is hard to read, there is a way to get a simplified printout:

real_roots(x^2 + x - 1) |> jprint # simplified output
2-element Array{ASCIIString,1}:
 "-1/2 + sqrt(5)/2"
 "-sqrt(5)/2 - 1/2"

More generally, the solve command can be used to solve equations -- and not just polynomial equations -- symbolically. The same problem above becomes:

solve(x^2 + x -1, x)
2-element Array{Sym,1}:
  1   ╲╱ 5 
- ─ + ─────
  2     2  
  ╲╱ 5    1
- ───── - ─
    2     2

Again, these are symbolic answers. We can coerce values back into julia with float to get numeric representations:

solve(x^2 + x -1, x) |> float
## [0.6180339887498949,-1.618033988749895]

Trigonometric functions have infinitely many solutions, with the sin function SymPy solves only within the range \([-90, 90]\) degrees (the range of the asin function).

solve(sin(x)^2 - (1/2)^2) * (180 /pi)
2-element Array{Sym,1}:

To solve an expression in another variable, we specify it through the second argument:

out = solve(x^2 + y^2 - 1, y)
2-element Array{Sym,1}:
   ╱    2     
-╲╱  - x  + 1 
  ╱    2     
╲╱  - x  + 1    

This returns a vector of two symbolic answers, as the expression being solved is a quadratic. We can then evaluate these two values at a point, say \(x=1/2\), as before:

out |> replace(x, 1/2) |> float
## [-0.8660254037844386,0.8660254037844386]

You can even do systems of equations. For this you specify the system and the variables to solve for as vectors:

solve([x + y - 1, x - y - (-1)], [x,y])



What is the coefficient of \(x^3\) in \((x-1)(x-2)(x-3)(x-4)(x-5)\)?

Question Horner's method

Expand \((((1x + 2)x + 3)x + 4)\). What do you get?

(Horner's method is an alternate way of evaluating \(a_n x^n + a_{n-1}x^{n-1} + \cdots a_1 x + a_0\) expressed in terms of the coefficients, which are 1,2,3,4 in the example.)


Find the largest root of \(x^3 - 6x^2 + 11x - 6\).


Embedded Image

Solve for \(x\) using solve(x^2 - (3^2 + 4^2)) not the picture. What do you get


The Soviet historian I. Y. Depman claimed that in 1486, Spanish mathematician Valmes was burned at the stake for claiming to have solved the quartic equation. Here we don't face such consequences.

Find the largest root of \(x^4 - 10x^3 + 32x^2 - 38x + 15\).


Solve \(e^x = x^4\) using julia. (You will need to convert to numeric with float.)

Graphing expressions

Plotting symbolic expressions can be done directly without bringing them back into julia. For example, here we plot the function \(f(x) = \sin(x^4) - \cos(x^4)\) over \([0,1]\):

f(x) = sin(x^4) - cos(x^4)
plot(f(x), 0, 1)

This is very similar to how we can plot a function object, f, using plot, but note that f(x) is a symbolic expression, not a function object:

(typeof(f(x)), typeof(f))
## (Sym,Function)

Graphing more than one function at a time is possible. We just specify an array of symbolic expressions. For example, to graph a function and its tangent line at a point we can do this.

f(x) = sqrt(x);
fp = diff(f(x));        # diff finds derivative, fp an expression (not function)
c = 2;
m = subs(fp, x, c);     # at c=2
plot([f(x), f(c) + m*(x-c)], 1, 3)



Plot the expression x^4 - 3x^3 + 3x^2 - 3x +2 over \([-1, 3]\). How many real roots are there?


The Limit function can find the limit of an expression. Let's see how well it does. The basic question

\[ \lim_{x \rightarrow c} f(x) = L \]

has three inputs: A function the limit is being taken of, a dummy variable \(x\), and a value where the limit is taken, \(c\). The output is the limit, when it exists. The Limit function is similar. Here we find an old classic:

limit(sin(x)/x, x, 0)

We can do other similar questions:

limit((1-cos(x))/x^2, x, 0)

Limits can be taken at infinity as well. We specify that value using oo:

limit(sin(x)/x, x, oo)

We can even compute derivatives using limits. Here we do so symbolically:

h = Sym("h")
f(x) = x^10
limit( (f(x+h) - f(x))/h, h, 0)

We can see that some of the more complicated formulas for derivatives give the same answer. Here we see the central difference gives the same answer:

central_difference(f, x, h) = ( f(x+h) - f(x-h) ) / (2h)
a = limit(central_difference(f, x, h), h, 0)

Even this more complicated expression works as expected:

central_4th_difference(f, x, h) = (-f(x + 2h) + 8f(x+h) - 8f(x-h) + f(x-2h))/(12h)
limit(central_4th_difference(f, x, h), h, 0)

This limit matches the chain rule:

g(x) = sin(x)
limit( (f(g(x+h)) - f(g(x)))/h, h, 0) |> replace(x, 1)
10⋅sin (1)⋅cos(1)




\[ \lim_{x \rightarrow 4} \frac{(3 - \sqrt(x + 5))}{(x-4)} \]



\[ lim_{x \rightarrow 1} \frac{x^{1/4} - 1}{x^{1/3} - 1}. \]


Define a symbolic variable h and let \(f(x) = sin(x)\):

h = Sym("h");
f(x) = sin(x);

Compute the following:

\[ \lim_{h->0} \frac{f(x+h) - 2f(x) + f(x-h)}{h^2} \]

Based on your answer, what do you think the above expression finds for any f?


Define a symbolic quantity a. Compute the limit

\[ \lim_{x \rightarrow 0} \frac{\cos(ax) - 1}{\cos(x) - 1} \]

What is the answer?



f(x) = 1/(x^(log(log(log(1/x)))-1))

We can investigate the limit \(\lim_{x \rightarrow 0} f(x)\) numerically:

xs = (0.1).^(2:10);
fxs = [f(x) for x in xs];
[xs fxs]
## 9x2 Array{Any,2}:
##  0.01      0.0702822
##  0.001     0.0947694
##  0.0001    0.155102 
##  1.0e-5    0.293154 
##  1.0e-6    0.619862 
##  1.0e-7    1.43553  
##  1.0e-8    3.58731  
##  1.0e-9    9.56746  
##  1.0e-10  27.0054   

Rather than do that, use limit to compute the limit. The answer is:

(This is an example in Gruntz's thesis from which SymPy's limit function is based upon.)



f(x) = log(log(x*exp(x*exp(x)) + 1)) - exp(exp(log(log(x)) + 1/x))

This function has a limit at \(\infty\). It can't be found numerically though (well, not easily). You can verify this by trying to evaluate the function at \(10\), let alone for really large values. However, it can be computed symbolically. What is the answer?

(Again, this is from Gruntz.)


We just saw we can take derivatives with a limit, though just as with pen and paper, it is better to use the rules of derivatives. These rules are implemented by the diff operator, which takes a symbolic expression, a variable to differentiate in (the default is x) and an optional integer for the number of derivatives and returns a symbolic derivative.

For example

f(x) = exp(exp(x))
diff(f(x))             ## not f, but the symbolic expression f(x)
diff(f(x), x)                  ## optional x is necessary if other symbols involved
diff(f(x), x, 2)               ## Finds f''(x), not f'(2).
      ⎛ x⎞       ⎛ x⎞
 2⋅x  ⎝ℯ ⎠    x  ⎝ℯ ⎠
ℯ   ⋅ℯ     + ℯ ⋅ℯ    

Here we find and plot the derivative of \(x^x\) avoiding the asymptote at \(x=0\).

f(x) = x^x
plot([ f(x), diff(f(x)) ], 1/10, 2)

Here we plot \(f(x)\) and it's tangent line at \(c=1\):

f(x) = x^x; c = 1
m = diff(f(x)) |> replace(x, c);
plot([ f(x), f(c) + m * (x - c)], .5, 1.5)


We can combine solve with diff to find extrema of a differentiable function over a closed interval, by locating the critical points where \(f'(x) = 0\).

For example, consider the problem enclosing the maximum amount of area with 3 sides of a rectangle where the length of the 3 sides is fixed.

For this, we have

\[ L = 2x + y, \quad A = xy, \quad\text{or}\quad A(x) = x \cdot (L-2x) \]

We can solve this with \(L\) as a symbolic value, by looking at critical points of \(A\) or when \(A'(x) = 0\):

L = Sym("L")
A = x*y
A = subs(A, y, L - 2x)
out = solve(diff(A, x), x)[1]   ## solve returns an array, we need its first component
subs(L - 2x, x, out)        

So \(x = L/4\) and \(y = L/2\).


Find the derivative of \(\tan^{-1}(x)\). What do you get?


Let \(f(x) = \exp(-\cos(x))\). Evaluate \(f'(3)\). (You'll need to substitute 3 for x and convert to floating point.)


Find the single critical point of \(f(x) = 1/x^2 + x^3\). (Solve, will return many answers but only one real one. You can also try nsolve(diff(f(x)), 1).)


Let \(f(x) = \tan(x)\). Newton's method finds the zero of the tangent line \(f(c) + f'(c)(x-c)\). You can do this with julia via:

f(x) = tan(x)
c = Sym("c");
solve(f(c) + diff(f(c)) * (x - c), x)
## 1-element Array{Sym,1}
## ⎡    sin(2⋅c)⎤
## ⎢c - ────────⎥
## ⎣       2    ⎦

What do you get when \(c = \pi/4\)?


The SymPy program can do symbolic integration, though it is not as effective at it, as say Mathematica, it can do many things. The standard function is integrate.

One can find general antiderivatives:

x, a = Sym("x", "a")
f(x) = cos(x) - sin(x)
integrate(f(x), x)
sin(x) + cos(x)

One can parameterize an integral using a constant:

integrate(1/cos(x + a), x)
     ⎛   ⎛a   x⎞    ⎞      ⎛   ⎛a   x⎞    ⎞
- log⎜tan⎜─ + ─⎟ - 1⎟ + log⎜tan⎜─ + ─⎟ + 1⎟
     ⎝   ⎝2   2⎠    ⎠      ⎝   ⎝2   2⎠    ⎠

Integrate can do many of the integrals thrown at it.

integrate(x^3 * (1-x)^4, x)
 8      7           5    4
x    4⋅x     6   4⋅x    x 
── - ──── + x  - ──── + ──
8     7           5     4 
integrate(x/(x^2+1), x)
   ⎛ 2    ⎞
log⎝x  + 1⎠

Definite integrals

By specifying a range to integrate over, the definite integral \(\int_a^b f(x) dx\) can be found.

integrate(x^2, (x, 0, 1))                                

More interestingly, we can find the median value of an integral by solving with the resulting symbolic answer from integrate:

f(x) = 4x^3
b = Sym("b")
eq = integrate(f(x), (x,  0, b)) - 1/2 * integrate(f(x), (x, 0, 1))
solve(eq, b)
4-element Array{Sym,1}:

We can see that the answer is \(0.840896415253715\), though the mathematical expression has 3 other possibilities for solution.



Integrate \(f(x) = 2/(x \cdot \log(x))\). What do you get:


Integrate \(f(x) = (2x+5)(x^2 + 5x)^7\) from \(0\) to \(1\).


Integrate f(x) = sqrt(4 - sqrt(x)) over \([0,9]\). What value do you get? Does it make any sense?