Instant SymPy Starter

May 2013


Learn to use SymPy's symbolic engine to simplify Python calculations

(For more resources related to this topic, see here.)


To compute the derivative of a function, create the corresponding expression and use diff(). Its first argument is the expression and the second is the variable with regard to which you want to differentiate. The result is the expression for the derivative:

>>> diff(exp(x**2), x)
>>> diff(x**2 * y**2, y)

Higher-order derivatives can also be computed with a single call to diff():

>>> diff(x**3, x, x)
>>> diff(x**3, x, 2)
>>> diff(x**2 * y**2, x, 2, y, 2)

Due to SymPy's focus on expressions rather than functions, the derivatives for symbolic functions can seem a little surprising, but LaTeX rendering in the notebook should make their meaning clear.

>>> f = Function('f')
>>> diff(f(x**2), x)
2*x*Subs(Derivative(f(_xi_1), _xi_1), (_xi_1,), (x**2,))

Let's take a look at the following screenshot:


Limits are obtained through limit(). The syntax for the limit of expr when x goes to some value x0 is limit(expr, x, x0). To specify a limit towards infinity, you need to use SymPy's infinity object, named oo. This object will also be returned for infinite limits:

>>> limit(exp(-x), x, oo)
>>> limit(1/x**2, x, 0)

There is also a fourth optional parameter, to specify the direction of approach of the limit target. "+" (the default) gives the limit from above, and "-" is from below. Obviously, this parameter is ignored when the limit target is infinite:

>>> limit(1/x, x, 0, "-")
>>> limit(1/x, x, 0, "+")

Let's take a look at the following screenshot:


SymPy has powerful algorithms for integration, and, in particular, can find most integrals of logarithmic and exponential functions expressible with special functions, and many more besides, thanks to Meijer G-functions.

The main function for integration is integrate(). It can compute both antiderivatives (indefinite integrals) and definite integrals. Note that the value of an antiderivative is only defined up to an arbitrary constant but the result does not include it.

>>> integrate(sin(x), x)
>>> integrate(sin(x), (x, 0, pi))

Unevaluated symbolic integrals and antiderivatives are represented by the Integral class. integrate() may return these objects if it cannot compute the integral.

It is also possible to create Integral objects directly, using the same syntax as integrate(). To evaluate them, call their .doit() method:

>>> integral = Integral(sin(x), (x, 0, pi))
>>> integral
Integral(sin(x), (x, 0, pi))
>>> integral.doit()

Let's take a look at the following screenshot:

Taylor series

A Taylor series approximation is an approximation of a function obtained by truncating its Taylor series. To compute it, use series(expr, x, x0, n), where x is the relevant variable, x0 is the point where the expansion is done (defaults to 0), and n is the order of expansion (defaults to 6):

>>> series(cos(x), x)
1 - x**2/2 + x**4/24 + O(x**6)
>>> series(cos(x), x, n=10)
1 - x**2/2 + x**4/24 - x**6/720 + x**8/40320 + O(x**10)

The O(x**6) part in the result is a "big-O" object. Intuitively, it represents all the terms of order equal to or higher than 6. This object automatically absorbs or combines with powers of the variable, which makes simple arithmetic operations on expansions convenient:

>>> O(x**2) + 2*x**3
>>> O(x**2) * 2*x**3
>>> expand(series(sin(x), x, n=6) * series(cos(x), x, n=4))
x - 2*x**3/3 + O(x**5)
>>> series(sin(x)*cos(x), x, n=5)
x - 2*x**3/3 + O(x**5)

If you want to use the expansion as an approximation of the function, the O() term prevents it from behaving like an ordinary expression, so you need to remove it. You can do so by using the aptly named .removeO() method:

>>> series(cos(x), x).removeO()
x**4/24 - x**2/2 + 1

Taylor series look better in the notebook, as shown in the following screenshot:

Solving equations

This section will teach you how to solve the different types of equations that SymPy handles.

The main function to use for solving equations is solve(). Its interface is somewhat complicated as it accepts many different kinds of inputs and can output results in various forms depending on the input.

In the simplest case, univariate equations, use the syntax solve(expr, x) to solve the equation expr = 0 for the variable x. If you want to solve an equation of the form A = B, simply put it under the preceding form, using solve(A - B, x). This can solve algebraic and transcendental equations involving rational fractions, square roots, absolute values, exponentials, logarithms, trigonometric functions, and so on. The result is then a list of the values of the variables satisfying the equation.

The following commands show a few examples of equations that can be solved:

>>> solve(x**2 - 1, x)
[-1, 1]
>>> solve(x*exp(x) - 1, x)
>>> solve(abs(x**2-4) - 3, x)
[-1, 1, -sqrt(7), sqrt(7)]

Note that the form of the result means that it can only return a finite set of solutions. In cases where the true solution is infinite, it can therefore be misleading. When the solution is an interval, solve() typically returns an empty list. For periodic functions, usually only one solution is returned:

>>> solve(0, x) # all x are solutions
>>> solve(x - abs(x), x) # all positive x are solutions
>>> solve(sin(x), x) # all k*pi with k integer are solutions

The domain over which the equation is solved depends on the assumptions on the variable. Hence, if the variable is a real Symbol object, only real solutions are returned, but if it is complex, then all solutions in the complex plane are returned (subject to the aforementioned restriction on returning infinite solution sets). This difference is readily apparent when solving polynomials, as the following example demonstrates:

>>> solve(x**2 + 1, x)
>>> solve(z**2 + 1, z)
[-I, I]

There is no restriction on the number of variables appearing in the expression. Solving a multivariate expression for any of its variables allows it to be expressed as a function of the other variables, and to eliminate it from other expressions. The following example shows different ways of solving the same multivariate expression:

>>> solve(x**2 - exp(a), x)
[-exp(a/2), exp(a/2)]
>>> solve(x**2 - exp(a), a)
>>> solve(x**2 - exp(a), x, a)
[{x: -exp(a/2)}, {x: exp(a/2)}]
>>> solve(x**2 - exp(a), x, b)
[{x: -exp(a/2)}, {x: exp(a/2)}]

To solve a system of equations, pass a list of expressions to solve(): each one will be interpreted, as in the univariate case, as an equation of the form expr = 0. The result can be returned in one of two forms, depending on the mathematical structure of the input: either as a list of tuples, where each tuple contains the values for the variables in the order given to solve, or a single dictionary, suitable for use in subs(), mapping variables to their values.

As you can see in the following example, it can be hard to predict what form the result will take:

>>> solve([exp(x**2) - y, y - 3], x, y)
[(-sqrt(log(3)), 3), (sqrt(log(3)), 3)]
>>> solve([x**2 - y, y - 3], x, y)
[(-sqrt(3), 3), (sqrt(3), 3)]
>>> solve([x - y, y - 3], x, y)
{y: 3, x: 3}

This variability in return types is fine for interactive use, but for library code, more predictability is required. In this case, you should use the dict=True option. The output will then always be a list of mappings of variables to value. Compare the following example to the previous one:

>>> solve([x**2 - y, y - 3], x, y, dict=True)
[{y: 3, x: -sqrt(3)}, {y: 3, x: sqrt(3)}]
>>> solve([x - y, y - 3], x, y, dict=True)
[{y: 3, x: 3}]


We successfully computed the various mathematical operations using the SymPy application, Calculus.

Resources for Article :

Further resources on this subject:

Books to Consider

comments powered by Disqus

An Introduction to 3D Printing

Explore the future of manufacturing and design  - read our guide to 3d printing for free