一元高次方程的求解

    最近在写程序的时候遇到一些公式需要求解一元三次、一元四次方程，在google上搜索到不少公式和公式后面的故事。以下将逐步讨论各阶次的一元方程的解法及程序上的实现。

 

1、一次方程

a₁x+a₀=0

一元一次方程实际上也可以叫做一元线性方程，这个求解很简单，编程上也不需做特出处理。 

x=-a₀/a₁

 

2、二次方程

a₂x²+a₁x+a₀=0

一元二次方程的求根公式推导很简单，这里不赘述。编程的时候要注意处理一种比较特殊的情况。例如求解方程x²-20000x+1=0的时候，如果直接使用求根公式编写程序的话，求得跟为20000和0，而实际的根约为19999.99995和0.00005。使用下面的办法可以避免上面这种方程的较小根为0的情况。

 

判别式Δ=a₁-4a₂a₀

令p1=-0.5*a₁/a₂，p2=0.5*SQRT(|Δ|)

当Δ>=0时，x1=p1+SGN(p1)*p2，若|x1|>0，x2=a₀/(a₂*x1)，否则，x2=p1-SGN(p1)*p2

 

3、三次方程

a₃x³+a₂x²+a₁x+a₀=0

令y=3a₃x+a₂，方程可变形为y³+py+q=0

p=9a₃a₁-3a₂²

q=2a₂³-9a₃a₂a₁+27a₃²a₀

 

关于方程y³+py+q=0的求解过程主要参考了百度百科关于“解一元三次方程”的条目，下面简单复述一下。

首先假定方程y³+py+q=0的根的形式为y=A^1/3+B^1/3，也即两数开立方之和。等式两边立方，整理可得

y³－3(AB)^1/3y－(A+B)＝0，与原方程比较则有

A+B=-q

AB=(-q/3)³

这是可以根据韦达定理建立一元二次方程W²-(A+B)W+AB=0

解次方程可得

A=W1=-(q/2)-SQRT((q/2)²+(p/3)³)

B=W2=-(q/2)+SQRT((q/2)²+(p/3)³)

这样方程的一个根y1=[-(q/2)-SQRT((q/2)²+(p/3)³)]^1/3+[-(q/2)+SQRT((q/2)²+(p/3)³)]^1/3

再根据韦达定理y1+y2+y3=0和y1y2y3=-q，可求得y2+y3和y2y3，再次使用韦达定理建立一元二次方程即可求解。

 

