4 ms·
What you've done here is tell SymPy to use extra precision for the intermediate (and final) output. This doesn't truly fix the problem of cancellation and loss
by Majromax 3y ago
What you've done here is tell SymPy to use extra precision for the intermediate (and final) output. This doesn't truly fix the problem of cancellation and loss of precision, but for many practical purposes it can postpone the problem long enough to give you a useful result.
Internally, SymPy uses mpmath (https://mpmath.org/ https://mpmath.org/) for representation of numbers to arbitrary precision. You could install and use the latter library directly, gaining extra precision without going through symbolic manipulation.
All that being said, it's still good practice to avoid loss of precision at the outset. Arbitrary-precision calculations are slow compared to hardware-native floating point operations. Using the example from mpmath's homepage in iPython:
In [1]: import mpmath as mp; import scipy as sp; import numpy as np
In [2]: mp.mp.dps=50 # set extended precision
In [3]: %%timeit
...: mp.quad(lambda x: mp.exp(-x**2), [-mp.inf, mp.inf]) ** 2
40.5 ms ± 4.74 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
In [4]: mp.quad(lambda x: mp.exp(-x**2), [-mp.inf, mp.inf]) ** 2
Out[4]: mpf('3.1415926535897932384626433832795028841971693993751015')
In [5]: %%timeit
...: sp.integrate.quad(lambda x: np.exp(-x**2), -np.inf, np.inf)[0]**2
570 µs ± 78.9 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
In [6]: sp.integrate.quad(lambda x: np.exp(-x**2), -np.inf, np.inf)[0]**2
Out[6]: 3.1415926535897927