Numerics Module

Aaron Meurer edited this page Mar 12, 2011 · 2 revisions
Clone this wiki locally

Table of Contents


The module `sympy.numerics` implements arbitrary-precision arithmetic. It has not yet been fully integrated into SymPy.

_Note: the `numerics` module is currently being reimplemented as a standalone Python library, mpmath._

Floats and ComplexFloats

The classes `Float` and `ComplexFloat` behave roughly like Python's `float` and `complex`, except for being capable of arbitrary precision. `Float`s and `ComplexFloat`s can be created from regular Python numbers (`int`, `float`, `complex`), decimal strings, and also from SymPy `Rational`s.

>>> from sympy import *
>>> from sympy.numerics import *
>>> Float(1) - Float(0.25) - Float('0.25') - Float(Rational(1,4))
>>> ComplexFloat(2, 3)
ComplexFloat(real='2', imag='3')
>>> ComplexFloat(1j) + 2
ComplexFloat(real='2', imag='1')

The working precision can be controlled through the static `setdps` method on the `Float` class (the default is 15 decimal places, which corresponds to 53 bits -- the same as Python's `float`).

>>> 2**0.5
>>> Float(2)**0.5

>>> Float.setdps(10)
>>> Float(2)**0.5

>>> Float.setdps(100)
>>> Float(2)**0.5

`Float` also supports various rounding modes: `ROUND_DOWN`, `ROUND_UP`, `ROUND_FLOOR`, `ROUND_CEILING`, `ROUND_HALF_UP`, `ROUND_HALF_DOWN`, and `ROUND_HALF_EVEN` (default). However, these rounding modes are not yet correctly supported in arithmetic operations. For more details on managing the working precision and rounding modes, refer to the docstrings for the `Float` class.

The method `ae` can be used for approximate equality testing:

>>> Float.setdps(15)
>>> Float(0).ae(1e-14)
>>> Float(0).ae(1e-16)

Evaluating SymPy expressions

The `evalf` function evaluates exact SymPy expressions into `Float`s or `ComplexFloat` approximations. It currently supports regular algebraic expressions and simple functions and constants like `pi`, `exp`, etc.

>>> evalf((sqrt(pi) + 2*I)**Rational(-1,3))
ComplexFloat(real='0.69217119929576842', imag='-0.20044721328786239')

>>> Float.setdps(50)
>>> evalf((1+sqrt(5))/2)

Special functions

The modules `numerics.functions` and `numerics.functions2` provide mathematical functions. `functions` contains standard functions like those found in the Python `math` library while `functions2` contains higher transcendental functions. Some mathematical constants are also available in the `numerics.constants` module.

The functions are incompatible with SymPy's symbolic functions (so be careful with the imports).

All available functions and constants can be computed with arbitrary precision. Most functions also support both real and complex numbers. For example, this computes a complex exponential with standard precision:

>>> print exp(1j)
(0.540302305868140 + 0.841470984807897*I)

As another example, the Riemann zeta function value zeta(3) can easily be computed to, say, 1000 digits:

>>> from sympy.numerics.functions2 import *
>>> Float.setdps(1000)
>>> print zeta(3)

A complex value of the gamma function:

>>> Float.setdps(100)
>>> print gamma(4+4j)
56803721515698482878208 + -0.496739083997423706378466662866211303078212882210124

The constant pi is available through a function `pi_float()` that returns an approximation of pi accurate to within the current working precision. The `examples` directory also contains a script `` for computing pi with extremely high precision (millions of digits for the patient).

Numerical optimization

The module `numerics.optimize` contains functions for root-finding and optimization.

The functions `bisect` and `secant` functions can be used to find isolated roots of decently well-behaved functions. As a simple example; using the secant method to calculate pi as the root of sin(x) closest to x = 3:

>>> from sympy.numerics.optimize import secant
>>> from sympy.numerics.functions import sin
>>> Float.setdps(50)
>>> print secant(sin, 3)

With details:

>>> print secant(sin, 3, verbose=True)
  secant: x_1=3.000000  delta=-0.03
  secant: x_2=3.142264  delta=0.142264
  secant: x_3=3.141590  delta=-0.0006733
  secant: x_4=3.141593  delta=2.23687e-006
  secant: x_5=3.141593  delta=-1.67327e-013
  secant: x_6=3.141593  delta=1.39539e-025
  secant: x_7=3.141593  delta=0

More examples and details on usage are available in the function docstrings.

The function `polyroots` finds all roots of a polynomial using complex arithmetic.

>>> from sympy.numerics.optimize import polyroots
>>> x = Symbol('x')
>>> roots, error = polyroots(x**3 - 2*x**2 + x + 2)
>>> for r in roots: print r
(-0.695620769559862 + 6.78546596234906E-24*I)
(1.34781038477993 + -1.02885225413669*I)
(1.34781038477993 + 1.02885225413669*I)
>>> print error

This polynomial is very ill-conditioned:

>>> x = Symbol('x')
>>> W = 1
>>> for i in range(1, 14):
...     W = W * (x-i)
>>> W = W.expand()
>>> roots, error = polyroots(W)
>>> for r in roots: print r
(1 + 0*I)
(2.00000000000052 + -6.78706583304582E-13*I)
(3.00000000091542 + 2.03402685728453E-9*I)
(3.99967891066583 + -0.000901013608686283*I)
(5.00026649719202 + -0.0000589434738259911*I)
(6.07300960794678 + 0.788765735425133*I)
(6.56514846614468 + 0.231277402207583*I)
(8.05562929719922 + -0.512541233266305*I)
(8.52882344788961 + 1.04323424522302*I)
(10.1764629806671 + -0.824843997843889*I)
(11.1117004516702 + 1.06172142342246*I)
(12.1969006408498 + -0.338168434918473*I)
(13.0738393202625 + 0.116644127674871*I)
>>> print error

In this case, calling `polyroots` with a higher `maxsteps` setting will give better results. It may also be necessary to increase the working precision.

Numerical integration

The function `nintegrate` can be used to compute integrals numerically. It takes as input a regular Python function that operates on `Float`s (or `ComplexFloat`s), rather than a SymPy expression.

A simple integral, sin(x) between 0 and pi (`sin` and `pi_float` must be imported from `numerics.functions` and `numerics.constants`):

>>> nintegrate(lambda x: sin(x), 0, pi_float())

The `verbose` option shows what `nintegrate` does internally:

>>> nintegrate(lambda x: sin(x), 0, pi_float(), verbose=True)
calculating nodes for degree-13 Gauss-Legendre quadrature...
  node 0 of 6
  node 4 of 6
calculating nodes for degree-6 Gauss-Legendre quadrature...
  node 0 of 3
0 0 3.14159265358979312
1 0 1.57079632679489656
2 0 0.785398163397448279
5 0.785398163397448279 1.57079632679489656
8 1.57079632679489656 3.14159265358979312
9 1.57079632679489656 2.35619449019234484
12 2.35619449019234484 3.14159265358979312

By default, `nintegrate` uses Gauss-Legendre quadrature. The degree of the quadrature rule depends on the precision level. Since computing nodes for the quadrature rule is a relatively expensive operation (often more expensive than evaluating the integral), computed weights are cached to make repeated calls to `nintegrate` faster.

The interval is adaptively divided into smaller parts to obtain high accuracy. In the example above, we see a count of the number of times the quadrature rule has been applied and the boundaries of each subinterval. In this case, a dozen applications were sufficient (a small piece of sine wave is an example of a very well-behaved integrand). If the integrand is less well-behaved and/or the working precision is set very high, `nintegrate` may require thousands of evaluations, or it might fail altogether. The `nintegrate` docstrings contain more information on how to deal with difficult integrals.

Here is an integral for Euler's constant (infinite integration intervals are supported):

>>> nintegrate(lambda x: -exp(-x)*log(x), 0, oo)

It turns out that the standard Gaussian quadrature rule performs poorly on this integral when the precision is increased. If we choose the tanh-sinh method (`method=1`) instead, it only takes a few seconds to evaluate it with 100 digits:

>>> Float.setdps(100)
>>> nintegrate(lambda x: -exp(-x)*log(x), 0, oo, method=1)

As verification, there is a function that computes Euler's constant in the `numerics.constants` module (it uses a more efficient algorithm):

>>> gamma_float()

If `nintegrate` fails to reach full accuracy, it prints a warning message along with an estimate of the magnitude of the error.

Here is a very badly behaved integral where `nintegrate` fails altogether (the correct value is pi/2):

>>> nintegrate(lambda x: sin(x)/x, 0, oo, maxsteps=100)
Warning: failed to reach full accuracy. Estimated magnitude of error:  2

Increasing the number of steps only makes things worse:

>>> nintegrate(lambda x: sin(x)/x, 0, oo, maxsteps=1000)
Warning: failed to reach full accuracy. Estimated magnitude of error:  6

A better approach in this case is to choose a finite upper limit:

>>> nintegrate(lambda x: sin(x)/x, 0, 100, maxsteps=1000)

Here is a similar but slightly easier integral (the correct value is roughly 1.253).

>>> nintegrate(lambda x: sin(x**2)/x**2, 0, oo, maxsteps=1000)
Warning: failed to reach full accuracy. Estimated magnitude of error:  0.003