Symbolic Computations - SymPy

In this chapter, we will give a brief introduction to using Python for symbolic computations. There is powerful software on the market for performing symbolic computations, for example, Maple™ or Mathematica™. But sometimes, it might be favorable to make symbolic calculations in the language or framework you are used to. At this stage of the book, we assume that this language is Python, so we seek a tool in Python—the module SymPy.

A complete description of SymPy, if even possible, would fill an entire book, and that is not the purpose of this chapter. Instead, we will stake out a path into this tool by examining some guiding examples, giving a flavor of the potential of this tool as a complement to NumPy and SciPy.

16.1 What are symbolic computations?

All computations we did so far in this book were so-called numeric computations. These were a sequence of operations mainly on floating-point numbers. It is the nature of numeric computations that the result is an approximation of the exact solution.

Symbolic computations operate on formulas or symbols by transforming them as taught in algebra or calculus into other formulas. The last step of these transformations might then require that numbers are inserted and a numeric evaluation is performed.

We illustrate the difference by computing this definite integral:

Symbolically this expression can be transformed by considering the primitive function of the integrand:

We now obtain a formula for the definite integral by inserting the integral bounds:

This is called a closed-form expression for the integral. Very few mathematical problems have a solution that can be given in a closed-form expression. It is the exact value...

16.1.1 Elaborating an example in SymPy

To begin with, let's elaborate on the previous example in SymPy and explain the steps.

First, we have to import the module:

from sympy import *
init_printing()

The second command makes sure that formulas are presented in a graphical way, if possible. Then, we generate a symbol and define the integrand:

x = symbols('x')
f = Lambda(x, 1/(x**2 + x + 1))

x is now a Python object of type Symbol and f is a SymPy Lambda function (note the command starting with a capital letter).

Now we start with the symbolic computation of the integral:

integrate(f(x),x)    

Depending on your working environment, the result is presented in different ways; see the following screenshot (Figure 16.2), which represents two different results of a SymPy formula in different environments:

Figure 16.2: Two screenshots of a SymPy presentation of formula in two different environments

We can check by differentiation...

16.2 Basic elements of SymPy

Here we introduce the basic elements of SymPy. You will find it favorable to be already familiar with classes and data types in Python.

16.2.1 Symbols – the basis of all formulas

The basic construction element to build a formula in SymPy is the symbol. As we saw in the introductory example, a symbol is created by the command symbols. This SymPy command generates symbol objects from a given string:

x, y, mass, torque = symbols('x y mass torque')

It is actually a short form of the following command:

symbol_list=[symbols(l) for l in 'x y mass torque'.split()]

Followed by an unpacking step to obtain variables:

 x, y, mass, torque = symbol_list

The arguments of the command define the string representation of the symbol. The variable name of the symbol chosen is often identical to its string representation, but this is not required by the language:

row_index=symbols('i',integer=True)
print(row_index**2)  # returns i**2

Here, we also defined that the symbol is assumed to be an integer.

An entire set of symbols can be defined in a very compact way:

integervariables = symbols...

16.2.2 Numbers

Python evaluates operations on numbers directly and introduces unavoidable rounding errors. These would obstruct all symbolic calculations. This is avoided when we sympify numbers:

1/3  # returns 0.3333333333333333
sympify(1)/sympify(3)  # returns '1/3'

The sympify command converts an integer to an object of type sympy.core.numbers.Integer.

Instead of writing 1/3 as an operation of two integers, it can also be represented directly as a rational number by Rational(1,3).

16.2.3 Functions

SymPy distinguishes between defined and undefined functions. The term undefined functions (which might be a bit misleading) refers to well-defined Python objects for generic functions that have no special properties.

An example of a function with special properties is atan or the Lambda function used in the introductory example of this chapter.  

Note the different names for the different implementations of the same mathematical function: sympy.atan and scipy.arctan.

Undefined functions

A symbol for an undefined function is created by giving the symbols command an extra class argument:

f, g = symbols('f g', cls=Function)

The same can be achieved by using the constructor Function :

f = Function('f')
g = Function('g')

With undefined functions, we can evaluate the general rules of calculus.

