一元三次方程求解(求根) - 盛金公式法
一、引言
只含有一个未知数（即“元”），并且未知数的最高次数为3（即“次”）的整式方程叫做一元三次方程（英文名：cubic equation in one unknown）。一元三次方程的标准形式（即所有一元三次方程经整理都能得到的形式）是ax3+bx2+cx+d=0（a，b，c，d为常数，x为未知数，且a≠0）。一元三次方程的公式解法有卡尔丹公式法。(--引自百度百科)；
由于卡尔丹公式法和多项式求解的方法不容易理解，不容易用程序实现，经过文献的查找了解到盛金公式求解一元三次方程，经验证求根正确可以引用。
二、盛金公式简介
80年代，中国中学数学教师范盛金对解一元三次方程问题进行了深入的研究和探索，发明了比卡尔丹公式更适用的新根公式 - 盛金公式，并建立了简明的、直观的、使用的新判别法 - 盛金判别法，同时提出了盛金定理；盛金定理清晰地回答了解一元三次方程的疑惑问题，并且很有趣味。
盛金公式的特点是由最简重根判别式：, , ；和总判别式来构成，体现了数学的有序、对称、和谐与简洁美。
三、盛金公式表达式
一元三次方程式为：
(1)
重根判别式为：
(2)
总判别式为：
(3)
盛金判别法的结论为：
1. 条件一： 当A = B = 0时，方程有一个三重实根
(4)
2. 条件二： 当 时 ，方程有一个实根和一个共轭虚根
(5)
其中:
;
3. 条件三：当时，方程有三个实根，其中有一个两重根
(6)
其中：
;
4. 条件四： 当时，方程有三个不相等的实根
(7)
其中：
四、盛金公式求解程序流程
五、java和matlab实例代码
matlab代码
function x = solve3Polynomial(a, b, c, d)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% funtion x = solve3Polynomial(a, b, c, d)
%
%> @brief 利用盛金公式求解三阶多项式的解.
%>
%> @details 输入三阶多项式系数，求解ax^3 + bx^2 + cx^1 + d = 0的根
%> 参考文献：范盛金. 一元三次方程的新求根公式与新判别法[J]. 海南师范学院学报,1989,2(2):91-98.
%>
%> @param[out] x 求解完成的三个根x1，x2，x3
%> @param[in] a 三次项系数
%> @param[in] b 二次项系数
%> @param[in] c 一次项系数
%> @param[in] d 零次项系数
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 重根判别式
A = b*b - 3*a*c; if abs(A) < 1e-14; A = 0; end
B = b*c - 9*a*d; if abs(B) < 1e-14; B = 0; end
C = c*c - 3*b*d; if abs(C) < 1e-14; C = 0; end
%% 总判别式
DET = B*B - 4*A*C; if abs(DET) < 1e-14; DET = 0; end
%% 条件一
if (A == 0) && (B == 0)
x1 = -c/b; x2 = x1 ; x3 = x1;
end
%% 条件二
if DET > 0
Y1 = A*b + 1.5*a*(-B + sqrt(DET));
Y2 = A*b + 1.5*a*(-B - sqrt(DET));
y1 = nthroot(Y1,3); y2 = nthroot(Y2,3);
x1 = (-b-y1-y2)/(3*a);
vec1 = (-b + 0.5*(y1 + y2))/(3*a);
vec2 = 0.5*sqrt(3)*(y1 - y2)/(3*a);
x2 = complex(vec1, vec2);
x3 = complex(vec1, -vec2);
clear Y1 Y2 y1 y2 vec1 vec2;
end
%% 条件三
if DET == 0 && (A ~= 0) && (B ~= 0)
K = (b*c-9*a*d)/(b*b - 3*a*c); K = round(K,14);
x1 = -b/a + K; x2 = -0.5*K; x3 = x2;
end
%% 条件四
if DET < 0
sqA = sqrt(A);
T = (A*b - 1.5*a*B)/(A*sqA);
theta = acos(T);
csth = cos(theta/3);
sn3th = sqrt(3)*sin(theta/3);
x1 = (-b - 2*sqA*csth)/(3*a);
x2 = (-b + sqA*(csth + sn3th))/(3*a);
x3 = (-b + sqA*(csth - sn3th))/(3*a);
clear sqA T theta csth sn3th;
end
x = [x1; x2; x3];
end
java代码
public double solve3Polynomial(double a, double b, double c, double d)
{
double x1 = 0.0, x2 = 0.0, x3 = 0.0;
double A = 0.0, B = 0.0 , C =0.0 , DET = 0.0;
// 1. 计算重根判别式
A = b*b - 3*a*c;
if(Math.abs(A) < 1e-14)
{
A = 0.0;
}
B = b*c - 9*a*d;
if(Math.abs(B) < 1e-14)
{
B = 0.0;
}
C = c*c - 3*b*d;
if(Math.abs(C) < 1e-14)
{
C = 0.0;
}
// 2. 计算总判别式
DET = B*B - 4*A*C;
if(Math.abs(DET) < 1e-14)
{
DET = 0;
}
// 3. 条件一，计算根
if((A == 0) && (B == 0))
{
x1 = (-1*c)/b;
x2 = x1;
x3 = x1;
//Log.i("roots", "条件一:" + x1 );
}
// 4. 条件二，计算根
if(DET > 0)
{
double Y1 = A*b + 1.5*a*(-1*B + Math.sqrt(DET));
double Y2 = A*b + 1.5*a*(-1*B - Math.sqrt(DET));
//Log.i("SQSddd", "Y1: " + Y1 + " Y2: " + Y2);
double y1 = getCubeRoot(Y1);
double y2 = getCubeRoot(Y2);
x1 = (-1.0*b-(y1+y2))/(3.0*a); // 一个实根
double vec1 = (-1*b + 0.5*(y1 + y2))/(3*a);
double vec2 = 0.5*Math.sqrt(y1 - y2)/(3*a);
// x2 = Math.
double x3_real = 0.0, x2_real = 0.0; // 实部
x2_real = (-b + getCubeRoot(Y1)) / (3 * a);
x3_real = x2_real;
double x2_virtual = 0.0, x3_virtual = 0.0; // 虚部
x2_virtual = ((Math.sqrt(3) / 2) * (y1 - y2 )) / (3 * a);
x3_virtual = -x2_virtual;
//Log.i("roots", "条件二:" + x1 );
}
// 5. 条件三，计算根
if(DET == 0 && A != 0 && B != 0)
{
double K = (b*c - 9*a*d)/(b*b -3*a*c);
K = Math.round(K);
x1 = (-1.0*b)/a + K;
x2 = -1*0.5*K;
x3 = x2;
//Log.i("roots", "条件三:" + x1 );
}
// 6.条件四，计算根
if(DET < 0)
{
double sqA = Math.sqrt(A);
double T = (A*b - 1.5*a*B)/(A*sqA);
double theta = Math.acos(T);
double csth = Math.cos(theta/3);
double sn3th = Math.sqrt(3)*Math.sin(theta/3);
x1 = (-1*b - 2*sqA*csth)/(3*a);
x2 = (-1*b + sqA*(csth + sn3th))/(3*a);
x3 = (-1*b + sqA*(csth - sn3th))/(3*a);
//Log.i("SQSd", "条件四:" + x1 );
}
// 7. 返回计算结果
return x1;
}
public double getCubeRoot(double value) {
if (value < 0) {
return -Math.pow(-value, 1.0/3.0);
} else if (value == 0) {
return 0;
} else {
return Math.pow(value, 1.0/3.0);
}
}