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

`Pkg.add("SymPy")`

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

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

`SymPy`

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

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

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}`

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

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

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`

.)

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

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

Find

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

Compute

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

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:

(This is an example in Gruntz's thesis from which `SymPy`

's `limit`

function is based upon.)

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

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

```
L
─
2
```

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 `subs`

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

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.

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?