# Intro to SymPy for Calculus (Jupyter)#

By Raheem Mir

Python has many powerful librairies that add convenience and functionality. One of them being SymPy, a feature rich library designed for symbolic mathematics. We can use it to simplify algebraic expressions, solve equations, and perform calculus tasks like evaluating limits or computing integrals. SymPy gives us exact symbolic answers instead of numerical approximations, handling math like we would on paper, which makes it a great tool for intuitively exploring calculus concepts.

## Getting Started With SymPy#

To begin using SymPy we start by “importing” the library:

```
import sympy as sy
sy.init_printing()
```

`import sympy as sy`

loads in the SymPy library for us to use. We give it the shorthand name `sy`

, so every time we call something from SymPy, we prefix it with `sy.`

. This is a good way to keep things organized when we are working with multiple libraries.

`sy.init_printing()`

calls SymPy’s `init_printing()`

function, which gets the output to be displayed as nicely formatted mathematics.

A good next step is telling SymPy the variables we would like to use with the `symbols()`

function:

```
x, y = sy.symbols('x y') # this makes x and y symbols
```

When creating a symbol using the `symbols()`

function, we can specify additional arguments that can modify SymPy’s assumptions when it’s dealing with the symbol. These assumptions determine how SymPy handles a symbol in situations like square roots, solving, simplifications, and so on. For instance, `positive=True`

tells SymPy that a particular symbol is a positive number, while using `real=True`

indicates that a symbol represents a real number.

A common task in SymPy is creating symbolic expressions to represent functions. The syntax is no different from assigning a variable.

Let’s represent \(f(x) = x^2 - 2x + 1\) in SymPy!

```
f = x**2 -2*x + 1 # defining a symbolic expression
```

```
f # outputting f(x)
```

It is important to note that we are not actually creating a Python function here, which would be defined on an input/output basis, rather we are creating a symbolic representation of a function, that can be easily manipulated with algebra and/or calculus. As such, trying something like `f(5)`

would return an error. To substitute values into `f`

, we would use the `subs()`

function:

```
f.subs(x, 4) #substituting 4 in for x
```

SymPy also gives us a function for factoring a polynomial, `factor()`

:

```
sy.factor(f)
```

SymPy implements many of the common mathematical functions as well:

```
sy.cos(x)
```

```
sy.tan(x)
```

```
sy.pi, sy.exp(x) # constants pi and e
```

```
sy.Abs(x) # absolute value
```

## Limits#

With SymPy, we can use its `limit()`

function to compute and evaluate limits. The syntax is as follows:

```
sy.limit(expression, variable, value)
```

`expression`

is the SymPy representation of the function whose limit we are trying to find.

`variable`

represents the variable that is approaching a specific value, and must be a SymPy symbol.

`value`

represents the value the variable is approaching. If we are dealing with limits at infinity, we can use `sy.oo`

and `-sy.oo`

to dentote positive and negative infinity.

To compute one-sided limits, we can add a fourth argument to the function, either a `'+'`

or `'-'`

for a right or left hand limit.

Let’s evaluate \(\displaystyle \lim_{x\to 0}\frac{\sin(x)}{x}\) with the `limit()`

function:

```
f = sy.sin(x) / x # representing the function
f
```

```
sy.limit(f, x, 0) # calling the limit() function
```

Since we set our symbolic representation equal to a variable, `f`

, we could have also called `limit()`

like this:

```
f.limit(x, 0)
```

Note that we can also directly pass in a SymPy expression to the `limit()`

function instead of using a variable.

Let’s use the `limit()`

function to find the value of \(\displaystyle \lim_{x\to 0}\dfrac{\cos(x) - 1}{x}\).

```
sy.limit((sy.cos(x) - 1)/x, x, 0)
```

Let’s try adding a fourth argument in our call to `limit()`

, and work with one-sided limits:

\(\displaystyle \lim_{x\to 0^-}\dfrac{1}{x}\)

```
sy.limit(1/x, x, 0, '-')
```

\(\displaystyle \lim_{x\to 0^+}\dfrac{1}{x}\)

```
sy.limit(1/x, x, 0, '+')
```

Now let’s explore a limit where \(x\) is approaching \(\infty\):

\(\displaystyle \lim_{x\to \infty}\dfrac{x^2 - 1}{3 - x}\)

```
sy.limit((x**2 - 1) / (3 - x), x, sy.oo)
```

## Derivatives#

To take derivatives we can use SymPy’s `diff()`

function:

```
sy.diff(function,variable, ...)
```

`function`

is the function we are differentiating. Note that this is a SymPy represenation of a function, not a traditional Python function.

`variable`

is the variable with respect to which we are differentiating. This must be defined as a SymPy symbol.

`...`

represents the possibility of adding further arguments, namely for higher order derivatives.

