# DOING MATH WITH PYTHON Use Programming to Explore Algebra, Statistics, Calculus, and More! (2015)

### 4

Algebra and Symbolic Math with SymPy

The mathematical problems and solutions in our programs so far have all involved the manipulation of numbers. But there’s another way math is taught, learned, and practiced, and that’s in terms of symbols and the operations between them. Just think of all the *x*s and *y*s in a typical algebra problem. We refer to this type of math as *symbolic math*. I’m sure you remember those dreaded “factorize *x*^{3} + 3*x*^{2} + 3*x* + 1” problems in your math class. Fear no more, for in this chapter, we learn how to write programs that can solve such problems and much more. To do so, we’ll use*SymPy*—a Python library that lets you write expressions containing symbols and perform operations on them. Because this is a third-party library, you’ll need to install it before you can use it in your programs. The installation instructions are described in *Appendix A*.

**Defining Symbols and Symbolic Operations**

*Symbols* form the building blocks of symbolic math. The term *symbol* is just a general name for the *x*s, *y*s, *a*s, and *b*s you use in equations and algebraic expressions. Creating and using symbols will let us do things differently than before. Consider the following statements:

>>> **x = 1**

>>> **x + x + 1**

3

Here we create a label, x, to refer to the number 1. Then, when we write the statement x + x + 1, it’s evaluated for us, and the result is 3. What if you wanted the result in terms of the symbol *x*? That is, if instead of 3, you wanted Python to tell you that the result is 2*x* + 1? You couldn’t just write x + x + 1 *without* the statement x = 1 because Python wouldn’t know what x refers to.

SymPy lets us write programs where we can express and evaluate mathematical expressions in terms of such symbols. To use a symbol in your program, you have to create an object of the Symbol class, like this:

>>> **from sympy import Symbol**

>>> **x = Symbol('x')**

First, we import the Symbol class from the sympy library. Then, we create an object of this class passing 'x' as a parameter. Note that this 'x' is written as a string within quotes. We can now define expressions and equations in terms of this symbol. For example, here’s the earlier expression:

>>> **from sympy import Symbol**

>>> **x = Symbol('x')**

>>> **x + x + 1**

2*x + 1

Now the result is given in terms of the symbol *x*. In the statement x = Symbol('x'), the x on the left side is the Python label. This is the same kind of label we’ve used before, except this time it refers to the symbol *x* instead of a number—more specifically, a Symbol object representing the symbol 'x'. This label doesn’t necessarily have to match the symbol either— we could have used a label like a or var1 instead. So, it’s perfectly fine to write the preceding statements as follows:

>>> **a = Symbol('x')**

>>> **a + a + 1**

2*x + 1

Using a non-matching label can be confusing, however, so I would recommend choosing a label that’s the same letter as the symbol it refers to.

**FINDING THE SYMBOL REPRESENTED BY A SYMBOL OBJECT**

For any Symbol object, its name attribute is a string that is the actual symbol it represents:

>>> **x = Symbol('x')**

>>> **x.name**

'x'

>>> **a = Symbol('x')**

>>> **a.name**

'x'

You can use .name on a label to retrieve the symbol that it is storing.

Just to be clear, the symbol you create has to be specified as a string. For example, you can’t create the symbol *x* using x = Symbol(x)—you must define it as x = Symbol('x').

To define multiple symbols, you can either create separate Symbol objects or use the symbols() function to define them more concisely. Let’s say you wanted to use three symbols—*x*, *y*, and *z*—in your program. You could define them individually, as we did earlier:

>>> **x = Symbol('x')**

>>> **y = Symbol('y')**

>>> **z = Symbol('z')**

But a shorter method would be to use the symbols() function to define all three at once:

>>> **from sympy import symbols**

>>> **x,y,z = symbols('x,y,z')**

First, we import the symbols() function from SymPy. Then, we call it with the three symbols we want to create, written as a string with commas separating them. After this statement is executed, x, y, and z will refer to the three symbols 'x', 'y', and 'z'.

Once you’ve defined symbols, you can carry out basic mathematical operations on them, using the same operators you learned in *Chapter 1* (+, -, /, *, and **). For example, you might do the following:

>>> **from sympy import Symbol**

>>> **x = Symbol('x')**

>>> **y = Symbol('y')**

>>> **s = x*y + x*y**

>>> **s**

2*x*y

Let’s see whether we can find the product of x(x + x):

>>> **p = x*(x + x)**

>>> **p**

2*x**2

SymPy will automatically make these simple addition and multiplication calculations, but if we enter a more complex expression, it will remain unchanged. Let’s see what happens when we enter the expression (x + 2)*(x + 3):

>>> **p = (x + 2)*(x + 3)**

>>> **p**

(x + 2)*(x + 3)

You may have expected SymPy to multiply everything out and output x**2 + 5*x + 6. Instead, the expression was printed exactly how we entered it. SymPy automatically simplifies only the most basic of expressions and leaves it to the programmer to explicitly require simplification in cases such as the preceding one. If you want to multiply out the expression to get the expanded version, you’ll have to use the expand() function, which we’ll see in a moment.

**Working with Expressions**

