Python-平方Equatin极限数

Question

我已经编写了代码来解决Python中的二次方程式。但是它有很大的错误，例如“ 1e10”和“ -1e10”。有数字解决方案还是python解决方案？

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

输出：（-10000000000.0，0.0）

Answer 1

这是一个非常古老且经常重复的主题。让我们再解释一次。使用二项式定理可以避免数值不稳定。

作为第一个根，您选择了-b的符号和平方根之前的符号相同的那个。

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

然后使用公式的第二个根

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

获取

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

在其他情况下，不会发生通过此过程避免的灾难性取消。

这完全避免了由于灾难性抵消而造成的所有数值不稳定。如果您想走得更远，也可以尝试避免判别式计算中可能出现的溢出。

Answer 2

您可能会遇到浮点精度问题。您可以使用Decimal或其他类似的库来获得更高的精度并规避此问题：

from decimal import *

def solve_quad(b, c):
    a = 1
    D = b*b - 4*a*c

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

    return x1, x2

print(solve_quad(Decimal(1e10), Decimal(4)))