一元高次方程的求解

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

 

1、一次方程

a1x+a0=0

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

x=-a0/a1

 

2、二次方程

a2x2+a1x+a0=0

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

 

判别式Δ=a1-4a2a0

令p1=-0.5*a1/a2,p2=0.5*SQRT(|Δ|)

当Δ>=0时,x1=p1+SGN(p1)*p2,若|x1|>0,x2=a0/(a2*x1),否则,x2=p1-SGN(p1)*p2

 

3、三次方程

a3x3+a2x2+a1x+a0=0

令y=3a3x+a2,方程可变形为y3+py+q=0

p=9a3a1-3a22

q=2a23-9a3a2a1+27a32a0

 

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

首先假定方程y3+py+q=0的根的形式为y=A1/3+B1/3,也即两数开立方之和。等式两边立方,整理可得

y3-3(AB)1/3y-(A+B)=0,与原方程比较则有

A+B=-q

AB=(-q/3)3

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

解次方程可得

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

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

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

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

 

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

 

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

 

为了避免进行除法计算,判别式变形如下Δ=3(27q2+4p3)

当Δ<=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)2=(-p/3)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)3/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-a2)/(3a3)

 

4、四次方程

a4x4+a3x3+a2x2+a1x+a0=0

令y=4a4x+a3,方程可变形为y4+py2+qy+r=0

p=-6a32+16a2a4

q=8a32-32a2a3a4+64a1a42

r=-3a34+16a2a32a4-64a1a3a42+256a43a0

方程两侧加上2ty2+t2,整理得(y2+t)2=(2t-p)y2-qy+t2-r

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

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

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

 

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

 

关于t的一元三次方程b3t3+b2t2+b1t+b0=0的系数

b3= -8

b2= 4p

b1= 8r

b0= q2-4pr

关于y的一元二次方程c2y2+c1y+c0=0的系数

 

第2组 c2= 1,c1= -SQRT(2t-p),c0= t+0.5q/SQRT(2t-p)

 

第1组 c2= 1,c1= SQRT(2t-p),c0= t-0.5q/SQRT(2t-p)

 

 

xi=(yi-a3)/(4a4)

 

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

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

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

  

 

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

 

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

 

劈因子法:

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

 

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

假设Q(x)=bn-2xn-2+bn-3xn-3+...+b1x+b0,并取bn=bn-1=0,可得递推公式bi=ai+2-ubi+1-vbi+2,i=n-2, n-1, ...,0。

 

在等式Pn(x)=U(x)Q(x)+r(x)中提取关于x的一次项的系数和常数项 a1x+a0=(ux+v)(b1x+b0)+(r1x+r0),可得r1=a1-ub0-vb1,r0=a0-vb0

设cui=[dbi/ du],cvi=[dbi/ dv],可知

cun=0,cvn=0

cun-1,cvn-1=0

cun-2=0,cvn-2=0

当i=n-3, n-4, ..., 1, 0时,cui=-bi+1-u*cui+1-v*cui+2,cvi=-u*cvi+1-bi+2-v*cvi+2

 

[∂r1/∂u] = -b0-u*cu0-v*cu1

[∂r1/∂v] = -u*cv0-b1-v*cv1

[∂r0/∂u] = -v*cu0

[∂r0/∂v] = -b0-v*cv0

 

这样就可以求解关于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;



补充:

QR方法求多项式方程的根
    设多项式方程为
        pn(x)=x^n+a1*x^(n-1)+......+an-1*x+an=0        (1)
构造一个方阵如下
    -a1   -a2     ... -(an-1)  -an
     1     0      ...    0      0
     0     1      ...    0      0
     .     .             .      .   
     .     .             .      .
     .     .             .      .
     0     0      ...    0      0
     0     0      ...    1      0
用QR方法求出该方阵的所有特征值(正好n个,包括重复的),它们即为方程(1)的n个根。
    QR方法比较麻烦,你可以从数值代数类的书中找到该算法的具体实现。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值