一元高次方程的求解
2014-08-29 03:32
337 查看
最近在写程序的时候遇到一些公式需要求解一元三次、一元四次方程，在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
= b[n-1] = 0.0;
cu
= cu[n-1] = cu[n-2] = 0.0;
cv
= 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方法比较麻烦，你可以从数值代数类的书中找到该算法的具体实现。
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
= b[n-1] = 0.0;
cu
= cu[n-1] = cu[n-2] = 0.0;
cv
= 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方法比较麻烦，你可以从数值代数类的书中找到该算法的具体实现。
内容来自用户分享和网络整理，不保证内容的准确性，如有侵权内容，可联系管理员处理 
相关文章推荐
- 一元两次方程求解
- 一元方程求解
- nyoj 1178 && hdu 5105 Math Problem 求解一元二次三次方程
- EXCEL 单变量求解 解一元三次方程
- 使用二分法求解一元N次方程的近似值
- 二分逼近/牛顿迭代——一元高次非线性方程求解
- 用类方法求解一元二次实系数方程
- TF随笔-14-二分法求解一元方程
- 求解一元多次方程(迭代法)
- Matlab练习——多项式和一元方程求解
- 高斯消元法求解一元多次方程 C/C++实现
- 问题六十五：二叉查找树的一个应用实例——求解一元十次方程时单实根区间的划分
- 求解一元多次方程的两种方法：牛顿迭代法和二分法
- 一元高次方程的求解
- 求解一元多次方程 牛顿迭代法
- 弦截法求一元三次方程的近似解
- 求一元二次方程，并且同时判断解的个数
- 编程求一元二次方程。ax²+bx+c=0（a≠0）
- C#一元二次方程的程序；
- POJ2115——C Looooops(扩展欧几里德+求解模线性方程)