Not only symbolics: numerical computing

Symbolic mathematics can’t exist without numerical methods. Most “symbolic” modules in SymPy take at least some advantage of numerical computing. SymPy uses the mpmath library for this purpose.

Let’s start from something simple and find numerical approximation to \(\pi\). Normally SymPy represents \(\pi\) as a symbolic entity:

>>> pi
>>> type(_)
<class 'sympy.core.numbers.Pi'>

To obtain numerical approximation of \(\pi\) we can use either the evalf() method or N(), which is a simple wrapper over the former method:

>>> pi.evalf()

The default precision is 15 digits. We can change this using the n parameter:

>>> pi.evalf(n=30)

The mpmath library implements arbitrary precision floating point arithmetics (limited only by available memory), so we can set n to a very big value, e.g. one million:

>>> million_digits = pi.evalf(n=1000000)
>>> str(million_digits)[-1]

evalf() can handle much more complex expressions than \(\pi\), for example:

>>> exp(sin(1) + E**pi - I)
 sin(1) + ℯ  - ⅈ

>>> _.evalf()
14059120207.1707 - 21895782412.4995⋅ⅈ


>>> zeta(S(14)/17)

>>> zeta(S(14)/17).evalf()

Note that in SymPy, exp(1) is denoted by capital :exp:E

>>> E.evalf()

Symbolic entities are ignored:

>>> pi*x
>>> _.evalf()

Built-in functions float() and complex() take advantage of evalf():

>>> float(pi)
>>> type(_)
<type 'float'>

>>> float(pi*I)
Traceback (most recent call last):
ValueError: Symbolic value, can't compute

>>> complex(pi*I)
>>> type(_)
<type 'complex'>

The base type for computing with floating point numbers in SymPy is Float. It allows for several flavors of initialization and keeps track of precision:

>>> 2.0
>>> type(_)
<type 'float'>

>>> Float(2.0)
>>> type(_)
<class 'sympy.core.numbers.Float'>

>>> sympify(2.0)
>>> type(_)
<class 'sympy.core.numbers.Float'>

>>> Float("3.14")
>>> Float("3.14e-400")

Notice that the last value is out of range for float:

>>> 3.14e-400

We expected a very small value but not zero. This raises an important issue, because if we try to construct a Float this way, we will still get zero:

>>> Float(3.14e-400)

The only way to fix this is to pass a string argument to Float.

When symbolic mathematics matter?

Consider a univariate function:

\[f(x) = x^{(1 - \log(\log(\log(\log(\frac{1}{x})))))}\]

We would like to compute:

\[\lim_{x \to 0^{+}} f(x)\]

Let’s define the function \(f\) in SymPy:

>>> f = x**(1 - log(log(log(log(1/x)))))
>>> f
      ⎛   ⎛   ⎛   ⎛1⎞⎞⎞⎞
 - log⎜log⎜log⎜log⎜─⎟⎟⎟⎟ + 1
      ⎝   ⎝   ⎝   ⎝x⎠⎠⎠⎠

A very straight forward approach is to “see” how \(f\) behaves on the right hand side of zero. We can try to read the solution from the graph of \(f\):

import matplotlib.pyplot as plt
from sympy.mpmath import plot, log

fig = plt.figure()
axes = fig.add_subplot(111)
axes.set_title(r"Plot of $f(x)$ in $[0, 0.01]$.")

f = lambda x: x**(1 - log(log(log(log(1/x)))))
plot(f, xlim=[0, 0.01], axes=axes)

(Source code, png, hires.png, pdf)


This gives us first hint that the limit might be zero. Of course reading a graph of a function isn’t a very precise method for computing limits. Instead of analyzing the graph of \(f\), we can improve this approach a little by evaluating \(f(x)\) for sufficiently small arguments.

Let’s start with arguments of the form \(x = 10^{-k}\):