For example, let's evaluate the following expression:

This is symbolically computed in Python by using the following command:

x = symbols('x')
f, g = symbols('f g', cls=Function)
diff(f(x*g(x)),x)

When executed, the previous code returns the following as output:

This example shows how the product rule and the chain rule were applied.

We can even use an undefined function as a function in several variables, for example:

x = symbols('x:3')
f(*x)

which returns the following output:

Note the use of the star operator to unpack a tuple to form f with arguments;&...

16.2.4 Elementary functions

Examples for elementary functions in SymPy are trigonometric functions and their inverses. The following example shows how simplify acts on an expression which includes elementary functions:

x = symbols('x')
simplify(cos(x)**2 + sin(x)**2)  # returns 1

Here is another example of the use of elementary functions:

atan(x).diff(x) - 1./(x**2+1)  # returns 0

If you use SciPy and SymPy together, we strongly recommend that you use them in different namespaces:

import numpy as np
import sympy as sym
# working with numbers
x=3
y=np.sin(x)
# working with symbols
x=sym.symbols('x')
y=sym.sin(x)   

16.2.5 Lambda functions

In Section 7.7: Anonymous functions, we saw how to define so-called anonymous functions in Python. The SymPy counterpart is the command Lambda. Note the difference; lambda is a keyword while Lambda is a constructor.

The command Lambda takes two arguments, the symbol of the function's independent variable, and a SymPy expression to evaluate the function.

Here is an example that defines air resistance (also called drag) as a function of speed:

C,rho,A,v=symbols('C rho A v')
# C drag coefficient, A coss-sectional area, rho density
# v speed
f_drag = Lambda(v,-Rational(1,2)*C*rho*A*v**2)

f_drag is displayed as a graphical expression:

This function can be evaluated in the usual way by providing it with an argument:

x = symbols('x')
f_drag(2)
f_drag(x/3)

Which results in the expressions:

It is also possible to create functions in several variables by just providing the first parameter...

16.3 Symbolic linear algebra

Symbolic linear algebra is supported by SymPy's data type matrix which we will introduce first. Then we will present some linear algebra methods as examples for the broad spectrum of possibilities for symbolic computations in this field.

16.3.1 Symbolic matrices

We briefly met the matrix data type when we discussed vector-valued functions. There, we saw it in its simplest form, which converts a list of lists into a matrix. To see an example, let's construct a rotation matrix:

phi=symbols('phi')
rotation=Matrix([[cos(phi), -sin(phi)],
                 [sin(phi), cos(phi)]])

When working with SymPy matrices we have to note that the operator * performs matrix multiplications and is not acting as an elementwise multiplication, which is the case for NumPy arrays. 

The previously defined rotation matrix can be checked for orthogonality by using this matrix multiplication and the transpose of a matrix:

simplify(rotation.T*rotation -eye(2))  # returns a 2 x 2 zero matrix

The previous example shows how a matrix is transposed and how the identity matrix is created. Alternatively, we could have checked whether its inverse is its transpose, which can be done as follows:

simplify...

16.3.2 Examples for linear algebra methods in SymPy

The basic task in linear algebra is to solve linear equation systems:

Let's do this symbolically for a  matrix:

A = Matrix(3,3,symbols('A1:4(1:4)'))
b = Matrix(3,1,symbols('b1:4'))
x = A.LUsolve(b)

The output of this relatively small problem is already merely readable, which can be seen in the following graphical expression:

               

Again, the use of simplify command helps us to detect canceling terms and to collect common factors:

simplify(x)

Which will result in the following output, which looks much better:

             

Symbolic computations become very slow when increasing matrix dimensions. For dimensions bigger than 15, memory problems might even occur.

The next figure (Figure 16.3) illustrates the differences in CPU time between symbolically and numerically solving a linear system...

16.4 Substitutions

Let's first consider a simple symbolic expression:

x, a = symbols('x a')
b = x + a

What happens if we set x = 0 ? We observe that b did not change. What we did was that we changed the Python variable x. It now no longer refers to the symbol object but to the integer object 0. The symbol represented by the string 'x' remains unaltered, and so does b.