仔细观察y1的表达式会发现存在一个判别式的问题，也即Δ=((q/2)²+(p/3)³是否大于零，但是比较巧妙的是，无论Δ是否大于零，y1始终为实数根，Δ的符号对其他两个根有影响，Δ>0时，y2、y3为共轭的虚数根，Δ<=0是，y2、y3为实数根，Δ=0时，y2=y3。

 

三次方程同样存在方程的病态的问题，也需要使用韦达定理来消除相近大数相减造成较小的根为0的情况。编写程序的时候可按下面的步骤操作。

 

为了避免进行除法计算，判别式变形如下Δ=3(27q²+4p³)

当Δ<=0时， 方程有三个实数根，此时p<=0

令R=2*SQRT(-p/3)，θ=arccos(-0.5q(-p/3)^-1.5)/3

y1=Rcos(θ)

y2=Rcos(θ+2π/3)

y3=Rcos(θ-2π/3)

 

当Δ=0时，(q/2)²=(-p/3)³

若q<0，则θ=0，y1=R,y2=y3=-0.5R

若q>0，则θ=π/3，有y2=-R，y1=y3=0.5R

 

当Δ>0时，方程只有一个实数根，也即上面推导出的公式

令R1=(-108q+12SQRT(Δ))^1/3，R2=(-108q-12SQRT(Δ))^1/3

R1和R2表达式里面根号里面的部分是一个二次方程的根，所以也存在处理病态条件的问题

令K1=-108q-SGN(q)*12SQRT(Δ)，K2=(-12p)³/K1

R1=(K1)^1/3，R2=(K2)^1/3

y1=(R1+R2)/6

y2=-(R1+R2)/12+ i (R1-R2)p/6

y2=-(R1+R2)/12- i (R1-R2)p/6

 

xi=(yi-a₂)/(3a₃)

 

4、四次方程

a₄x⁴+a₃x³+a₂x²+a₁x+a₀=0

令y=4a₄x+a₃，方程可变形为y⁴+py²+qy+r=0

p=-6a₃²+16a₂a₄

q=8a₃²-32a₂a₃a₄+64a₁a₄²

r=-3a₃⁴+16a₂a₃²a₄-64a₁a₃a₄²+256a₄³a₀

方程两侧加上2ty²+t²，整理得(y²+t)²=(2t-p)y²-qy+t²-r

方程左侧为完全平方式，右侧也应为完全平方式，则可令右侧多项式的Δ=0

也即(-q)2-4(2t-p)(t²-r)=0，这是关于t的一元三次方程，可按上述办法求解，得到满足t>=0.5p的实数根(避免求解虚系数的方程)

进而方程转化为关于y的两个一元二次方程，继续求解可得到方程的4个根

 

四次方程的求解是建立在三次和二次方程的基础上的，虽然同样存在方程的病态的问题，但并不需要额外处理。编写程序的时候可按下面的步骤操作。

 

关于t的一元三次方程b₃t³+b₂t²+b₁t+b₀=0的系数

b₃= -8

b₂= 4p

b₁= 8r

b₀= q²-4pr

关于y的一元二次方程c₂y²+c₁y+c₀=0的系数

 

第2组 c₂= 1，c₁= -SQRT(2t-p)，c₀= t+0.5q/SQRT(2t-p)

 

第1组 c₂= 1，c₁= SQRT(2t-p)，c₀= t-0.5q/SQRT(2t-p)

 

 

xi=(yi-a₃)/(4a₄)

 

-----------------------------------------------------------------------------

    以上讨论的都是使用方程的解析解来求解高次方程的根的办法。对三次及更高次的方程，可以使用通用的方式来求解方程的根。求解方程的目的就是获取高精度的方程的根，数值计算的办法虽然需要多次迭代，但是通用性较好，精度也比较高。

-----------------------------------------------------------------------------

  

 

5、劈因子法求解的高次方程

 

    “劈因子法”参考了“互动百科”中“高次代数方程求根”这个条目。“CSDN博客”上我也找到了对应的文章有完整的程序来实现这个算法，参见“一元高次方程求解”。这里对“劈因子法”简单介绍一下，然后推导一下计算出方程两个根之后，其余的根的求解方法。

 

劈因子法：

    多项式的方程 P_n(x)=a_nxⁿ+a_n-1x^n-1+...+a₁x+a₀=0，称为n次代数方程，多项式的零点就是对应代数方程的根。
 　用x的二次式 U(x)=x²+ux+v 除P_n(x)则得商Q(x)及余式 r(x)=r₁(x)+r₀，因而有 P_n(x)=U(x)Q(x)+r(x).........(1)
 　设U(x)是一个近似的二次因式，问题是怎样修改u和v使对应的余式更接近于零。为此，作线性近似，取U(x)=x²+(u+du)*x+(v+dv)，是r(x)=(r₁+dr₁)x+(r₀+dr₀)=0，则修正量du、dv应满足方程组 { dr₁+r₁=0，dr₀+r₀=0 }，
进一步可以写成 { [∂r₁/∂u]du + [∂r₁/ ∂v]dv + r₁=0，[∂r₀/∂v]du + [∂r₀/ ∂v]dv + r₀=0 } .........(2)
    利用已知关系可求出 [∂r₁/∂u]、[∂r₁/ ∂v]、[∂r₀/∂u]、[∂r₀/ ∂v] 带入(2)后，就能求出u和v的校正量du和dv。而u+du、v+dv就是更好的二次因式的两个系数。

 

下面将公式中未明确的项加以推导：

假设Q(x)=b_n-2x^n-2+b_n-3x^n-3+...+b₁x+b₀，并取b_n=b_n-1=0，可得递推公式b_i=a_i+2-ub_i+1-vb_i+2，i=n-2, n-1, ...,0。

 

在等式P_n(x)=U(x)Q(x)+r(x)中提取关于x的一次项的系数和常数项 a₁x+a₀=(ux+v)(b₁x+b₀)+(r₁x+r₀)，可得r₁=a₁-ub₀-vb₁，r₀=a₀-vb₀。

设cu_i=[db_i/ du]，cv_i=[db_i/ dv]，可知

cu_n=0，cv_n=0

cu_n-1，cv_n-1=0

cu_n-2=0，cv_n-2=0

当i=n-3, n-4, ..., 1, 0时，cu_i=-b_i+1-u*cu_i+1-v*cu_i+2，cv_i=-u*cv_i+1-b_i+2-v*cv_i+2。

 

[∂r₁/∂u] = -b₀-u*cu₀-v*cu₁

[∂r₁/∂v] = -u*cv₀-b₁-v*cv₁

[∂r₀/∂u] = -v*cu₀

[∂r₀/∂v] = -b₀-v*cv₀

 

这样就可以求解关于du和dv的2阶线性方程组，解出du和dv

反复进行上述计算，知道du和dv趋于0。

求解U(x)=0就可以得到方程的两个根x1和x2。

继续按上述办法求解Q(x)=0的到余下的根。

 

 

#include "stdafx.h"
#include <afxtempl.h>
#include <math.h>
typedef CArray< long double, long double > CLongDoubleArray;

int PolynomialEquation(CLongDoubleArray & A, CLongDoubleArray & X)
{
 int r = 0;

 int n = A.GetSize() - 1;//方程的最高次数

 if (n < 1) return r;
 CLongDoubleArray a, x;
 a.SetSize(n+1);
 TRACE(_T("Pn(x)="));
 for (i = n; i >= 0; i--)
 {
  a[i] = static_cast<long double>(A[i]);
  TRACE(_T("(%f)*x^%d%s"), a[i], i, i==0?_T("/r/n"):_T("+"));
 }
 x.SetSize(2 * n);
 int m = (n + (n % 2)) / 2;//需要计算的次数
 long double u = 0.0;//U(x)的一次项
 long double v = 0.0;//U(X)的常数项
 long double du = 0.0;
 long double dv = 0.0;
 int count = 0;
 
 CLongDoubleArray b, cu, cv;
 b.SetSize(n + 1);//Q(x)的系数
 cu.SetSize(n + 1);//偏导数u
 cv.SetSize(n + 1);//偏导数v

 long double r0, r1;
 long double pr0u, pr0v;
 long double pr1u, pr1v;

 long double det, p1, p2;
 for (k = 0; k < m; k++)
 {
  u = 4.0;
  v = 3.0;//这里固定了取值，对于一般计算满足要求
  if (n >= 3)//三次及三次以上的方程
  {
   count = 0;
   b[n] = b[n-1] = 0.0;
   cu[n] = cu[n-1] = cu[n-2] = 0.0;
   cv[n] = cv[n-1] = cv[n-2] = 0.0;
   do
   {
    TRACE(_T("-----/r/nk = %d, count = %d/r/n-----/r/n"), k, count);
    TRACE(_T("Q(x)="));
    for (i = n - 2; i >= 0; i--)
    {
     b[i] = a[i+2] - u * b[i+1] - v * b[i+2];
     TRACE(_T("(%f)*x^%d%s"), b[i], i, i==0?_T("/r/n"):_T("+"));
    }
    r1 = a[1] - u * b[0] - v * b[1];
    r0 = a[0] - v * b[0];
    TRACE(_T("r(x) = (%.10f)*x^1+(%.10f)/r/n"), r1, r0);
    
    for (i = n - 3; i >= 0; i--)
    {
     cu[i] = - b[i+1] - u * cu[i+1] - v * cu[i+2];
     cv[i] = - u * cv[i+1] - b[i+2] - v * cv[i+2];
    }
    pr1u = - b[0] - u * cu[0] - v * cu[1];
    pr1v = - u * cv[0] - b[1] - v * cv[1];
    pr0u = - v * cu[0];
    pr0v = - b[0] - v * cv[0];
    
    du = (r1 * pr0v - r0 * pr1v) / (pr1v * pr0u - pr1u * pr0v);
    dv = (r0 * pr1u - r1 * pr0u) / (pr1v * pr0u - pr1u * pr0v);
    
    u += du;
    v += dv;

    count++;
    if (count > 1000) break;//这里可以适当调整,u、v取值适当的话，收敛很快
   } while (fabs(du) > 1.0e-10 || fabs(dv) > 1.0e-10);
   //确定下一步计算的数据
   n -= 2;
   a.SetSize(n+1);
   for (i = n; i >= 0; i--) a[i] = b[i];
  }
  else
  {
   if (n == 2)//剩下一个二次式
   {
    if (fabs(a[2]) > 0.0)
    {
     u = a[1] / a[2];
     v = a[0] / a[2];
    }
    else
    {
     //不应该出现这种情况
    }
   }else
   if (n == 1)//剩下一个一次式
   {
    r += 11;
    if (fabs(a[1]) > 0.0)
    {
     x[4 * k] = - a[0] / a[1];
    }
    else
    {
     //不应该出现这种情况
    }
    break;
   }
  }

  //下面求解一元二次方程
  det = u * u - 4.0 * v;
  p1 = - 0.5 * u;
  p2 = 0.5 * sqrt(fabs(det));
  if (det >= 0.0)
  {
   r += 22;
   x[4 * k] = p1 + (p1 > 0.0 ? 1.0 : p1 < 0.0 ? -1.0 : 0.0) * p2;
   if (fabs(X[4 * k]) > 0.0)
   {
    x[4 * k + 2] = v / x[4 * k];
   }
   else
   {
    x[4 * k + 2] = p1 - (p1 > 0.0 ? 1.0 : p1 < 0.0 ? -1.0 : 0.0) * p2;
   }
  }
  else
  {
   r += 20;
   x[4 * k    ] =  p1;
   x[4 * k + 1] =  p2;
   x[4 * k + 2] =  p1;
   x[4 * k + 3] = -p2;
  }
 }
 
 X.SetSize(x.GetSize());
 for (i = 0; i < x.GetSize(); i++) X[i] = x[i];
 return r;

}