>>> f.subs(x, 10**-1).evalf()
0.00114216521536353 + 0.00159920801047526⋅ⅈ
>>> f.subs(x, 10**-2).evalf()
>>> f.subs(x, 10**-3).evalf()
>>> f.subs(x, 10**-4).evalf()
>>> f.subs(x, 10**-5).evalf()
>>> f.subs(x, 10**-6).evalf()
>>> f.subs(x, 10**-7).evalf()
>>> f.subs(x, 10**-8).evalf()

We obtained a decreasing sequence values which suggests that the limit is zero. Let’s now try points of the form \(x = 10^{-10^k}\):

>>> f.subs(x, 10**-10**1).evalf()
>>> f.subs(x, 10**-10**2).evalf()
>>> f.subs(x, 10**-10**3).evalf()

For \(x = 10^{-10^3}\) we got a very peculiar value. This happened because:

>>> 10**-10**3

and the reason for this is that we used Python’s floating point values. Instead we can use either exact numbers or SymPy’s floating point numbers:

>>> Integer(10)**-10**3 != 0
>>> Float(10.0)**-10**3 != 0

Let’s continue with SymPy’s floating point numbers:

>>> f.subs(x, Float(10.0)**-10**1).evalf()
>>> f.subs(x, Float(10.0)**-10**2).evalf()
>>> f.subs(x, Float(10.0)**-10**3).evalf()
>>> f.subs(x, Float(10.0)**-10**4).evalf()
>>> f.subs(x, Float(10.0)**-10**5).evalf()
>>> f.subs(x, Float(10.0)**-10**6).evalf()
>>> f.subs(x, Float(10.0)**-10**7).evalf()
>>> f.subs(x, Float(10.0)**-10**8).evalf()

This time the sequence of values is rapidly decreasing, but only until a sufficiently small numer where \(f\) has an inflexion point. After that, values of \(f\) increase very rapidly, which may suggest that the actual limit is +\inf. It seems that our initial guess is wrong. However, for now we still can’t draw any conclusions about behavior of \(f\), because if we take even smaller numbers we may reach other points of inflection.

The mpmath library implements a function for computing numerical limits of function, we can try to take advantage of this:

>>> from sympy.mpmath import limit as nlimit
>>> F = lambdify(x, f, modules='mpmath')

>>> nlimit(F, 0)
(2.23372778188847e-5 + 2.28936592344331e-8j)

This once again suggests that the limit is zero. Let’s use an exponential distribution of points in nlimit():

>>> nlimit(F, 0, exp=True)
(3.43571317799366e-20 + 4.71360839667667e-23j)

This didn’t help much. Still zero. The only solution to this problem is to use analytic methods. For this we will use limit():

>>> limit(f, x, 0)

which shows us that our initial guess was completely wrong. This nicely shows that solving ill conditioned problems may require assistance of symbolic mathematics system. More about this can be found in Dominic Gruntz’s PhD tesis (, where this problem is explained in detail and an algorithm shown, which can solve this problem and which is implemented in SymPy.


  1. Compute first 55 digits of numerical approximation of \(f(\pi)\).


  2. Read this webcomic. What is the first digit of \(e\) to contain \(999999\)? What is the first digit of \(\pi\) to contain \(789\)?


  3. In addition to the above example, Gruntz gives another example of ill conditioned function in his thesis to show why symbolic computation of limits can be preferred to numerical computation:

    \[\lim_{x \to \infty}{\left(\operatorname{erf}\left(x - {e^{-e^{x}}}\right) - \operatorname{erf}\left(x\right)\right) e^{e^{x}}} e^{x^{2}}\]

    (in SymPy, (erf(x - exp(-exp(x))) - erf(x))*exp(exp(x))*exp(x**2)). Compute the above limit in SymPy using methods similar to the ones presented in this section. What are the drawbacks of computing this limit numerically? What is the limit, exactly?


Table Of Contents

Previous topic

Deriving trigonometric identities

Next topic

Summing roots of polynomials

This Page