在编写张学文老师的“时间胎”演示程序中，卡壳了。

卡在求解根据屏幕坐标反投影求模型坐标时，要解一个一元四次方程。

看起来并不复杂的一元四次方程，用迭代法很容易求解。但为了提高效率，想找个公式解法。

 

网上查到所谓的通用公式解法有“费拉里”法。

基本想法是将一元四次方程的表达式配成两个完全平方表达式的等式，为了实现这个目的，需要在此之前解一个一元三次方程，得到配方的系数。方法是不错，程序也编出来了，就是实际使用中怎么也得不到正确的解。解法如下：

一般一元四次方程为：

x^4+c1x^3+c2x^2+c3x+c4=0

移项后：

左边为一个完全平方式，右边为一个x的二次多项式(系数中含有y),

为了使右边也成为一个完全平方式，必须使其判别式为零，即：

DAT = [(c1y/2-c3)]^2 - 4 * [(c1^2/4-c2+y)] * [(y/2)^2-c4] = 0

于是，得到一个关于y的一元三次方程。

-y^3 + c2y^2 + (4c4-c1c3)y+c3^2+c4c1^2-4c4c2=0

 

只要解这个一元三次方程得到一个解y0，代入（1）右边，并令右边=0解出配方常量x0;

便可得到两边都为完全平方的方程，原一元四次方程被降为一元二次方程。可轻松求解。