Now that we know how to define our own symbolic expressions, let’s learn more about using them in our programs.

*Factorizing and Expanding Expressions*

The factor() function decomposes an expression into its factors, and the expand() function expands an expression, expressing it as a sum of individual terms. Let’s test out these functions with the basic algebraic identity *x*^{2} - *y*^{2} = (*x* + *y*)(*x* - *y*). The left side of the identity is the expanded version, and the right side depicts the corresponding factorization. Because we have two symbols in the identity, we’ll create two Symbol objects:

>>> **from sympy import Symbol**

>>> **x = Symbol('x')**

>>> **y = Symbol('y')**

Next, we import the factor() function and use it to convert the expanded version (on the left side of the identity) to the factored version (on the right side):

>>> **from sympy import factor**

>>> **expr = x**2 - y**2**

>>> **factor(expr)**

(x - y)*(x + y)

As expected, we get the factored version of the expression. Now let’s expand the factors to get back the original expanded version:

>>> **factors = factor(expr)**

>>> **expand(factors)**

x**2 - y**2

We store the factorized expression in a new label, factors, and then call the expand() function with it. When we do this, we receive the original expression we started with. Let’s try it with the more complicated identity *x*^{3} + 3*x*^{2}*y* + 3*xy*^{2} + *y*^{3} = (*x* + *y*)^{3}:

>>> **expr = x**3 + 3*x**2*y + 3*x*y**2 + y**3**

>>> **factors = factor(expr)**

>>> **factors**

(x + y)**3

>>> **expand(factors)**

x**3 + 3*x**2*y + 3*x*y**2 + y**3

The factor() function is able to factorize the expression, and then the expand() function expands the factorized expression to return to the original expression.

If you try to factorize an expression for which there’s no possible factorization, the original expression is returned by the factor() function. For example, see the following:

>>> **expr = x + y + x*y**

>>> **factor(expr)**

x*y + x + y

Similarly, if you pass in an expression to expand() that can’t be expanded further, it returns the same expression.

*Pretty Printing*

If you want the expressions we’ve been working with to look a bit nicer when you print them, you can use the pprint() function. This function will print the expression in a way that more closely resembles how we’d normally write it on paper. For example, here’s an expression:

>>> **expr = x*x + 2*x*y + y*y**

If we print it as we’ve been doing so far or use the print() function, this is how it looks:

>>> **expr**

x**2 + 2*x*y + y**2

Now, let’s use the pprint() function to print the preceding expression:

>>> **from sympy import pprint**

>>> **pprint(expr)**

x^{2} + 2·x·y + y^{2}

The expression now looks much cleaner—for example, instead of having a bunch of ugly asterisks, exponents appear above the rest of the numbers.

You can also change the order of the terms when you print an expression. Consider the expression 1 + 2*x* + 2*x*^{2}:

>>> **expr = 1 + 2*x + 2*x**2**

>>> **pprint(expr)**

2·x^{2} + 2·x + 1

The terms are arranged in the order of powers of *x*, from highest to lowest. If you want the expression in the opposite order, with the highest power of *x* last, you can make that happen with the init_printing() function, as follows:

>>> **from sympy import init_printing**

>>> **init_printing(order='rev-lex')**

>>> **pprint(expr)**

1 + 2·x + 2·x^{2}

The init_printing() function is first imported and called with the keyword argument order='rev-lex'. This indicates that we want SymPy to print the expressions so that they’re in *reverse lexicographical order*. In this case, the keyword argument tells Python to print the lower-power terms first.

**NOTE**

*Although we used the* *init_printing()* *function here to set the printed order of the expressions, this function can be used in many other ways to configure how an expression is printed. For more options and to learn more about printing in SymPy, see the documentation at**http://docs.sympy.org/latest/tutorial/printing.html*.

Let’s apply what we’ve learned so far to implement a series printing program.

**Printing a Series**

Consider the following series:

Let’s write a program that will ask a user to input a number, *n*, and print this series for that number. In the series, *x* is a symbol and *n* is an integer input by the program’s user. The *n*th term in this series is given by

We can print this series using the following program:

'''

Print the series:

x + x**2 + x**3 + ... + x**n

____ _____ _____

2 3 n

'''

from sympy import Symbol, pprint, init_printing

def print_series(n):

# Initialize printing system with reverse order

init_printing(order='rev-lex')

x = Symbol('x')

➊ series = x

➋ for i in range(2, n+1):

➌ series = series + (x**i)/i

pprint(series)

if __name__ == '__main__':

n = input('Enter the number of terms you want in the series: ')

➍ print_series(int(n))

The print_series() function accepts an integer, n, as a parameter that is the number of terms in the series that will be printed. Note that we convert the input to an integer using the int() function when calling the function at ➍. We then call the init_printing() function to set the series to print in reverse lexicographical order.

At ➊, we create the label, series, and set its initial value as x. Then, we define a for loop that will iterate over the integers from 2 to n at ➋. Each time the loop iterates, it adds each term to series at ➌, as follows:

i = 2, series = x + x**2 / 2

i = 3, series = x + x**2/2 + x**3/3

--*snip*--

