[SciPy-User] Numerical precision: Scipy Vs Numpy

[Evgeni Burovski]

On Wed, Nov 4, 2015 at 12:58 PM, Sergio Rojas <sergio_r at mail.com> wrote:
> """
>  Having fun with numerical precision/errors, a solution of the quadratic
> equation [a*x**2 + b*x + c == 0] with small values for the
> constant a is given an unexpected result when comparing
> SciPy and SymPy numerical results, which I hope an earthly kind soul could
> explain. The output is as follows (the bothering result
> is on the second root property):
>
>          SciPy
> a*x1**2 + b*x1 + c =  16777217.0
> a*x2**2 + b*x2 + c =  1.11022302463e-16
>          Root properties:
> (x1 + x2) + b/a =  0.0
>   x1*x2   - c/a =  0.0
>          ----
>          SymPy
> a*x1**2 + b*x1 + c =  16777217.0000000
> a*x2**2 + b*x2 + c =  0
>          Root properties:
> (x1 + x2) + b/a =  0
>   x1*x2   - c/a =  0.125000000000000
> """
>
> import scipy
>
> def myprint(a, b, c, x1, x2):
>   print('a*x1**2 + b*x1 + c = '),
>   print(a*x1**2 + b*x1 + c)
>   print('a*x2**2 + b*x2 + c = '),
>   print(a*x2**2 + b*x2 + c)
>   print("\t Root properties:   ")
>   print('(x1 + x2) + b/a = '),
>   print((x1+x2) + float(b)/float(a))
>   print('  x1*x2   - c/a = '),
>   print(x1*x2 - float(c)/float(a))
>
> a = 1.0e-15
> b = 10000.0
> c = 1.0
>
> coeff = [a, b, c]
> x1, x2 = scipy.roots(coeff)
>
> print("\t SciPy   ")
> myprint(a, b, c, x1, x2)
>
> from sympy import *
> var('x')
> x1, x2 = solve(a*x**2 + b*x + c, x)
>
> print(" \t ---- \n \t SymPy   ")
> myprint(a, b, c, x1, x2)
>
> # Thanks in advance,
> # Sergio
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> https://mail.scipy.org/mailman/listinfo/scipy-user

This is a catastrophic loss of precision.
For a way out, see, for example,
http://math.stackexchange.com/questions/866331/how-to-implement-a-numerically-stable-solution-of-a-quadratic-equation

Demo:

In [1]: a, b, c = 1e-15, 1e4, 1.0

In [2]: D = b**2 - 4.*a*c

In [3]: import numpy as np

In [4]: 2.*c / (-b - np.sqrt(D))
Out[4]: -0.0001

In [5]: import sympy as sy

In [6]: x = sy.var('x')

In [7]: sy.solve(a*x**2 + b*x + c)
Out[7]: [-1.00000000000000e+19, -0.000100000000000000]