Symbolic math with julia

Pkg.add("SymPy")

About

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")

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

(h,y)

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

   2
 -x 
 ───
  2 
ℯ   

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

 -1/2
ℯ    

The subs function has an alternate syntax:

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

 -1/2
ℯ    

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.

equations

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

1

One can factor polynomials:

factor(x^2 + 2x + 1)

       2
(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}:
 -30.0000000000000
 30.0000000000000 

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])

{"x"=>0,"y"=>1}

Practice

Question

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?

4x^3 + 3x^2 + 2x + 1

1x^3 + 2x^2 + 3x + 4

(1+2+3+4)x

(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.)

Question

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

Question

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

Same as the picture

two answers -5 and 5

The value 5

Question

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\).

Question

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)

Practice

Question

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

Limits

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)

1

We can do other similar questions:

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

1/2

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

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

0

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)

    9
10⋅x 

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)

    9
10⋅x 

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)

    9
10⋅x 

This limit matches the chain rule:

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

      9          
10⋅sin (1)⋅cos(1)

Practice

Question

Find

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

Question

Compute

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

Question

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?

f'(x)

f''(x)

f'''(x)

Question

Define a symbolic quantity a. Compute the limit

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

What is the answer?

1

a

a^2

No limit exists

Question

Let

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:

0

1

e = 2.7182818284590455

Inf

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

Question

Let

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.)

Derivatives

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)

Extrema

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)        

L
─
2

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

Practice

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

1/sqrt(-x^2 + 1)

1/(x^2 + 1)

1/x

Questions

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

Questions

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).)

Question

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\)?

Integration

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⎠
───────────
     2     

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))                                

1/3

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}:
 -0.840896415253715  
 0.840896415253715   
 -0.840896415253715⋅ⅈ
 0.840896415253715⋅ⅈ 

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

Practice

Question

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

3*log(x)

3*log(x)^2/2

-log(x) + log(x*log(x))

-3log(x) + 3log(x*log(x))

Question

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

Questions

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

12.5333, and yes it makes sense.

-8.884, and no it makes no sense as f(x) > 0, but the integral is negative

Inf, and yes it makes sense as the vertical asymptote at 0 makes this a problem integral.