The value of series starts off as just plain x, but with each iteration, x**i/i gets added to the value of series until the series we want is completed. You can see SymPy addition put to good use here. Finally, the pprint() function is used to print the series.

When you run the program, it asks you to input a number and then prints the series up to that term:

Enter the number of terms you want in the series: **5**

x^{2} x^{3} x^{4} x^{5}

x + -- + -- + -- + --

2 3 4 5

Try this out with a different number of terms every time. Next, we’ll see how to calculate the sum of this series for a certain value of *x*.

*Substituting in Values*

Let’s see how we can use SymPy to plug values into an algebraic expression. This will let us calculate the value of the expression for certain values of the variables. Consider the mathematical expression *x*^{2} + 2*xy* + *y*^{2}, which can be defined as follows:

>>> **x = Symbol('x')**

>>> **y = Symbol('y')**

>>> **x*x + x*y + x*y + y*y**

x**2 + 2*x*y + y**2

If you want to evaluate this expression, you can substitute numbers in for the symbols using the subs() method:

➊ >>> **expr = x*x + x*y + x*y + y*y**

>>> **res = expr.subs({x:1, y:2})**

First, we create a new label to refer to the expression at ➊, and then we call the subs() method. The argument to the subs() method is a Python *dictionary*, which contains the two symbol labels and the numerical values we want to substitute in for each symbol. Let’s check out the result:

>>> **res**

9

You can also express one symbol in terms of another and substitute accordingly, using the subs() method. For example, if you knew that *x* = 1 - *y*, here’s how you could evaluate the preceding expression:

>>> **expr.subs({x:1-y})**

y**2 + 2*y*(-y + 1) + (-y + 1)**2

**PYTHON DICTIONARIES**

A dictionary is another type of data structure in Python (lists and tuples are other examples of data structures, which you’ve seen earlier). Dictionaries contain key-value pairs inside curly braces, where each key is matched up with a value, separated by a colon. In the preceding code listing, we entered the dictionary {x:1, y:2} as an argument to the subs() method. This dictionary has two key-value pairs—x:1 and y:2, where x and y are the keys and 1 and 2 are the corresponding values. You can retrieve a value from a dictionary by entering its associated key in brackets, much as we would retrieve an element from a list using its index. For example, here we create a simple dictionary and then retrieve the value corresponding to key1:

>>> **sampledict = {"key1": 5, "key2": 20}**

>>> **sampledict["key1"]**

5

To learn more about dictionaries, see *Appendix B*.

If you want the result to be simplified further—for example, if there are terms that cancel each other out, we can use SymPy’s simplify() function, as follows:

➊ >>> **expr_subs = expr.subs({x:1-y})**

>>> **from sympy import simplify**

➋ >>> **simplify(expr_subs)**

1

At ➊, we create a new label, expr_subs, to refer to the result of substituting *x* = 1 - *y* in the expression. We then import the simplify() function from SymPy and call it at ➋. The result turns out to be 1 because the other terms of the expression cancel each other.

Although there was a simplified version of the expression in the preceding example, you had to ask SymPy to simplify it using the simplify() function. Once again, this is because SymPy won’t do any simplification without being asked to.

The simplify() function can also simplify complicated expressions, such as those including logarithms and trigonometric functions, but we won’t get into that here.

**Calculating the Value of a Series**

Let’s revisit the series-printing program. In addition to printing the series, we want our program to be able to find the value of the series for a particular value of *x*. That is, our program will now take two inputs from the user—the number of terms in the series and the value of *x* for which the value of the series will be calculated. Then, the program will output both the series and the sum. The following program extends the series printing program to include these enhancements:

'''

Print the series:

x + x**2 + x**3 + ... + x**n

____ _____ _____

2 3 n

'''

from sympy import Symbol, pprint, init_printing

def print_series(n, x_value):

# Initialize printing system with reverse order

init_printing(order='rev-lex')

x = Symbol('x')

series = x

for i in range(2, n+1):

series = series + (x**i)/i

pprint(series)

# Evaluate the series at x_value

➊ series_value = series.subs({x:x_value})

print('Value of the series at {0}: {1}'.format(x_value, series_value))

if __name__ == '__main__':

n = input('Enter the number of terms you want in the series: ')

➋ x_value = input('Enter the value of x at which you want to evaluate the series: ')

print_series(int(n), float(x_value))

The print_series() function now takes an additional argument, x_value, which is the value of x for which the series should be evaluated. At ➊, we use the subs() method to perform the evaluation and the label series_value to refer to the result. In the next line, we display the result.

The additional input statement at ➋ asks the user to enter the value of x using the label x_value to refer to it. Before we call the print_series() function, we convert this value into its floating point equivalent using the float() function.

If you execute the program now, it will ask you for the two inputs and print out the series and the series value:

Enter the number of terms you want in the series: **5**

Enter the value of x at which you want to evaluate the series: **1.2**

x^{2} x^{3} x^{4} x^{5}

x + -- + -- + -- + --

2 3 4 5

Value of the series at 1.2: 3.51206400000000

In this sample run, we ask for five terms in the series, with x set to 1.2, and the program prints and evaluates the series.

*Converting Strings to Mathematical Expressions*