Instead, altering an expression by replacing symbols with numbers, other symbols, or expressions is done by a special substitution method, which can be seen in the following code:

x, a = symbols('x a')
b = x + a
c = b.subs(x,0)   
d = c.subs(a,2*a)  
print(c, d)   # returns (a, 2a)

This method takes one or two arguments. The following two statements are equivalent:

b.subs(x,0)
b.subs({x:0})  # a dictionary as argument

Dictionaries as arguments allow us to make several substitutions in one step:

b.subs({x:0, a:2*a})  # several substitutions...

16. 5 Evaluating symbolic expressions

In the context of scientific computing, there is often the need to first make symbolic manipulations and then convert the symbolic result into a floating-point number.

The central tool for evaluating a symbolic expression is evalf. It converts symbolic expressions to floating-point numbers by using the following:

pi.evalf()   # returns 3.14159265358979

The data type of the resulting object is Float (note the capitalization), which is a SymPy data type that allows floating-point numbers with an arbitrary number of digits (arbitrary precision).

The default precision corresponds to 15 digits, but it can be changed by giving evalf an extra positive integer argument

specifying the desired precision in terms of the numbers of digits:

pi.evalf(30)   # returns  3.14159265358979323846264338328

A consequence of working with arbitrary precision is that numbers can be arbitrarily small, that is, the limits of the classical...

16.5.1 Example: A study on the convergence order of Newton's method

An iterative method that iterates  is said to converge with order  with _, if there exists a positive constant  such that

Newton's method, when started with a good initial value, has order , and for certain problems, even . Newton's method when applied to the problem  gives the following iteration scheme:

Which converges cubically; that is, q = 3.

This implies that the number of correct digits triples from iteration to iteration. To demonstrate cubic convergence and to numerically determine the constant  is hardly possible with the standard 16-digit float data type.

The following code uses SymPy together with high-precision evaluation instead and takes the study on cubic convergence to the extreme:

import sympy as sym
x = sym.Rational(1,2)
xns=[x]

for i in range(1,9):
    x = (x - sym.atan(x)*(1+x**2)).evalf(3000)
    xns.append(x)

The result...

16.5.2 Converting a symbolic expression into a numeric function

As we have seen, the numerical evaluation of symbolic expressions is done in three steps: first, we do some symbolic computations and then we substitute values by numbers and do an evaluation to a floating-point number with evalf.

The reason for symbolic computations is often that we want to make parameter studies. This requires that the parameter is modified within a given parameter range. This requires that a symbolic expression is eventually turned into a numeric function.

A study on the parameter dependency of polynomial coefficients

We demonstrate a symbolic/ numeric parameter study by an interpolation example to introduce the SymPy command lambdify.

Let's consider the task to interpolate the data  and . Here,  is a free parameter, which we will vary over the interval.

The quadratic interpolation polynomial has coefficients depending on this parameter:

Using SymPy and the monomial approach described in Exercise 3 in Section 4.11: Exercises gives us the closed formula for these coefficients:

t=symbols('t')
x=[0,t,1]
# The Vandermonde Matrix
V = Matrix([[0, 0, 1], [t**2, t, 1], [1, 1,1]])
y = Matrix([0,1,-1]) # the data vector
a = simplify(V.LUsolve(y)) # the coefficients
# the leading coefficient as a function of the parameter
a2 = Lambda(t,a[0])

We obtain a symbolic function for the leading coefficient  of the interpolation polynomial:

Now it is time to turn the expression into a numeric function...

16.6 Summary

In this chapter, you were introduced to the world of symbolic computations and you got a glimpse of the power of SymPy. By following the examples, you learned how to set up symbolic expressions, how to work with symbolic matrices, and you saw how to make simplifications. Working with symbolic functions and transforming them into numerical evaluations, finally, established the link to scientific computing and floating-point results. You experienced the strength of SymPy as you used its full integration into Python with its powerful constructs and legible syntax.

Consider this last chapter as an appetizer rather than a complete menu. We hope you became hungry for future fascinating programming challenges in scientific computing and mathematics.