# [SciPy-User] Numerical precision: Scipy Vs Numpy Classic List Threaded 4 messages Open this post in threaded view
|

## [SciPy-User] Numerical precision: Scipy Vs Numpy

 """  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 [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user
Open this post in threaded view
|

## Re: Numerical precision: Scipy Vs Numpy

 On Wed, Nov 4, 2015 at 12:58 PM, Sergio Rojas <[hidden email]> 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 > [hidden email] > https://mail.scipy.org/mailman/listinfo/scipy-userThis 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-equationDemo: In : a, b, c = 1e-15, 1e4, 1.0 In : D = b**2 - 4.*a*c In : import numpy as np In : 2.*c / (-b - np.sqrt(D)) Out: -0.0001 In : import sympy as sy In : x = sy.var('x') In : sy.solve(a*x**2 + b*x + c) Out: [-1.00000000000000e+19, -0.000100000000000000] _______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user