Given the function \(f(x) = 4x^2 + 6x - 2\) let’s find \(f'(x)\) using the `diff()`

function:

```
f = 4*x**2 + 6*x - 2 # defining a symbolic representation of f(x)
f
```

```
f_prime = sy.diff(f, x) # differentiating f with respect to x
f_prime
```

Like before, we have some freedom with the syntax when calling the `diff()`

function:

```
f_prime = f.diff(x)
f_prime
```

Now let’s differentiate \(g(x) = x^2\sin(x)\) with `diff()`

:

```
g = sy.sin(x)*x**2
g_prime = g.diff(x)
g_prime
```

The `diff()`

function can also compute higher order derivatives, as we simply add more arguments to our function call! This is not unlike how we add a tick (or prime) to our notation when taking the next higher order derivative when working on paper.

Given \(f(x) = \sin(x)\) let’s find its second and third derivatives:

```
f1 = sy.sin(x)
fpp = f1.diff(x, x) # second derivative of f(x)
fppp = f1.diff(x, x, x) # third derivative of f(x)
print(fpp, fppp)
```

Instead of writing out a particular variable multiple times, we can also compute higher order derivatives like this:

```
fpp = f1.diff(x, 2)
fppp = f1.diff(x, 3)
print(fpp, fppp)
```

## Plotting#

SymPy also provides functionality for visualizing / plotting functions:

We simply use the `plot()`

function, like so:

```
sy.plot(f, domain, ...)
```

`f`

represents the function (or functions) we are visualizing, which can be written explicitly or passed in as a variable.

`domain`

refers to the interval we are the plotting the function(s) over, the default is \([-10, 10]\).

`...`

alludes to the many optional arguments we can use to change the plot’s appearance and behaviour.

An example would be the argument `show=False`

, which stops the graph from being displayed until we use the command, `.show()`

.

Let’s plot \(f(x) = x^2\) over the interval \([-5, 5]\):

```
p = sy.plot(x**2, (x, -5, 5)) # notice how we input the domain
```

Interestingly, we don’t need to type `p`

on the next line for the plot to output, like we’d have to when using other SymPy functions. Also, we save our plot to a variable as it’s more efficient when we are manipulating/adjusting an existing plot, or wanting more control over customization.

Now let’s plot \(f(x) = \sin(5x)\) on the interval \([0, 2\pi]\) and add an optional argument to change the color of the line!

```
p1 = sy.plot(sy.sin(5*x), (x, 0, sy.pi*2), line_color='purple') # using SymPy's built in functions for sin and pi
```

Here we plot \(f(x) = \ln(x)\) on the interval \([-10, 10]\). Like with the trigonomtric functions, we’ll be using SymPy’s built in function for natural logarithms, `sy.ln()`

. Remember that \([-10, 10]\) is the default domain so we don’t need to write it.

```
p2 = sy.plot(sy.ln(x), line_color="green")
```

SymPy’s plotting functionality also gives us the ability to plot multiple functions together, allowing us to combine different plots.

Let’s illustrate this with an example, by plotting the following equations:

\(y = x\), \(y = x^2\) , \(y = x^3\), and \(y = x^4\) over the interval \([0, 2]\):

```
multi_plot = sy.plot(x, x**2, x**3, x**4, (x, 0, 2), legend=True) # adding a legend to our plot
```

SymPy also lets us combine existing plots, with the `extend()`

function:

```
plot_a = sy.plot(-x**2, show=False)
plot_b = sy.plot(-x, show=False)
plot_a.extend(plot_b)
plot_a.show()
```

Notice how we added the `show=False`

parameter to our plots, so they do not output until we use `.show()`

.

## Integration#

To compute integrals we can use SymPy’s `integrate`

function:

For indefinite integrals (antiderivatives), the syntax is quite similar to the `diff()`

function:

```
sy.integrate(expression, variable)
```

`expression`

is a SymPy representation of the function we want to integrate.

`variable`

is the variable of integration and must be defined as a SymPy symbol.

One thing to note is that SymPy does not include a constant of integration (the “\(+ c\)”).

For definite integrals the syntax changes a bit:

```
sy.integrate(expression, (variable, lower_limit, upper_limit))
```

`expression`

is the same as with indefinite integrals, a SymPy representation of the function we want to integrate.

`(variable, lower_limit, upper_limit)`

is known as a tuple in Python, a container of more arguments, where:

`variable`

is the same as with indefinite integrals, the variable of integration that must be defined as a SymPy symbol.

`lower_limit`

and `upper_limit`

are the lower and upper bounds of integration.

Let’s evaluate the indefinite integral \(\int \sin(x) \ dx\) with the `integrate()`

function:

```
indef_integral = sy.integrate(sy.sin(x), x)
indef_integral
```

Now let’s evaluate the indefinite integral \(\int x^3 \ dx\) with `integrate()`

:

```
indef_integral1 = sy.integrate(x**3, x)
indef_integral1
```

Next, let’s look at evaluating definite integrals with the `integrate()`

function:

\(\int_{-3}^{3} x^7 - x^3 \, dx\)

```
def_integral = sy.integrate(x**7 - x**3, (x, -3, 3))
def_integral
```

Now, let’s evaluate the definite integral \(\int_{1}^{2} e^x \, dx\)

```
def_integral1 = sy.integrate(sy.exp(x), (x, 1, 2))
def_integral1
```

After computing \(\int_{1}^{2} e^x \, dx\) symbolically (the result is kept as an exact value), we can also find its numerical value with SymPy’s `evalf()`

function:

```
numerical_result = def_integral1.evalf()
numerical_result
```