3次方程式の実数解を求める.js

複素数型がなくてもなんとかなるのよ

方針

実数解3つの場合は「ビエトの方法」

3重解持ちの場合は微分した2次方程式を解く

重解+実数解の場合は微分した2次方程式の解のどちらかが重解、もう一個の解は解の公式から求めれる

実数解1つ+共役複素数の場合は「カルダノの方法」

cf.

code:script.js

/**

* 3次以下の方程式の実数解をすべて求める

* @param {number} a - 3次の係数

* @param {number} b - 2次の係数

* @param {number} c - 1次の係数

* @param {number} d - 定数項

* @returns {number[]}

*/

const solveAtMostCubic = (a, b, c, d) => {

/** @description - a = 0 の場合は2次方程式になるので、それを解く */

if (a === 0) {

/** @description - b = 0 の場合は1次方程式になるので、それを解く */

if (b === 0) {

/** @description - 1次方程式 cx + d = 0 の実数解は x = -d / c */

}

/** @description - 2次方程式 bx^2 + cx + d = 0 の実数解は x = (-c ± √(c^2 - 4bd)) / 2b */

// 判別式 : D = c² - 4bd → D > 0 なら2つの実数解、D = 0 なら実数重解、D < 0 なら実数解なし

const D = c ** 2 - 4 * b * d;

if (D === 0) {

} else if (D > 0) {

} else {

return [];

}

}

/** @description - 3次方程式 ax^3 + bx^2 + cx + d = 0 → 解の様子に応じて手法を分岐 */

// 判別式(1) : D₁ = b²c² - 4ac³ - 4b³d - 27a²d² + 18abcd

const D1 = b ** 2 * c ** 2 - 4 * a * c ** 3 - 4 * b ** 3 * d - 27 * a ** 2 * d ** 2 + 18 * a * b * c * d;

if (D1 === 0) {

/** @description - パターン①② : D₁ = 0 → 少なくとも1つの実数重解を持つ */

// 重解は微分した 3ax² + 2bx + c = 0 の解にもなるので、一旦そっちを求めておく

const differentiated = solveAtMostCubic(0, 3 * a, 2 * b, c);

// 判別式(2) : D₂ = -2b³ + 9abc - 27a²d

const D2 = -2 * b ** 3 + 9 * a * b * c - 27 * a ** 2 * d;

if (D2 === 0) {

/** @description - パターン① : D₂ = 0 → 実数の三重解を持つ */

// このとき、さっき求めた"微分した方程式"は重解になっているのでそのまま返せばOK

return differentiated;

} else {

/** @description - パターン② : D₂ ≠ 0 → 1つの実数重解と1つの実数解を持つ */

// "微分した方程式"は重解になっていない。どっちかが元の3次方程式の重解なので両方試す

const multiRoot = differentiated.find(x_ => {

// 実際に計算してみる

const left_side = a * x_ ** 3 + b * x_ ** 2 + c * x_ + d;

// たぶん誤差が出るので、ある程度0に近ければ0とみなして正解とする

const epsilon = Math.abs(a * Number.EPSILON * 8);

return Math.abs(left_side) < epsilon;

});

if (multiRoot === undefined) throw new Error("???????");

// もう一つの解は「-(b/3a) + 2׳√(D₂/54a)」になる

const anotherRoot = -(b / (3 * a)) + 2 * Math.cbrt(D2 / (54 * a));

}

}

/** @description パターン③④ : D₁≠0 → どちらにせよ立体完成が必要なのでやる */

// 各係数をaで割って x₃ + Bx² + Cx + D = 0 の形にする

// x = y - B/3 として変数変換をすると二次の項が消えて y³ + py + q = 0 になる

// このときの p, q の値を求める

const p = C - B ** 2 / 3;

const q = D - (C * B / 3) + (2 * B ** 3 / 27);

// 分岐内で求めたyを入れる配列を先に用意

const y_ = [];

if (D1 > 0) {

/** @description - パターン③ : D₁ > 0 → 相異なる3つの実数解を持つ */

// 方針 : このパターンの場合は"ビエトの方法"が有効なので、立法完成してからそれを使う

// A = 2 * √(-p/3)、θ = 1/3 * acos((3q/2p)×√(-p/3)) とする

const A = 2 * Math.sqrt(-p / 3);

const theta = (1 / 3) * Math.acos((3 * q / (2 * p)) * Math.sqrt(-3 / p));

// y³ + py + q = 0 の解は y = A * cos(θ), A * cos(θ + 2π/3), A * cos(θ + 4π/3) になる

const y1 = A * Math.cos(theta);

const y2 = A * Math.cos(theta + 2 * Math.PI / 3);

const y3 = A * Math.cos(theta + 4 * Math.PI / 3);

// yを配列に入れる

y_.push(y1, y2, y3);

} else {

// y = u + v と置いて、u³ + v³ = -q, 3uv = p となるようなu, vについて考える

// u³, v³ = -(q / 2) ± √((q/2)² + (p/3)²) となる

const u_3 = -(q / 2) + Math.sqrt((q / 2) ** 2 + (p / 3) ** 3);

const v_3 = -(q / 2) - Math.sqrt((q / 2) ** 2 + (p / 3) ** 3);

// u, vを求める (立方根をとるだけ)

const u = Math.cbrt(u_3);

const v = Math.cbrt(v_3);

// y = u + v

y_.push(u + v);

}

// x = y - B/3 としてxを求める

const x = y_.map(y => y - B / 3);

// return

return x;

}