3

我正在使用 Python 2 和维基百科文章“三次函数”中给出的相当简单的方法。这也可能是我必须定义的立方根函数的问题，以便创建标题中提到的函数。

# Cube root and cubic equation solver
#
# Copyright (c) 2013 user2330618
#
# This Source Code Form is subject to the terms of the Mozilla Public
# License, v. 2.0. If a copy of the MPL was not distributed with this
# file, you can obtain one at http://www.mozilla.org/MPL/2.0/.

from __future__ import division
import cmath
from cmath import log, sqrt

def cbrt(x):
    """Computes the cube root of a number."""
    if x.imag != 0:
        return cmath.exp(log(x) / 3)
    else:
        if x < 0:
            d = (-x) ** (1 / 3)
            return -d
        elif x >= 0:
            return x ** (1 / 3)

def cubic(a, b, c, d):
    """Returns the real roots to cubic equations in expanded form."""
    # Define the discriminants
    D = (18 * a * b * c * d) - (4 * (b ** 3) * d) + ((b ** 2) * (c ** 2)) - \
    (4 * a * (c ** 3)) - (27 * (a ** 2) * d ** 2)
    D0 = (b ** 2) - (3 * a * c)
    i = 1j  # Because I prefer i over j
    # Test for some special cases
    if D == 0 and D0 == 0:
        return -(b / (3 * a))
    elif D == 0 and D0 != 0:
        return [((b * c) - (9 * a * d)) / (-2 * D0), ((b ** 3) - (4 * a * b * c)
        + (9 * (a ** 2) * d)) / (-a * D0)]
        else:
            D1 = (2 * (b ** 3)) - (9 * a * b * c) + (27 * (a ** 2) * d)
            # More special cases
            if D != 0 and D0 == 0 and D1 < 0:
                C = cbrt((D1 - sqrt((D1 ** 2) - (4 * (D0 ** 3)))) / 2)
            else:
                C = cbrt((D1 + sqrt((D1 ** 2) - (4 * (D0 ** 3)))) / 2)
                u_2 = (-1 + (i * sqrt(3))) / 2
                u_3 = (-1 - (i * sqrt(3))) / 2
                x_1 = (-(b + C + (D0 / C))) / (3 * a)
                x_2 = (-(b + (u_2 * C) + (D0 / (u_2 * C)))) / (3 * a)
                x_3 = (-(b + (u_3 * C) + (D0 / (u_3 * C)))) / (3 * a)
                if D > 0:
                    return [x_1, x_2, x_3]
                else:
                    return x_1

我发现这个函数能够解决一些简单的三次方程：

print cubic(1, 3, 3, 1)
-1.0

不久前，我已经达到了可以求解具有两个根的方程的程度。我刚刚做了一个重写，现在它变得混乱了。例如，这些系数是 (2x - 4)(x + 4)(x + 2) 的扩展形式，它应该返回 [4.0, -4.0, -2.0] 或类似的东西：

print cubic(2, 8, -8, -32)
[(-4+1.4802973661668753e-16j), (2+2.9605947323337506e-16j), (-2.0000000000000004-1.1842378929335002e-15j)]

这更像是我正在犯的数学或编程错误吗？

更新：谢谢大家的回答，但是这个函数的问题比我迄今为止迭代的要多。例如，我经常收到与立方根函数有关的错误：

print cubic(1, 2, 3, 4)  # Correct solution: about -1.65
...
    if x > 0:
TypeError: no ordering relation is defined for complex numbers
print cubic(1, -3, -3, -1)  # Correct solution: about 3.8473
    if x > 0:
TypeError: no ordering relation is defined for complex numbers
4

3 回答 3

7

Wolfram Alpha 确认你最后一个三次方的根确实是

(-4, -2, 2)

而不是你说的

...它应该返回[4.0, -4.0, -2.0]

尽管有那个（我认为）错字，你的程序给出了

[(-4+1.4802973661668753e-16j), (2+2.9605947323337506e-16j), (-2.0000000000000004-1.1842378929335002e-15j)]

准确度与正确解决方案的10**(-15)根完全相同。正如其他人所说，微小的虚部可能是由于四舍五入。

请注意，如果您使用像Cardano's这样的解决方案，您必须使用精确的算法才能始终正确取消。MAPLE这是程序喜欢或Mathematica存在的原因之一，通常从公式到实现之间存在脱节。

要在纯 python 中只获取数字的实部，请调用.real. 例子：

a = 3.0+4.0j
print a.real
>> 3.0
于 2013-04-29T05:16:51.043 回答
3

如果您想以数字方式执行此操作，那么Hooked 的答案就是要走的路。您也可以使用sympy象征性地执行此操作：

>>> from sympy import roots
>>> roots('2*x**3 + 8*x**2 - 8*x - 32')
{2: 1, -4: 1, -2: 1}

这为您提供了根源及其多样性。

于 2013-04-29T15:18:30.813 回答
-5

您正在使用整数值 - Python 不会自动将其转换为浮点数。更通用的解决方案是将函数中的系数写为浮点数 - 18.0 而不是 18 等。这将起到一个插图 - 来自代码：

>>> 2**(1/3)
1
>>> 2**(1/3.)
1.2599210498948732
>>>
于 2013-04-29T03:58:34.467 回答