c# - 对三次方程解实现 Cardano 方法

Question

我正在尝试通过使用卡尔达诺方法找到由一组四个系数定义的三次方程的实根，如 here 所述.问题是，我的实现找到的根实际上不起作用——通过将它们插入等式进行测试会产生重大错误(超过所需的 10^-6)。是算法实现错误，还是其他原因导致错误，例如舍入精度？

static double CubicRoot(double n)
    {
        return Math.Pow(Math.Abs(n), 1d / 3d) * Math.Sign(n);
    }


public static List<double> SolveCubic(double A, double B = 0, double C = 0, double D = 0)
    {
        List<double> output = new List<double>();
        if (A != 0)
        {

                double A1 = B / A;
                double A2 = C / A;
                double A3 = D / A;
                double P = -((A1 * A1) / 3) + A2;
                double Q = ((2.0 * A1 * A1 * A1) / 27.0) - ((A1 * A2) / 3.0) + A3;
                double cubeDiscr = Q * Q / 4.0 + P * P * P / 27.0;
                if (cubeDiscr > 0)
                {
                    double u = CubicRoot(-Q / 2.0 + Math.Sqrt(cubeDiscr));
                    double v = CubicRoot(-Q / 2.0 - Math.Sqrt(cubeDiscr));
                    output.Add(u + v - (A1 / 3.0));
                    return output;
                }
                else if (cubeDiscr == 0)
                {
                    double u = CubicRoot(-Q / 2.0);
                    output.Add(2u - (A1 / 3.0));
                    output.Add(-u - (A1 / 3.0));
                }
                else if (cubeDiscr < 0)
                {
                    double r = CubicRoot(Math.Sqrt(-(P * P * P / 27.0)));
                    double alpha = Math.Atan(Math.Sqrt(-cubeDiscr) / (-Q / 2.0));
                    output.Add(r * (Math.Cos(alpha / 3.0) + Math.Cos((6 * Math.PI - alpha) / 3.0)) - A1 / 3.0);
                    output.Add(r * (Math.Cos((2 * Math.PI + alpha) / 3.0) + Math.Cos((4 * Math.PI - alpha) / 3.0)) - A1 / 3.0);
                    output.Add(r * (Math.Cos((4 * Math.PI + alpha) / 3.0) + Math.Cos((2 * Math.PI - alpha) / 3.0)) - A1 / 3.0);
                }
        }
        return output;
    }

Answer 1

一些事情

• Math.Sign将以零返回零，在这种情况下这恰好是您想要的，但也许您对代码或算法更改没有那么幸运。
• 您将遇到舍入问题并且不会执行 cubeDiscr == 0在你应该的时候分支。您可能有舍入问题并执行错误的 > 0和 < 0出于同样的原因分支。改为在零增量内进行测试(见下文)。
• 但是cubeDiscr == 0分支是错误的，因为 1) 你没有计算 v和 2) 2u是一个 UInt32值为 2，而不是 2*u .
• 计算 alpha 错误(见下文)
• (可能还有更多，但我很快就看到了这些)

关于计算 alpha:

double alpha = Math.Atan(Math.Sqrt(-cubeDiscr) / (-Q / 2.0));

不一样

double alpha = Math.Atan(Math.Sqrt(-d) / q * 2.0);
if (q > 0)                         // if q > 0 the angle becomes  PI + alpha
    alpha = Math.PI + alpha;

使用该页面中包含的代码有什么问题？

public double Xroot(double a, double x)
{
    double i = 1;
    if (a < 0)
        i = -1;
    return (i * Math.Exp( Math.Log(a*i)/x));
}

public int Calc_Cardano()  // solve cubic equation according to cardano
{
    double p, q, u, v;
    double r, alpha;
    int res;
    res = 0;
    if (a1 != 0)
    {
        a = b / a1;
        b = c / a1;
        c = d / a1;

        p = -(a * a / 3.0) + b;
        q = (2.0 / 27.0 * a * a * a) - (a * b / 3.0) + c;
        d = q * q / 4.0 + p * p * p / 27.0;
        if (Math.Abs(d) < Math.Pow(10.0, -11.0))
            d = 0;
        // 3 cases D > 0, D == 0 and D < 0
        if (d > 1e-20)
        {
            u = Xroot(-q / 2.0 + Math.Sqrt(d), 3.0);
            v = Xroot(-q / 2.0 - Math.Sqrt(d), 3.0);
            x1.real = u + v - a / 3.0;
            x2.real = -(u + v) / 2.0 - a / 3.0;
            x2.imag = Math.Sqrt(3.0) / 2.0 * (u - v);
            x3.real = x2.real;
            x3.imag = -x2.imag;
            res = 1;
        }
        if (Math.Abs(d) <= 1e-20)
        {
            u = Xroot(-q / 2.0, 3.0);
            v = Xroot(-q / 2.0, 3.0);
            x1.real = u + v - a / 3.0;
            x2.real = -(u + v) / 2.0 - a / 3.0;
            res = 2;
        }
        if (d < -1e-20)
        {
            r = Math.Sqrt(-p * p * p / 27.0);
            alpha = Math.Atan(Math.Sqrt(-d) / q * 2.0);
            if (q > 0)                         // if q > 0 the angle becomes  PI + alpha
                alpha = Math.PI + alpha;

            x1.real = Xroot(r, 3.0) * (Math.Cos((6.0 * Math.PI - alpha) / 3.0) + Math.Cos(alpha / 3.0)) - a / 3.0;
            x2.real = Xroot(r, 3.0) * (Math.Cos((2.0 * Math.PI + alpha) / 3.0) + Math.Cos((4.0 * Math.PI - alpha) / 3.0)) - a / 3.0;
            x3.real = Xroot(r, 3.0) * (Math.Cos((4.0 * Math.PI + alpha) / 3.0) + Math.Cos((2.0 * Math.PI - alpha) / 3.0)) - a / 3.0;
            res = 3;
        }
    }
    else
        res = 0;
    return res;
 }