So far, we’ve been writing out individual expressions each time we want to do something with them. However, what if you wanted to write a more general program that could manipulate any expression provided by the user? For that, we need a way to convert a user’s input, which is a string, into something we can perform mathematical operations on. SymPy’s sympify() function helps us do exactly that. The function is so called because it converts the string into a SymPy object that makes it possible to apply SymPy’s functions to the input. Let’s see an example:

➊ >>> **from sympy import sympify**

>>> **expr = input('Enter a mathematical expression: ')**

Enter a mathematical expression: **x**2 + 3*x + x**3 + 2*x**

➋ >>> **expr = sympify(expr)**

We first import the sympify() function at ➊. We then use the input() function to ask for a mathematical expression as input, using the label expr to refer to it. Next, we call the sympify() function with expr as its argument at ➋ and use the same label to refer to the converted expression.

You can perform various operations on this expression. For example, let’s try multiplying the expression by 2:

>>> **2*expr**

2*x**3 + 2*x**2 + 10*x

What happens when the user supplies an invalid expression? Let’s see:

>>> **expr = input('Enter a mathematical expression: ')**

Enter a mathematical expression: **x**2 + 3*x + x**3 + 2x**

>>> **expr = sympify(expr)**

Traceback (most recent call last):

File "<pyshell#146>", line 1, in <module>

expr = sympify(expr)

File "/usr/lib/python3.3/site-packages/sympy/core/sympify.py", line 180, in sympify

raise SympifyError('could not parse %r' % a)

sympy.core.sympify.SympifyError: SympifyError: "could not parse 'x**2 + 3*x + x**3 + 2x'"

The last line tells us that sympify() isn’t able to convert the supplied input expression. Because this user didn’t add an operator between 2 and x, SymPy doesn’t understand what it means. Your program should expect such invalid input and print an error message if it comes up. Let’s see how we can do that by catching the SympifyError exception:

>>> **from sympy import sympify**

>>> **from sympy.core.sympify import SympifyError**

>>> **expr = input('Enter a mathematical expression: ')**

Enter a mathematical expression: **x**2 + 3*x + x**3 + 2x**

>>> **try:**

**expr = sympify(expr)****except SympifyError:**

**print('Invalid input')**

Invalid input

The two changes in the preceding program are that we import the SympifyError exception class from the sympy.core.sympify module and call the sympify() function in a try...except block. Now if there’s a SympifyError exception, an error message is printed.

**Expression Multiplier**

Let’s apply the sympify() function to write a program that calculates the product of two expressions:

'''

Product of two expressions

'''

from sympy import expand, sympify

from sympy.core.sympify import SympifyError

def product(expr1, expr2):

prod = expand(expr1*expr2)

print(prod)

if __name__=='__main__':

➊ expr1 = input('Enter the first expression: ')

➋ expr2 = input('Enter the second expression: ')

try:

expr1 = sympify(expr1)

expr2 = sympify(expr2)

except SympifyError:

print('Invalid input')

else:

➌ product(expr1, expr2)

At ➊ and ➋, we ask the user to enter the two expressions. Then, we convert them into a form understood by SymPy using the sympify() function in a try...except block. If the conversion succeeds (indicated by the else block), we call the product() function at ➌. In this function, we calculate the product of the two expressions and print it. Note how we use the expand() function to print the product so that all its terms are expressed as a sum of its constituent terms.

Here’s a sample execution of the program:

Enter the first expression: **x**2 + x*2 + x**

Enter the second expression: **x**3 + x*3 + x**

x**5 + 3*x**4 + 4*x**3 + 12*x**2

The last line displays the product of the two expressions. The input can also have more than one symbol in any of the expressions:

Enter the first expression: **x*y+x**

Enter the second expression: **x*x+y**

x**3*y + x**3 + x*y**2 + x*y

**Solving Equations**

SymPy’s solve() function can be used to find solutions to equations. When you input an expression with a symbol representing a variable, such as *x*, solve() calculates the value of that symbol. This function always makes its calculation by assuming the expression you enter is equal to zero—that is, it prints the value that, when substituted for the symbol, makes the entire expression equal zero. Let’s start with the simple equation *x* - 5 = 7. If we want to use solve() to find the value of x, we first have to make one side of the equation equal zero (*x* - 5 - 7 = 0). Then, we’re ready to use solve(), as follows:

>>> **from sympy import Symbol, solve**

>>> **x = Symbol('x')**

>>> **expr = x - 5 - 7**

>>> **solve(expr)**

[12]

When we use solve(), it calculates the value of 'x' as 12 because that’s the value that makes the expression (*x* - 5 - 7) equal to zero.

Note that the result 12 is returned in a list. An equation can have multiple solutions—for example, a quadratic equation has two solutions. In that case, the list will have all the solutions as its members. You can also ask the solve() function to return the result so that each member is dictionary instead. Each dictionary is composed of the symbol (variable name) and its value (the solution). This is especially useful when solving simultaneous equations where we have more than one variable to solve for because when the solution is returned as a dictionary, we know which solution corresponds to which variable.

*Solving Quadratic Equations*

