[SciPy-User] Numerical precision: Scipy Vs Numpy

[Sergio Rojas]

"""
 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