[SciPy-User] Numerical precision: Scipy Vs Numpy

[Sergio Rojas]

 Following up on my question, thanks to
Evgeni Burovski [ https://mail.scipy.org/pipermail/scipy-user/2015-November/036753.html ] and
David [ https://mail.scipy.org/pipermail/scipy-user/2015-November/036754.html ]

I just want to add that on Fortran 90, via real*16 (also called
quad precision) the output for the roots are (not
sure about the rearrangement proposed by David in this case):

         FORTRAN 90 (via real*16)

  a =   1.00000000362749372553872184710144211E-0015
  b =    10000.000000000000000000000000000000
  c =    1.0000000000000000000000000000000000

  x1 =   -9999999963725062876.1997883398821330
  x2 =  -1.00000000000000000000001000000003626E-0004

 a*x1**2 + b*x1 + c =    0.0000000000000000000000000000000000
 (a*x1 + b)*x1 + c =  -1.05518034299496531738542345583594055E-0011

 a*x2**2 + b*x2 + c =    0.0000000000000000000000000000000000
 (a*x2 + b)*x2 + c =    0.0000000000000000000000000000000000

 (x1 + x2) + b/a =    0.0000000000000000000000000000000000
   x1*x2   - c/a =    0.0000000000000000000000000000000000

Sergio
 

Sent: Wednesday, November 04, 2015 at 1:58 PM
From: "Sergio Rojas" <sergio_r at mail.com>
To: scipy-user at scipy.org
Subject: 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