In *Chapter 1*, we found the roots of the quadratic equation *ax*^{2} + *bx* + *c* = 0 by writing the formulas for the two roots and then substituting the values of the constants *a*, *b*, and *c*. Now, we’ll learn how we can use SymPy’s solve() function to find the roots without needing to write out the formulas. Let’s see an example:

➊ >>> **from sympy import solve**

>>> **x = Symbol('x')**

➋ >>> **expr = x**2 + 5*x + 4**

➌ >>> **solve(expr, dict=True)**

➍ [{x: -4}, {x: -1}]

The solve() function is first imported at ➊. We then define a symbol, x, and an expression corresponding to the quadratic equation, x**2 + 5*x + 4, at ➋. Then, we call the solve() function with the preceding expression at ➌. The second argument to the solve() function (dict=True) specifies that we want the result to be returned as a list of Python dictionaries.

Each solution in the returned list is a dictionary using the symbol as a key matched with its corresponding value. If the solution is empty, an empty list will be returned. The roots of the preceding equation are -4 and -1, as you can see at ➍.

We found out in the first chapter that the roots of the equation

*x*^{2} + *x* + 1 = 0

are complex numbers. Let’s attempt to find those using solve():

>>> **x=Symbol('x')**

>>> **expr = x**2 + x + 1**

>>> **solve(expr, dict=True)**

[{x: -1/2 - sqrt(3)*I/2}, {x: -1/2 + sqrt(3)*I/2}]

Both the roots are imaginary, as expected with the imaginary component indicated by the I symbol.

*Solving for One Variable in Terms of Others*

In addition to finding the roots of equations, we can take advantage of symbolic math to use the solve() function to express one variable in an equation in terms of the others. Let’s take a look at finding the roots for the generic quadratic equation *ax*^{2} + *bx* + *c* = 0. To do so, we’ll define *x* and three additional symbols—*a*, *b*, and *c*, which correspond to the three constants:

>>> **x = Symbol('x')**

>>> **a = Symbol('a')**

>>> **b = Symbol('b')**

>>> **c = Symbol('c')**

Next, we write the expression corresponding to the equation and use the solve() function on it:

>>> **expr = a*x*x + b*x + c**

>>> **solve(expr, x, dict=True)**

[{x: (-b + sqrt(-4*a*c + b**2))/(2*a)}, {x: -(b + sqrt(-4*a*c + b**2))/(2*a)}]

Here, we have to include an additional argument, x, to the solve() function. Because there’s more than one symbol in the equation, we need to tell solve() which symbol it should solve for, which is what we indicate by passing in x as the second argument. As we’d expect, solve() prints the quadratic formula: the generic formula for finding the value(s) of *x* in a polynomial expression.

To be clear, when we use solve() on an equation with more than one symbol, we specify the symbol to solve for as the second argument (and now the third argument specifies how we want the results to be returned).

Next, let’s consider an example from physics. According to one of the equations of motion, the distance traveled by a body moving with a constant acceleration *a*, with an initial velocity *u*, in time *t*, is given by

Given *u* and *a*, however, if you wanted to find the time required to travel a given distance, *s*, you’d have to first express *t* in terms of the other variables. Here’s how you could do that using SymPy’s solve() function:

>>> **from sympy import Symbol, solve, pprint**

>>> **s = Symbol('s')**

>>> **u = Symbol('u')**

>>> **t = Symbol('t')**

>>> **a = Symbol('a')**

>>> **expr = u*t + (1/2)*a*t*t - s**

>>> **t_expr = solve(expr,t, dict=True)**

>>> **pprint(t_expr)**

The result looks like this:

Now that we have the expression for *t* (referred to by the label t_expr), we can use the subs() method to replace the values of *s*, *u*, and *a* to find the two possible values of *t*.

*Solving a System of Linear Equations*

Consider the following two equations:

2*x* + 3*y* = 6

3*x* + 2*y* = 12

Say we want to find the pair of values (*x*, *y*) that satisfies both the equations. We can use the solve() function to find the solution for a system of equations like this one.

First, we define the two symbols and create the two equations:

>>> **x = Symbol('x')**

>>> **y = Symbol('y')**

>>> **expr1 = 2*x + 3*y - 6**

>>> **expr2 = 3*x + 2*y - 12**

The two equations are defined by the expressions expr1 and expr2, respectively. Note how we’ve rearranged the expressions so they both equal zero (we moved the right side of the given equations to the left side). To find the solution, we call the solve() function with the two expressions forming a tuple:

>>> **solve((expr1, expr2), dict=True)**

[{y: -6/5, x: 24/5}]

As I mentioned earlier, getting the solution back as a dictionary is useful here. We can see that the value of x is 24/5 and the value of y is -6/5. Let’s verify whether the solution we got really satisfies the equations. To do so, we’ll first create a label, soln, to refer to the solution we got and then use the subs() method to substitute the corresponding values of x and y in the two expressions:

>>> **soln = solve((expr1, expr2), dict=True)**

>>> **soln = soln[0]**

>>> **expr1.subs({x:soln[x], y:soln[y]})**

0

>>> **expr2.subs({x:soln[x], y:soln[y]})**

0

The result of substituting the values of x and y corresponding to the solution in the two expressions is zero.

**Plotting Using SymPy**

