该函数解三次方程式有什么问题? [英] What's wrong with this function to solve cubic equations?

问题描述

我正在使用Python 2和Wikipedia的文章三次函数"中给出的相当简单的方法.这也可能是我必须定义才能创建标题中提到的功能的多维数据集根函数的问题.

I am using Python 2 and the fairly simple method given in Wikipedia's article "Cubic function". This could also be a problem with the cube root function I have to define in order to create the function mentioned in the title.

# 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

我发现此函数能够求解一些简单的三次方程式:

I've found that this function is capable of solving some simple cubic equations:

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

前一段时间，我已经了解到它可以求解具有两个根的方程式.我刚刚进行了重写，但现在已经一去不复返了.例如，这些系数是(2x-4)(x + 4)(x + 2)的扩展形式，并且应返回[4.0，-4.0，-2.0]或类似的值:

And a while ago I had gotten it to a point where it could solve equations with two roots. I've just done a rewrite and now it's gone haywire. For example, these coefficients are the expanded form of (2x - 4)(x + 4)(x + 2) and it should return [4.0, -4.0, -2.0] or something similar:

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

这是我正在犯的数学错误还是编程错误?

Is this more a mathematical or a programming mistake I'm making?

更新:谢谢大家，您的回答是正确的，但是到目前为止，此功能存在的问题比我迭代的更多.例如，我经常收到与多维数据集根函数有关的错误:

Update: Thank you, everyone, for your answers, but there are more problems with this function than I've iterated so far. For example, I often get an error relating to the cube root function:

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

推荐答案

Wolfram Alpha confirms that the roots to your last cubic are indeed

(-4, -2, 2)

而不是您所说的

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

尽管我猜错了，但您的程序却给出了

Not withstanding that (I presume) typo, your program gives

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

完全相同的根源是10**(-15)的正确解决方案.正如其他人所说，微小假想部分可能是由于四舍五入.

Which to accuracy of 10**(-15) are the exact same roots as the correct solution. The tiny imaginary part is probably due, as others have said, to rounding.

请注意，如果您使用的是 Cardano的解决方案，则必须始终使用精确算术来始终正确取消. .这就是为什么像MAPLE或Mathematica这样的程序存在的原因之一，通常是从公式到实现之间存在脱节.

Note that you'll have to use exact arithmetic to always correctly cancel if you are using a solution like Cardano's. This one of the reasons why programs like MAPLE or Mathematica exist, there is often a disconnect from the formula to the implementation.

要仅获取纯python中数字的实部，请调用.real.示例:

To get only the real portion of a number in pure python you call .real. Example:

a = 3.0+4.0j
print a.real
>> 3.0

这篇关于该函数解三次方程式有什么问题?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！