0

I'm pretty new to python and am trying to write some code to solve a given quadratic function. I'm having some trouble with rounding errors in floats, I think because I am dividing two numbers that are very large with a very small difference. (Also I'm assuming all inputs have real solutions for now.) I've put two different versions of the quadratic equation to show my problem. It works fine for most inputs, but when I try a = .001, b = 1000, c = .001 I get two answers that have a significant difference. Here is my code:

from math import sqrt

a = float(input("Enter a: "))
b = float(input("Enter b: "))
c = float(input("Enter c: "))

xp = (-b+sqrt(b**2-4*a*c))/(2*a)
xn = (-b-sqrt(b**2-4*a*c))/(2*a)

print("The solutions are: x = ",xn,", ",xp,sep = '')

xp = (2*c)/(-b-sqrt(b**2-4*a*c))
xn = (2*c)/(-b+sqrt(b**2-4*a*c))

print("The solutions are: x = ",xn,", ",xp,sep = '')
inspectorG4dget
  • 110,290
  • 27
  • 149
  • 241
user1934956
  • 1
  • 1
  • 1
  • 2
    What results do you get and how do they differ from what you expect? – inspectorG4dget Dec 28 '12 at 16:57
  • When run with the given values the program prints: The solutions are: x = -999999.999999, -9.999894245993346e-07 The solutions are: x = -1000010.5755125057, -1.000000000001e-06 I would expect the values to be the same... – user1934956 Dec 28 '12 at 17:18
  • In general, you cannot get the exact solutions due to rounding errors, as you noticed. But what is your question? How can you get the exact solutions? Or is it really rounding errors that mess up your results? – Ali Dec 28 '12 at 18:16

5 Answers5

1

I'm no expert in the maths field but I believe you should use numpy (a py module for maths), due to internal number representation on computers your calculus will not match real math. (floating point arithmetics)

http://docs.python.org/2/tutorial/floatingpoint.html

Check this is almost exaclty what you want.

http://www.annigeri.in/2012/02/python-class-for-quadratic-equations.html

MGP
  • 2,981
  • 35
  • 34
1

To get more precise results with floating point, be careful not to subtract similar quantities. For the quadratic x^2 + a x + b = 0 you know that the roots x1 and x2 make

b = x1 * x2

Compute the one with larger absolute value, and get the other one from this relation.

CassOnMars
  • 6,153
  • 2
  • 32
  • 47
vonbrand
  • 11,412
  • 8
  • 32
  • 52
0

Solutions:

Numpy as suggested by user dhunter is usually the best solution for math in python. The numpy libraries are capable of doing quick and accurate math in a number of different fields.

Decimal data types were added in python 2.4 If you do not want to download an external library and do not anticipate doing many long or complex equations, decimal datatypes may fit the bill. Simply add: 

from decimal import *

to the top of your code and then replace all instances of the word float with the word Decimal (note the uppercase "D".)

Ex: Decimal(1.1047262519) as opposed to float(1.1047262519)

Theory:

Float arithmetic is based off of binary math and is therefore not always exactly what a user would expect. An excelent description of the float Vs. decimal types is located Here

Dlinet
  • 1,193
  • 3
  • 15
  • 22
0

There are special cases that you should deal with:

  1. a == 0 means a linear equation and one root: x = -c/b
  2. b == 0 means two roots of the form x1, x2 = ±sqrt(-c/a)
  3. c == 0 means two roots, but one of them is zero: x*(ax+b) = 0
  4. If the discriminant is negative, you have two complex conjugate roots.

I'd recommend calculating the discriminant this way: 

discriminant = b*sqrt(1.0-4.0*a*c/b)

I'd also recommend reading this:

https://math.stackexchange.com/questions/187242/quadratic-equation-error

Community
  • 1
  • 1
duffymo
  • 305,152
  • 44
  • 369
  • 561
0

The previously-mentioned numpy module is not particularly relevant to the rounding error mentioned in the question. On the other hand, the decimal module can be used in a brute-force manner to get accurate computations. The following snippet from an ipython interpreter session illustrates its use (with default 28-digit accuracy), and also shows that the corresponding floating-point calculation only has 5 decimal places of accuracy.

In [180]: from decimal import Decimal
In [181]: a=Decimal('0.001'); b=Decimal('1000'); c=Decimal('0.001')
In [182]: (b*b - 4*a*c).sqrt()
Out[182]: Decimal('999.9999999979999999999980000')
In [183]: b-(b*b - 4*a*c).sqrt()
Out[183]: Decimal('2.0000000000020000E-9')
In [184]: a = .001; b = 1000; c = .001
In [185]: math.sqrt(b*b - 4*a*c)
Out[185]: 999.999999998
In [186]: b-math.sqrt(b*b - 4*a*c)
Out[186]: 1.999978849198669e-09
In [187]: 2*a*c/b
Out[187]: 1.9999999999999997e-09

Taylor series for the square root offers an alternative method to use when 4ac is tiny compared to b**2. In this case, √(b*b-4*a*c) ≈ b - 4*a*c/(2*b), whence b - √(b*b-4*a*c) ≈ 2*a*c/b. As can be seen in the line [187] entries above, Taylor series computation gives a 12-digits-accurate result while using floating point instead of Decimal. Using another Taylor series term might add a couple more digits of accuracy.

James Waldby - jwpat7
  • 8,593
  • 2
  • 22
  • 37