In *Chapter 2*, we learned to make graphs where we explicitly specified the numbers we wanted to plot. For example, to plot the graph of the gravitational force against the distance between two bodies, you had to calculate the gravitational force for each distance value and supply the lists of distances and forces to matplotlib. With SymPy, on the other hand, you can just tell SymPy the equation of the line you want to plot, and the graph will be created for you. Let’s plot a line whose equation is given by *y* = 2*x* + 3:

>>> **from sympy.plotting import plot**

>>> **from sympy import Symbol**

>>> **x = Symbol('x')**

>>> **plot(2*x+3)**

All we had to do was import plot and Symbol from sympy.plotting, create a symbol, x, and call the plot() function with the expression 2*x+3. SymPy takes care of everything else and plots the graph of the function, as shown in *Figure 4-1*.

*Figure 4-1: Plot of the line* y = *2*x + *3*

The graph shows that a default range of *x* values was automatically chosen: -10 to 10. You may notice that the graph window looks very similar to those you saw in *Chapters 2* and *3*. That’s because SymPy uses matplotlib behind the scenes to draw the graphs. Also note that we didn’t have to call the show() function to show the graphs because this is done automatically by SymPy.

Now, let’s say that you wanted to limit the values of 'x' in the preceding graph to lie in the range -5 to 5 (instead of -10 to 10). You’d do that as follows:

>>> **plot((2*x + 3), (x, -5, 5))**

Here, a tuple consisting of the symbol, the lower bound, and the upper bound of the range—(x, -5, 5)—is specified as the second argument to the plot() function. Now, the graph displays only the values of *y* corresponding to the values of *x* between -5 and 5 (see *Figure 4-2*).

*Figure 4-2: Plot of the line* y = *2*x + *3 with the values of* x *restricted to the range -5 to 5*

You can use other keyword arguments in the plot() function, such as title to enter a title or xlabel and ylabel to label the *x*-axis and the *y*-axis, respectively. The following plot() function specifies the preceding three keyword arguments (see the corresponding graph in *Figure 4-3*):

>>> **plot(2*x + 3, (x, -5, 5), title='A Line', xlabel='x', ylabel='2x+3')**

*Figure 4-3: Plot of the line* y = *2*x + *3 with the range of* x *and other attributes specified*

The plot shown in *Figure 4-3* now has a title and labels on the *x*-axis and the *y*-axis. You can specify a number of other keyword arguments to the plot() function to customize the behavior of the function as well as the graph itself. The show keyword argument allows us to specify whether we want the graph to be displayed. Passing show=False will cause the graph to not be displayed when you call the plot() function:

>>> **p = plot(2*x + 3, (x, -5, 5), title='A Line', xlabel='x', ylabel='2x+3', show=False)**

You will see that no graph is shown. The label p refers to the plot that is created, so you can now call p.show() to display the graph. You can also save the graph as an image file using the save() method, as follows:

>>> **p.save('line.png')**

This will save the plot to a file *line.png* in the current directory.

*Plotting Expressions Input by the User*

The expression that you pass to the plot() function must be expressed in terms of *x* only. For example, earlier we plotted *y* = 2*x* + 3, which we entered to the plot function as simply 2*x* + 3. If the expression were not originally in this form, we’d have to rewrite it. Of course, we could do this manually, outside the program. But what if you want to write a program that allows its users to graph any expression? If the user enters an expression in the form of 2*x* + 3*y* - 6, say, we have to first convert it. The solve() function will help us here. Let’s see an example:

>>> **expr = input('Enter an expression: ')**

Enter an expression: **2*x + 3*y - 6**

➊ >>> **expr = sympify(expr)**

➋ >>> **y = Symbol('y')**

>>> **solve(expr, y)**

➌ [-2*x/3 + 2]

At ➊, we use the sympify() function to convert the input expression to a SymPy object. At ➋, we create a Symbol object to represent 'y' so that we can tell SymPy which variable we want to solve the equation for. Then we solve the expression to find y in terms of x by specifying y as the second argument to the solve() function. At ➌, this returns the equation in terms of x, which is what we need for plotting.

Notice that this final expression is stored in a list, so before we can use it, we’ll have to extract it from the list:

>>> **solutions = solve(expr, 'y')**

➍ >>> **expr_y = solutions[0]**

>>> **expr_y**

-2*x/3 + 2

We create a label, solutions, to refer to the result returned by the solve() function, which is a list with only one item. Then, we extract that item at ➍. Now, we can call the plot() function to graph the expression. The next listing shows a full graph-drawing program:

'''

Plot the graph of an input expression

'''

from sympy import Symbol, sympify, solve

from sympy.plotting import plot

def plot_expression(expr):

y = Symbol('y')

solutions = solve(expr, y)

expr_y = solutions[0]

plot(expr_y)

if __name__=='__main__':

expr = input('Enter your expression in terms of x and y: ')

try:

expr = sympify(expr)

except SympifyError:

print('Invalid input')

else:

plot_expression(expr)

Note that the preceding program includes a try...except block to check for invalid input, as we’ve done with sympify() earlier. When you run the program, it asks you to input an expression, and it will create the corresponding graph.

*Plotting Multiple Functions*

