文章目录
一、一元四次方程求根算法
天珩公式:https://baike.baidu.com/item/%E4%B8%80%E5%85%83%E5%9B%9B%E6%AC%A1%E6%96%B9%E7%A8%8B%E6%B1%82%E6%A0%B9%E5%85%AC%E5%BC%8F/10721996?fr=aladdin
代码如下(示例):
/** * 沈天珩简化四次方求根公式算法。 * 返回根的数组,数组长度: 0 - 4 * @param {*} param0 a,b,c,d 分别为四、三、二、一次项系数;e为常数项 * @returns {Number[]} */
function cal_quartic_ik ([a, b, c, d, e]){
let a_Quadruple = 4 * a;
let D = 3 * pow(b, 2) - 8 * a * c,
E = - pow(b, 3) + a_Quadruple * b * c - 8 * pow(a, 2) * d,
F = 3 * pow(b, 4) + 16 * pow(a, 2) * pow(c, 2) - 16 * a * pow(b, 2) * c + 16 * pow(a, 2) * b * d - 64 * pow(a, 3) * e;
let A = D ** 2 - 3 * F,
B = D * F - 9 * pow(E, 2),
C = F ** 2 - 3 * D * pow(E, 2);
let delta = B ** 2 - 4 * A * C; //总判别式
//四重实根
if (isZero(D) && isZero(E) && isZero(F)) return [-b / a_Quadruple]
//两个实根,其中一个三重实根
if (isZero(A) && isZero(B) && isZero(C) && !isZero(D * E * F)) {
let m = a_Quadruple * D, n = -b * D;
return [(n + 9 * E) / m, (n - 3 * E) / m]
}
//一对二重根,有可能没有根
if (isZero(E) && isZero(F) && !isZero(D)) {
if (D > 0) return [(-b + sqrt(D)) / a_Quadruple, (-b - sqrt(D)) / a_Quadruple]
if (D < 0) return []
}
//四个根
if (!isZero(A * B * C) && isZero(delta)) {
let temp = sign(A * B * E) * sqrt(D - B / A),
sqrt_BA = sqrt(2 * B / A);
//其余两根为不等实根,四个实根
if (A * B > 0)
return [
(-b + temp + sqrt_BA) / a_Quadruple,
(-b + temp - sqrt_BA) / a_Quadruple,
(-b - temp) / a_Quadruple
]
//其余两根为共轭虚根
if (A * B < 0) return [(-b - temp) / a_Quadruple]
}
//两个不等实根和一对共轭虚根
if (delta > 0) {
let z1 = A * D + 3 * ((-B + sqrt(delta)) / 2.0);
let z2 = A * D + 3 * ((-B - sqrt(delta)) / 2.0);
let z = D ** 2 - D * (sign(z1) * pow(abs(z1), One_third) + sign(z2) * pow(abs(z2), One_third)) +
(sign(z1) * pow(abs(z1), One_third) + sign(z2) * pow(abs(z2), One_third)) ** 2 - 3 * A,
m = (-b + sign(E) * sqrt((D + sign(z1) * pow(abs(z1), One_third) + sign(z2) * pow(abs(z2), One_third)) / 3.0)),
n = sqrt((2 * D - sign(z1) * pow(abs(z1), One_third) - sign(z2) * pow(abs(z2), One_third) + 2 * sqrt(z)) / 3.0);
return [(m + n) / a_Quadruple, (m - n) / a_Quadruple]
}
if (delta < 0) {
if (!(D > 0) && (F > 0)) return []
//四个不等实根
if (isZero(E)) {
let sqrt_F_two = 2 * sqrt(F);
return [
(-b + sqrt(D + sqrt_F_two)) / a_Quadruple,
(-b - sqrt(D + sqrt_F_two)) / a_Quadruple,
(-b + sqrt(D - sqrt_F_two)) / a_Quadruple,
(-b - sqrt(D - sqrt_F_two)) / a_Quadruple
]
}
//四个不等实根
let sqrt_A = sqrt(A);
let theta_one_third = acos((3 * B - 2 * A * D) / (2 * A * sqrt_A)) / 3.0,
y1 = (D - 2 * sqrt_A * cos(theta_one_third)) / 3.0,
y2 = (D + sqrt_A * (cos(theta_one_third) + sqrt(3) * sin(theta_one_third))) / 3.0,
y3 = (D + sqrt_A * (cos(theta_one_third) - sqrt(3) * sin(theta_one_third))) / 3.0,
y4 = sign(E) * sqrt(y1);
return [
(-b + y4 + (sqrt(y2) + sqrt(y3))) / a_Quadruple,
(-b + y4 - (sqrt(y2) + sqrt(y3))) / a_Quadruple,
(-b - y4 + (sqrt(y2) - sqrt(y3))) / a_Quadruple,
(-b - y4 - (sqrt(y2) - sqrt(y3))) / a_Quadruple
]
}
return []
}
二、一元三次方程求根算法
盛金公式:https://baike.baidu.com/item/%E7%9B%9B%E9%87%91%E5%85%AC%E5%BC%8F#9
/** * 三次方求根公式算法。 * 返回根的数组,数组长度: 0 - 3 * @param {*} param0 a,b,c 分别为三、二、一次项系数;d为常数项 * @returns {Number[]} */
function cal_cubic_ik([a, b, c, d]){
let a_Threefold = 3 * a;
let A = b ** 2 - a_Threefold * c,
B = b * c - 9 * a * d,
C = c ** 2 - 3 * b * d,
delta = B ** 2 - 4 * A * C; //总判别式
// 一个三重实根
if (isZero(A) && isZero(B)){
return [-b / (a_Threefold)]
}
// 一个实根和一对共轭复根
if (delta > 0){
let y1 = A * b + a_Threefold * (( -B + sqrt(delta)) / 2.0),
y2 = A * b + a_Threefold * (( -B - sqrt(delta)) / 2.0);
return [(-b - (sign(y1) * pow(abs(y1), One_third) + sign(y2) * pow(abs(y2), One_third)) ) / a_Threefold ]
}
//三个实根,其中有一个二重根
if (isZero(delta)){
let K = B / float(A);
return [- b / float(a) + K, - K / 2.0]
}
//三个不等实根
if (delta < 0){
let theta = acos((2 * A * b - 3 * a * B) / (2 * sqrt(pow(A, 3))));
m = sqrt(A),
n = cos(theta / 3.0),
q = sqrt(3) * sin(theta / 3.0);
return [
(- b - 2 * m * n) / a_Threefold,
(- b + m * (n + q)) / a_Threefold,
(- b + m * (n - q)) / a_Threefold
]
}
}
总结
提示:以上算法均返回实数解,且去掉了重根。
需导入的头文件
const { pow, sqrt, sign, abs, sin, cos, acos} = Math;
const One_third = 1.0 / 3.0;
function isZero (value, ep = EP) {
return Math.abs(value) < ep;
}