Search code examples
pythonnumericnumerical-methodsequationquadratic

Python - Quadratic Equatin Extreme numbers

I have written code to solve quadratic equation in Python. But it has errors with very large numbers like "1e10" and "-1e10". Is there any solution numerically or as python solution?

enter code here
import math
import cmath

def solve_quad(b, c):

    a = 1
    D = b*b - 4*a*c

    if (D > 0):
        x1 = (-b - math.sqrt(D)) / (2*a)
        x2 = (-b + math.sqrt(D)) / (2*a)
    elif D == 0:
        x1 = (-b - math.sqrt(D)) / (2*a)
        x2 = x1
    else:
        x1 = (-b - cmath.sqrt(D)) / (2*a)
        x2 = (-b + cmath.sqrt(D)) / (2*a)


    return x1, x2
print(solve_quad(1e10, 4))

Output: (-10000000000.0, 0.0)

Solution

  • This is a very old and often repeated topic. Let's explain it one more time. You use the binomial theorem to avoid numerical instability.

    As the first root you chose the one where the sign of -b and the sign before the square root are the same.

    if D>0:
    SD = D**0.5;
    if b > 0: SD = -SD
    x1 = (-b+SD) / (2*a)

    Then for the second root you use the formula

    (-b-SD) / (2*a) 
= (b^2-SD^2) / (2*a*(-b+SD)) 
= 4*a*c / (2*a*(-b+SD)) 
= (2*c) / (-b+SD)

    to get

    x2 = (2*c) / (-b+SD)

    In the other cases the catastrophic cancellation that is avoided with this procedure does not occur.

    This avoids completely all numerical instability due to catastrophic cancellation. If you want go further you can also try to avoid the potential overflow in the computation of the discriminant.