You can enter multiple expressions when calling the SymPy plot function to plot more than one expression on the same graph. For example, the following code plots two lines at once (see *Figure 4-4*):

>>> **from sympy.plotting import plot**

>>> **from sympy import Symbol**

>>> **x = Symbol('x')**

>>> **plot(2*x+3, 3*x+1)**

*Figure 4-4: Plotting two lines on the same graph*

This example brings out another difference between plotting in matplotlib and in SymPy. Here, using SymPy, both lines are the same color, whereas matplotlib would have automatically made the lines different colors. To set different colors for each line with SymPy, we’ll need to perform some extra steps, as shown in the following code, which also adds a legend to the graph:

>>> **from sympy.plotting import plot**

>>> **from sympy import Symbol**

>>> **x = Symbol('x')**

➊ >>> **p = plot(2*x+3, 3*x+1, legend=True, show=False)**

➋ >>> **p[0].line_color = 'b'**

➌ >>> **p[1].line_color = 'r'**

>>> **p.show()**

At ➊, we call the plot() function with the equations for the two lines but pass two additional keyword arguments—legend and show. By setting the legend argument to True, we add a legend to the graph, as we saw in *Chapter 2*. Note, however, that the text that appears in the legend will match the expressions you plotted—you can’t specify any other text. We also set show=False because we want to set the color of the lines before we draw the graph. The statement at ➋, p[0], refers to the first line, 2*x* + 3, and we set its attribute line_color to 'b', meaning that we want this line to be blue. Similarly, we set the color of the second plot to red using the string 'r' ➌. Finally, we call the show() to display the graph (see *Figure 4-5*).

*Figure 4-5: Plot of the two lines with each line drawn in a different color*

In addition to red and blue, you can plot the lines in green, cyan, magenta, yellow, black, and white (using the first letter of the color in each case).

**What You Learned**

In this chapter, you learned the basics of symbolic math using SymPy. You learned about declaring symbols, constructing expressions using symbols and mathematical operators, solving equations, and plotting graphs. You will be learning more features of SymPy in later chapters.

**Programming Challenges**

Here are a few programming challenges that should help you further apply what you’ve learned. You can find sample solutions at * http://www.nostarch.com/doingmathwithpython/*.

*#1: Factor Finder*

You learned about the factor() function, which prints the factors of an expression. Now that you know how your program can handle expressions input by a user, write a program that will ask the user to input an expression, calculate its factors, and print them. Your program should be able to handle invalid input by making use of exception handling.

*#2: Graphical Equation Solver*

Earlier, you learned how to write a program that prompts the user to input an expression such as 3*x* + 2*y* - 6 and create the corresponding graph. Write a program that asks the user for two expressions and then graphs them both, as follows:

>>> **expr1 = input('Enter your first expression in terms of x and y: ')**

>>> **expr2 = input('Enter your second expression in terms of x and y: ')**

Now, expr1 and expr2 will store the two expressions input by the user. You should convert both of these into SymPy objects using the sympify() step in a try...except block.

All you need to do from here is plot these two expressions instead of one.

Once you’ve completed this, enhance your program to print the solution—the pair of *x* and *y* values that satisfies both equations. This will also be the spot where the two lines on the graph intersect. (Hint: Refer to how we used the solve() function earlier to find the solution of a system of two linear equations.)

*#3: Summing a Series*

We saw how to find the sum of a series in “*Printing a Series*” on *page 99*. There, we manually added the terms of the series by looping over all the terms. Here’s a snippet from that program:

for i in range(2, n+1):

series = series + (x**i)/i

SymPy’s summation() function can be directly used to find such summations. The following example prints the sum of the first five terms of the series we considered earlier:

>>> **from sympy import Symbol, summation, pprint**

>>> **x = Symbol('x')**

>>> **n = Symbol('n')**

➊ >>> **s = summation(x**n/n, (n, 1, 5))**

>>> **pprint(s)**

x^{5} x^{4} x^{3} x^{2}

-- + -- + -- + -- + x

5 4 3 2

We call the summation() function at ➊, with the first argument being the *n*th term of the series and the second argument being a tuple that states the range of *n*. We want the sum of the first five terms here, so the second argument is (n, 1, 5).

Once you have the sum, you can use the subs() method to substitute a value for *x* to find the numerical value of the sum:

>>> **s.subs({x:1.2})**

3.51206400000000

Your challenge is to write a program that’s capable of finding the sum of an arbitrary series when you supply the *n*th term of the series and the number of terms in it. Here’s an example of how the program would work:

Enter the nth term: **a+(n-1)*d**

Enter the number of terms: **3**

3·a + 3·d

In this example, the *n*th term supplied is that of an *arithmetic progression*. Starting with a and d as the *common difference*, the number of terms up to which the sum is to be calculated is 3. The sum turns out to be 3a + 3d, which agrees with the known formula for the same.

*#4: Solving Single-Variable Inequalities*

You’ve seen how to solve an equation using SymPy’s solve() function. But SymPy is also capable of solving single-variable inequalities, such as *x* + 5 > 3 and sin*x* - 0.6 > 0. That is, SymPy can solve relations besides equality, like >, <, and so on. For this challenge, create a function,isolve(), that will take any inequality, solve it, and then return the solution.