((x^2+c1x/2）+ y0/2)^2 =(x-x0)^2

即：

(x^2+c1x/2）+ y0/2 = x-x0

整理为：

原码核心部分如下：

 

Procedure T1P4Func.Calculate; //一元四次方程求解计算方法
  var
    Fa,Fb,Fc,Fd,Fe,C4 : Real; //系数计算中间变量
    Za,Zb,Zc,zd,Ze : Real;     //系数计算中间变量
    AFunc : T1P3Func;          //一元三次方程实例，用来计算配方系数
    BFunc,CFunc,DFunc : T1P2Func;  //一元二次方程实例，用来求解配方解和最终解
  begin
    Fa := Coefficients[0];  //按通式ax^4+bx^3+cx^2+dx+e=0获取系数
    Fb := Coefficients[1];
    Fc := Coefficients[2];
    Fd := Coefficients[3];
    Fe := Coefficients[4];
    Zb := Fb/Fa;                     //化简为通式：x^4+zbx^3+zcx^2+zdx+ze=0
    Zc := Fc/Fa;
    Zd := Fd/Fa;
    Ze := Fe/Fa;
    Za := Zb*Zb;
    AFunc := T1P3Func.Create(3);     //建立配方系数求解方程,通式：c1y^3+c2y^2+c3^y+c4=0;
    AFunc.Coefficients[0] := -1;     //-y^3 + Zc * y^2 + (4Ze-Zb*Ze)y + Zd*Zd + Ze*Zb*Zb - 4*Zc*Ze = 0
    AFunc.Coefficients[1] := Zc;
    AFunc.Coefficients[2] := -Zb*Zd + 4*Ze;
    AFunc.Coefficients[3] := Zd*Zd + Za*Ze-4*Zc*Ze;
    AFunc.Calculate;               //求解配方系数y0
    if AFunc.RootCount > 0 then
      begin         //将解出的y代入(c1^2/4-c2+y)x^2 + (c1y/2-c3)x +(y/2)^2-c4 = 0，以配出完全平方式
        C4 := AFunc.RSolutions[0]/2;     //y0/2    ?????????这里可能需要判断，要取实根。   

        BFunc := T1P2Func.Create(2);    //用配方解配方
        BFunc.Coefficients[0] := Za/4-Zc + AFunc.RSolutions[0];
        BFunc.Coefficients[1] := Zb*C4 - Zd;
        BFunc.Coefficients[2] := C4*C4 - Ze;
        BFunc.Calculate;                //求出配方的常数项x=x0,右边可配成(x-x0)^2
        CFunc := T1P2Func.Create(2);    //消去配方，化为两个对偶的一元二次方程
        DFunc := T1P2Func.Create(2);    // ((x^2+c1x/2）+ y/2)^2 =(x-x0)^2
        CFunc.Coefficients[0] := 1;     //x^2+（c1/2-1）x + y0/2 + x0 = 0

 

跟踪程序运行,常发现一元三次方程AFunc无解返回,导致求解失败.

进而查看一元三次方程求解程序问题.

一元三次方程求解用的是"盛金公式:

 

        盛金公式

　　  一元三次方程aX^3＋bX^2＋cX＋d=0，（a，b，c，d∈R，且a≠0）。
         重根判别式：
         A=b^2－3ac；
         B=bc－9ad；
         C=c^2－3bd，
         总判别式：
         Δ=B^2－4AC。
      

         当A=B=0时，盛金公式①：
         X1=X2=X3=－b/(3a)=－c/b=－3d/c。
     

        当Δ=B^2－4AC>0时，盛金公式②：
        X1=(－b－((Y1)^(1/3)＋(Y2)^(1/3)))/(3a)；
        X2，3=(－2b＋(Y1)^(1/3)＋(Y2)^(1/3)±3^(1/2)((Y1)^(1/3)－(Y2)^(1/3))i)/(6a)，
        其中Y1，2=Ab＋3a(－B±(B^2－4AC)^(1/2))/2，i^2=－1。
     

        当Δ=B^2－4AC=0时，盛金公式③：
        X1=－b/a＋K；
        X2=X3=－K/2，
        其中K=B/A，(A≠0)。
    

        当Δ=B^2－4AC<0时，盛金公式④：
        X1=(－b－2A^(1/2)cos(θ/3))/(3a)；
        X2，3=(－b＋A^(1/2)(cos(θ/3)±3^(1/2)sin(θ/3)))/(3a)，
        其中θ=arccosT，T= (2Ab－3aB)/(2A^(3/2))，(A>0，－1<T<1)。

 

关键原码如下:

{T1P3Func}
Procedure T1P3Func.Calculate;
  var
    Fa,Fb,Fc,Fd : Real;
    DA,DB,DC,DD : Real;
    Y1,Y2,K : Real;
  begin
    Fa := Coefficients[0];
    Fb := Coefficients[1];
    Fc := Coefficients[2];
    Fd := Coefficients[3];
    DA := Fb*Fb - 3*Fa*Fc;
    DB := Fb*Fc - 9*Fa*Fd;
    DC := Fc*Fc - 3*Fb*Fd;
    DD := DB*DB - 4 * DA * DC;
    if (DA = DB) and (DA =0) then
      begin
        RSolutions[0]:= -Fb/(3*Fa);
        RSolutions[1]:= -Fb/(3*Fa);
        RSolutions[2]:= -Fb/(3*Fa);
        VSolutions[0]:= 0;
        VSolutions[1]:= 0;
        VSolutions[2]:= 0;
        FRootCount := 1;
      end else if DD > 0 then
      begin
        Y1:= DA * Fb + 3*Fa*(-Fb+sqrt(DD))/2;
        Y2:= DA * Fb + 3*Fa*(-Fb-sqrt(DD))/2;
        if (Y1>0) and (Y2>0) then
          begin
            Y1:= power(Y1,1/3);
            Y2:= power(Y2,1/3);
          end else if (Y1>0) and (Y2<0) then
          begin
            Y1:= power(Y1,1/3);
            Y2:= -power(-Y2,1/3);
          end else if (Y1<0) and (Y2<0) then
          begin
            Y1:= -power(-Y1,1/3);
            Y2:= -power(-Y2,1/3);
          end else if (Y1<0) and (Y2>0) then
          begin
            Y1:= -power(-Y1,1/3);
            Y2:= power(Y2,1/3);
          end;
          RSolutions[0]:= (-Fb - Y1 - Y2)/(3*Fa);
          VSolutions[0]:= 0;
          RSolutions[1]:= (-2*Fb+Y1+Y2)/(6*Fa);
          VSolutions[1]:= Sqrt(3)*(Y1-Y2)/(6*Fa);
          RSolutions[2]:= RSolutions[1];
          VSolutions[2]:= -VSolutions[1];
          FRootCount := 3;
       end else if DD = 0 then
      begin
        K := DB/DA;
        RSolutions[0]:= -Fb/Fa + K;
        RSolutions[1]:= -K/2;
        RSolutions[2]:= RSolutions[1];
        VSolutions[0]:= 0;
        VSolutions[1]:= 0;
        VSolutions[2]:= 0;
        FRootCount := 3;
      end else if DD<0 then
      begin
        K := (2*DA*Fb-3*Fa*DB)/(2*Power(DA,3/2));
        Y1 := Arccos(K);
        RSolutions[0] := (-Fb-2*Sqrt(DA)*Cos(Y1/3))/(3*Fa);
        RSolutions[1] := (-Fb+Sqrt(DA)*(Cos(Y1/3)+ sqrt(3)*sin(Y1/3)))/(3*Fa);
        RSolutions[2] := (-Fb+Sqrt(DA)*(Cos(Y1/3)- sqrt(3)*sin(Y1/3)))/(3*Fa);
        VSolutions[0]:= 0;
        VSolutions[1]:= 0;
        VSolutions[2]:= 0;
        FRootCount := 3;
      end;
  end;

仔细对照算法说明,并未发现程序实现和算法说明的差异.用实例测试发现:

这段程序在解DD>0的情况时,不能得到正确的解.

如方程(x^2+1)(x+3)=0的解,一看就知道是:i,-i和-3.

展开为:x^3+3x^2+x+3=0

用程序解时:DD=1200>0(和手工验算符合)

得到的解是:

-1.21822935060556+0i
-0.890885324697222+2.13785617558042i
-0.890885324697222-2.13785617558042i

和i,-i,-3相距甚远,难道盛金公式用的有问题?