sympy

Computer Algebra Systems

Computer Algebra Systems (CAS), allow the user to perform symbolic calculations. What we have seen so far were only numerical calculations, because when evaluating a formula each variable had a value assigned to it so that it could return a numerical result. For example suppose we had a function like this one: $$f(x)=3x^3+4x^2-1$$ and we wanted the computer to calculate the derivative for us, keeping x as a symbol. Well, this is exactly what Computer Algebra Systems are for. In Python the most widely used CAS is SymPy. It can store formulas and variables as symbolic expressions, and perform symbolic operations on them. Given $\mathbf{f(x)}$, SymPy would be able to calculate $\mathbf{f'(x)}$ as an expression, returning $$f(x)=9x^2+8x$$ SymPy is a very vast and powerful library, and in this book we will only go through its most useful and powerful functionalities. At the end of this chapter you will have at your disposal a solid set of tools to perform symbolic calculations.

Importing the library

As we did when we first imported NumPy in the previous chapter, we're going to import SymPy specifying an abbreviation to for it. In an empty cell write the following code:

import sympy as sp
sp.init_printing(use_latex="mathjax")
display(2+sp.sqrt(3))
display(sp.sqrt(32))

$\sqrt{3}+2$
$4\sqrt{2}$

 In line 1 we import the library, and in line 2 we tell Jupyter how to print SymPy's expressions. More specifically, we want the expressions to be printed in LaTeX format, and we want to use Mathjax to render LaTeX. Lines 3 and 4 use the command display to print two SymPy expressions. We could use \textbf{print} instead of \textbf{display}, but then the output wouldn't be as pretty and easy to read. Notice how SymPy automatically simplifies $\sqrt{32}$ as $4\sqrt{2}$.

 

NOTE:

LaTeX  is a typesetting language, and is in fact the lingua franca for creating beautiful expressions and scientific documents. Don't worry if you have never used it, for now it's just a way to display prettier outputs in the cells.

 

 

Symbols

SymPy lets the user define symbols that will be used during symbolic calculations. You can think of them as undefined variables, created using the command symbols()

a, b=sp.symbols("a b")
expr1=a+b
expr2=expr1**2
display(expr1)
display(expr2)

$a+b$
$(a+b)^2$

 In line 1 we create two symbols, a and b, by passing the symbols() function a string that contains the names of the variables. Each symbol in the string is separated by a space and corresponds to one of the variables on the left of the equal sign. From this point on, the variables a and b will be symbolic, and any operation that involves them will return a symbolic result as well.

The .subs command

Suppose that you have created a symbolic expression, and you want to substitute its symbols with numbers. This kind of operation is done using the subs() method:

x, y=sp.symbols("x y")
f=x**2+y**2+sp.Rational(1,3)*x*y
display(f)
res1=f.subs(x,1)
display(res1)
res2=f.subs([(x,1), (y,5)])
display(res2)
display(res2.evalf(5))

$x^2+\cfrac{xy}{3}+y$
$y^2+\cfrac{y}{3}+1$
$\cfrac{83}{3}$
$27.667$

 After creating the symbols x and y, we define the expression f as seen in line 2. The function sp.Rational is used to specify a rational number, in this case $1/3$. Now we want to substitute $1$ to x: in order to do this we use the .subs() method, specifying the variable that we want to substitute and its value. This operation results in a new expression that is stored in res1. Note that when res1 is printed out in line 5, x has been substituted with $1$. In line 5 we do the same thing again, except this time we substitute both x and y. This is done by passing the .subs() command a list of tuples containing the variables and the values in pairs. Finally we use the .evalf() command to collapse the whole expression in one single float value. The argument of .evalf() specifies the number of digits of the output.

Calculus

Derivatives and integrals are very common operations in engineering, so let's see how SymPy deals with them:

x=sp.symbols("x")
expr=x**2+sp.Rational(1,2)*x-5
display(expr)
derivative=sp.diff(expr,x)
display(derivative)
integral1=sp.integrate(expr,x)
display(integral1)
integral2=sp.integrate(expr,(x, 0, 3))
display(integral2)

$x^2+\cfrac{x}{2}-5$
$2x+\cfrac{1}{2}$
$\cfrac{x^3}{3}+\cfrac{x^2}{4}-5x$
$-\cfrac{15}{4}$

 The code is pretty self-explanatory. To perform the derivative of expr with respect of x, simply pass them as arguments to to sp.diff(), as shown in line 4. The indefinite integral is performed using sp.integrate(), and the arguments are the same as sp.diff(). To calculate the definitive integral you need to group the variable you wish to integrate with respect to with the two extremes of integration in a tuple. Line 8 show an example of this kind of operation, where we integrate expr with respect to x between $0$ and $3$.

Solvers

Solvers are one of the most powerful aspects of SymPy.They are used to solve various types of equations, something that engineers often need to do while designing.

Equations

In SymPy, equations are defined using the Eq() command. These equations can then be used as inputs for the various solvers that SymPy provides in order to find the solutions. Note that it's wrong to specify equations using $=$ or $==$.

x=sp.symbols("x")
my_eq=sp.Eq(x**2-4,0)
display(my_eq)
sp.solveset(my_eq,x)

$x^2-4=0$
${-2,2}$

 In the previous code snippet we create an equation called my_eq,
and then we pass it to solveset(), which is the main solver of SymPy.

NOTE:

It is often superflous to use Eq() to define an equation. Most of the times it's faster to write the equation as $f(x)=0$ and pass $f(x)$ straight to the solver. By default, if the solver recieves an expression as input, it will automatically add $=0$ at the end of it, and solve the equation accordingly.

The solveset command

The main function that SymPy uses to solve algebraic equations is solveset, to which the user can pass an equation written in the form of an Eq instance or an expression that will automatically assumed to be equal to 0. The next code snippet gives two examples of its usage:

x=sp.symbols("x") res1=sp.solveset(sp.cos(x)*2, x) display(res1) res2=sp.solveset(sp.Rational(2,5)*x+3, x) display(res2)
display(res2.args[0])

$\left\{2n\pi+\cfrac{\pi}{2}|n\in\mathbb{Z}\right\}\cup\left\{2n\pi+\cfrac{3\pi}{2}|n\in\mathbb{Z}\right\}$
$\left\{-\cfrac{15}{2}\right\}$
$-\cfrac{15}{2}$

 Line 2 shows how to use solveset() to solve the simple equation $2cos(x)=0$. Notice how we specify the cosine we use the SymPy function $sp.cos()$, and how we don't need to add $=0$ at the end of the expression. The lines 5 through 6 show how to solve an equation and retrieve the resoult value using the .args argument. But why use .args instead of writing, say, res2[0]? Because the results of solveset() aren't lists, but rather data structures defined by the library itself. In the example above, the result what SymPy calls a FiniteSet, which isn't in itself an iterable object. To acces the values stored in a FiniteSet, you need to use the .args, and pass it the index of the result you want to retrieve.

Systems of equations: linsolve and nonlinsolve

SymPy can of course also solve systems of equations. It gives the user the ability to solve linear equations using linsolve(), and nonlinear equations using nonlinsolve(). The syntax is very similar to  solveset, only you need to pass a list of equations instead of just one and a tuple containing the unknown terms. Let's see how to use these two modules:

x, y=sp.symbols("x y")
res1=sp.linsolve([x+y-4, x-y-9], (x,y))
display(res1.args[0])
res2=sp.nonlinsolve([2*x**2+3*x-y-1, 3*x+2*y-5],(x,y))
display(res2)

$\left(\cfrac{13}{2},-\cfrac{5}{2}\right)$
$\left\{\left(-\cfrac{9}{8}+\cfrac{\sqrt{193}}{8},\; -\cfrac{3\sqrt{193}}{16}+\cfrac{67}{16}\right),\;\left(-\cfrac{\sqrt{193}}{8}-\cfrac{9}{8},\; \cfrac{3\sqrt{193}}{16}+\cfrac{67}{16}\right)\right\}$

And that is pretty much all there is to know to start using SymPy. If you followed this tutorial carefully, you should now be able to use it in your own projects.