First, let’s learn about the SymPy functions that will help you implement this. The inequality-solving functions are available as three separate functions for polynomial, rational, and all other inequalities. We’ll need to pick the right function to solve various inequalities, or we’ll get an error.

A *polynomial* is an algebraic expression consisting of a variable and coefficients and involving only the operations of addition, subtraction, and multiplication and only positive powers of the variable. An example of a polynomial inequality is *x*^{2} + 4 < 0.

To solve a polynomial inequality, use the solve_poly_inequality() function:

>>> **from sympy import Poly, Symbol, solve_poly_inequality**

>>> **x = Symbol('x')**

➊ >>> **ineq_obj = -x**2 + 4 < 0**

➋ >>> **lhs = ineq_obj.lhs**

➌ >>> **p = Poly(lhs, x)**

➍ >>> **rel = ineq_obj.rel_op**

>>> **solve_poly_inequality(p, rel)**

[(-oo, -2), (2, oo)]

First, create the expression representing an inequality, -*x*^{2} + 4 < 0, at ➊ and refer to this expression with the label ineq_obj. Then, extract the left side of the inequality—that is, the algebraic expression -*x*^{2} + 4—using the lhs attribute at ➋. Next, create a Poly object at ➌ to represent the polynomial we extracted at ➋. The second argument passed when creating the object is the symbol object that represents the variable, x. At ➍, extract the relational operator from the inequality object using the rel attribute. Finally, call the solve_poly_inequality() function with the polynomial object, p, and rel as the two arguments. The program returns the solution as a list of tuples, with each tuple representing a solution for the inequality as the lower limit and the upper limit of the range of numbers. For this inequality, the solution is all numbers less than -2 and all numbers greater than 2.

A *rational expression* is an algebraic expression in which the numerator and denominator are both polynomials. Here’s an example of a rational inequality:

For rational inequalities, use the solve_rational_inequalities() function:

>>> **from sympy import Symbol, Poly, solve_rational_inequalities**

>>> **x = Symbol('x')**

➊ >>> **ineq_obj = ((x-1)/(x+2)) > 0**

>>> **lhs = ineq_obj.lhs**

➋ >>> **numer, denom = lhs.as_numer_denom()**

>>> **p1 = Poly(numer)**

>>> **p2 = Poly(denom)**

>>> **rel = ineq_obj.rel_op**

➌ >>> **solve_rational_inequalities([[((p1, p2), rel)]])**

(-oo, -2) U (1, oo)

Create an inequality object representing our example rational inequality at ➊ and then extract the rational expression using the lhs attribute. Separate out the numerator and the denominator into the labels numer and denom using the as_numer_denom() method at ➋, which returns a tuple with the numerator and denominator as the two members. Then, create two polynomial objects, p1 and p2, representing the numerator and denominator, respectively. Retrieve the relational operator and call the solve_rational_inequalities() function, passing it the two polynomial objects—p1 andp2—and the relational operator.

The program returns the solution (-oo, -2) U (1, oo), where U denotes that the solution is a *union* of the two *sets* of solutions consisting of all numbers less than -2 and all numbers greater than 1. (We’ll learn about sets in *Chapter 5*.)

Finally, sin*x* - 0.6 > 0 is an example of an inequality that belongs to neither the polynomial nor rational expression categories. If you have such an inequality to solve, use the solve_univariate_inequality() function:

>>> **from sympy import Symbol, solve, solve_univariate_inequality, sin**

>>> **x = Symbol('x')**

>>> **ineq_obj = sin(x) - 0.6 > 0**

>>> **solve_univariate_inequality(ineq_obj, x, relational=False)**

(0.643501108793284, 2.49809154479651)

Create an inequality object representing the inequality sin(x) - 0.6 > 0 and then call the solve_univariate_inequality() function with the first two arguments as the inequality object, ineq_obj, and the symbol object, x. The keyword argument relational=False specifies to the function that we want the solution to be returned as a *set*. The solution for this inequality turns out to be all numbers lying between the first and second members of the tuple the program returns.

**Hints: Handy Functions**

Now remember—your challenge is (1) to create a function, isolve(), that will take any inequality and (2) to choose one of the appropriate functions discussed in this section to solve it and return the solution. The following hints may be useful to implement this function.

The is_polynomial() method can be used to check whether an expression is a polynomial or not:

>>> **x = Symbol('x')**

>>> **expr = x**2 - 4**

>>> **expr.is_polynomial()**

True

>>> **expr = 2*sin(x) + 3**

>>> **expr.is_polynomial()**

False

The is_rational_function() can be used to check whether an expression is a rational expression:

>>> **expr = (2+x)/(3+x)**

>>> **expr.is_rational_function()**

True

>>> **expr = 2+x**

>>> **expr.is_rational_function()**

True

>>> **expr = 2+sin(x)**

>>> **expr.is_rational_function()**

False

The sympify() function can convert an inequality expressed as a string to an inequality object:

>>> **from sympy import sympify**

>>> **sympify('x+3>0')**

x + 3 > 0

When you run your program, it should ask the user to input an inequality expression and print back the solution.