equation.pas

一个解一元多次方程的单元

unit Equation;

interface

type
  TValues = array of double;
  function Get3EquationValue(a, b, c: Double): TValues;
  function Get4EquationValue(a, b, c, d: Double): TValues;
implementation

uses math;

const
  Invert3 = 1 / 3;


//求一元三次方程的根，形如 X^3 + a * (X^2) + b * X + c = 0
function Get3EquationValue(a, b, c: Double): TValues;
var
  i: Integer;
  p, q, y, z, Judge: Double;
  RealPart, VirtualPart, Mould: Double;   //虚根的幅角的Cos,Sin，实部，虚部和模
  AngA: Double;
begin
  //将一元三次方程转换成形如 X^3 + p * X + q = 0, 令 X = Y - a / 3
  p := b - a * a / 3;
  q := (a * (2 * Sqr(a) - 9 * b)) / 27 + c;
  //判断一元三次方程有几个根的判别式，公式为 △ = q^2 / 4 + p^3 / 27;
  // △ > 0, 有一个实根和一对共轭复根；△ = 0，都是实根，其中有两个相等；
  //△ < 0, 有三个相异的实根, △ = Judge := Sqr(q) / 4 + p * p * p / 27;
  //y := Power(-q/2+sqrt(Judge), 1/3)
  Judge := Power(q / 2, 2) + Power(p / 3, 3);
  if Judge > 1E-10 then
  begin
    SetLength(Result, 1);
    y := -q / 2 + Sqrt(Judge);
    y := Sign(y) * Power(Abs(y), Invert3);
    z := -p / (3 * y);
    Result[0] := y + z;
  end
  else if abs(Judge) <= 1E-10 then
  begin
    SetLength(Result, 2);
    y := -q / 2;
    Result[0] := Sign(y) * 2 * Power(Abs(y), Invert3);
    Result[1] := -1;
  end else
  begin
    RealPart := -q / 2;
    VirtualPart := Sqrt(Abs(Judge));
    Mould := Sqrt(Sqr(RealPart) + Abs(Judge));
//    SinA := RealPart / Mould;
//    CosA := VirtualPart / Mould;
//    AngA := ArcTan2(CosA, SinA); //为什么不直接使用AngA := ArcTag(VirtualPart / RealPart)
    AngA := ArcTan2(VirtualPart, RealPart);
    AngA := AngA / 3;         //z := r*(cos(Ang) + i * sin(Ang)), z^n = r^n(cos(n*Ang) + i*sin(n*Ang))
    Mould := Power(Mould, Invert3);
    RealPart := Mould * Cos(AngA);
    VirtualPart := Mould * Sin(AngA);
    SetLength(Result, 3);
    Result[0] := 2 * RealPart;
    Result[1] := -RealPart - VirtualPart * Sqrt(3);
    Result[2] := -RealPart + VirtualPart * Sqrt(3);
  end;
  for i := 0 to High(Result) do
    Result[i] := Result[i] - a / 3;
end;

//求一元四次方程的根，形如 X^4 + a * (X^3) + b * (X^2) + c * X + d = 0
function Get4EquationValue(a, b, c, d: Double): TValues;
var
  i: Integer;
  Values: TValues;
  z, n1, n2: Double;
  HasRst: Boolean;
  Judge, p, q: Double;
  Len: Integer;
begin
  //转换成一元三次标准方程，结果为 z^3 - b * z^2 + (a*c-4*d)*z + d*(4b-a^2)-c^2=0
  Values := Get3EquationValue(-b, a * c - 4 * d, d * (4 * b - a * a) - c * c);
  HasRst := False;                              
  //计算出立方根后，取其中一个，计算N1,2=±Sqrt(4*z1+a^2-4*b)
  n1 := 0; z := 0;
  for i := 0 to High(Values) do
  begin
    z := Values[i];
    n1 := 4 * z + a * a - 4 * b;
    if n1 > 0 then
    begin
      n1 := Sqrt(n1);
      HasRst := True;
      Break;
    end
    else if abs(n1) <= 1E-10 then
    begin
      SetLength(Result, 1);
      Result[0] := 0;
      Exit;
    end;
  end;
  if not HasRst then
    Exit;
  n2 := -n1;
  //求出n1,n2后,解二次方程X^2+ (a+n1,2) * X/2 + 0.5*((z1 + (az1-2c)/n1,2) = 0
  //△ = Judge = b^2 - 4ac
  p := 0.5 * (a + n1);
  q := 0.5 * (z + (a * z - 2 * c) / n1);
  Judge := Sqr(p) - 4 * q;
  if Judge > 1E-10 then
  begin
    SetLength(Result, 2);
    Result[0] := 0.5 * (-p + Sqrt(Judge));
    Result[1] := 0.5 * (-p - Sqrt(Judge));
  end
  else if abs(Judge) < 1E-10 then
  begin
    SetLength(Result, 1);
    Result[0] := -0.5 * p;
  end;

  p := 0.5 * (a + n2);
  q := 0.5 * (z + (a * z - 2 * c) / n2);
  Judge := Sqr(p) - 4 * q;
  Len := Length(Result);
  if Judge > 1E-10 then
  begin
    SetLength(Result, Len + 2);
    Result[Len] := 0.5 * (-p + Sqrt(Judge));
    Result[Len + 1] := 0.5 * (-p - Sqrt(Judge));
  end
  else if abs(Judge) < 1E-10 then
  begin
    SetLength(Result, Len + 1);
    Result[Len] := -0.5 * p;
  end;
end;


end.