Root Finding デモ — 方程式の根を求める

calx の求根モジュールは、1D の非線形方程式から多項式の全根まで統一的に扱える。このページでは 3 つのデモで代表的なユースケースを紹介する。

Demo 1: Brent 法で √2 を求める

$f(x) = x^2 - 2$ の正の根を Brent 法で求める。 Brent 法は二分法の確実性と割線法の速さを組み合わせた実用上最良のアルゴリズムである。

=== Brent's Method: f(x) = x^2 - 2 === root = 1.4142135623730951 expected = 1.4142135623730951 converged: true iterations: 5 error: 4.44e-16

わずか 5 回の反復で機械精度 ($\approx 4.4 \times 10^{-16}$) に達している。二分法なら同じ精度に約 52 回かかる。

Demo 2: 3次方程式の因数分解を検証

$p(x) = x^3 - 6x^2 + 11x - 6 = (x-1)(x-2)(x-3)$ の全根を求め、因数分解が正しいことを数値的に検証する。 3次までは Cardano/Ferrari の公式解が使われる。

=== Cubic: x^3 - 6x^2 + 11x - 6 === root[0] = 1+0i root[1] = 2+0i root[2] = 3+0i Verification: p(root) = 0 p(1) = 0 p(2) = 0 p(3) = 0 ✓

findAllRoots は次数 1〜4 では公式解を使うため、反復なしで厳密に近い結果が得られる。全根が Complex<T> で返されるが、この例では虚部がすべて 0 (実根のみ)。

Demo 3: 5次方程式の複素根

$x^5 + 1 = 0$ は 1 つの実根 $x = -1$ と 4 つの複素根を持つ。 5次以上では公式解が存在しない (Abel-Ruffini の定理) ため、 findAllRoots は自動的に Jenkins-Traub 法を選択する。

=== Quintic: x^5 + 1 === root[0] = -1+0i root[1] = 0.309017+0.951057i root[2] = 0.309017-0.951057i root[3] = -0.809017+0.587785i root[4] = -0.809017-0.587785i All roots lie on the unit circle: |root| = 1 |root[0]| = 1.000000 |root[1]| = 1.000000 |root[2]| = 1.000000 |root[3]| = 1.000000 |root[4]| = 1.000000 ✓

$x^5 + 1 = 0$ の根は $x = e^{i(2k+1)\pi/5}$ ($k = 0, 1, 2, 3, 4$) であり、すべて単位円上に並ぶ。Jenkins-Traub 法は数値的に最も安定な多項式求根アルゴリズムとして知られている。

ソースコードと実行方法

example_roots.cpp (完全なソースコード)

// example_roots.cpp — Root Finding demo
#include <math/roots/root_finding.hpp>
#include <iostream>
#include <iomanip>
using namespace calx;

int main() {
    std::cout << std::setprecision(16);

    // --- Demo 1: Brent's method ---
    std::cout << "=== Brent's Method: f(x) = x^2 - 2 ===\n";
    auto result = brent_method(
        [](double x) { return x * x - 2.0; },
        1.0, 2.0
    );
    std::cout << "root     = " << result.value << "\n";
    std::cout << "expected = " << std::sqrt(2.0) << "\n";
    std::cout << "converged: " << (result.converged ? "true" : "false") << "\n";
    std::cout << "iterations: " << result.iterations << "\n";

    // --- Demo 2: Cubic roots ---
    std::cout << "\n=== Cubic: x^3 - 6x^2 + 11x - 6 ===\n";
    Polynomial<double> p = {-6.0, 11.0, -6.0, 1.0};
    auto roots = findAllRoots(p);
    for (size_t i = 0; i < roots.size(); ++i)
        std::cout << "root[" << i << "] = " << roots[i] << "\n";
    std::cout << "\nVerification:\n";
    for (auto& r : roots)
        std::cout << "  p(" << r.real() << ") = " << p(r.real()) << "\n";

    // --- Demo 3: Quintic ---
    std::cout << "\n=== Quintic: x^5 + 1 ===\n";
    Polynomial<double> q = {1.0, 0.0, 0.0, 0.0, 0.0, 1.0};
    auto roots5 = findAllRoots(q);
    for (size_t i = 0; i < roots5.size(); ++i)
        std::cout << "root[" << i << "] = " << roots5[i] << "\n";
    std::cout << "\nAll on unit circle:\n";
    for (size_t i = 0; i < roots5.size(); ++i)
        std::cout << "  |root[" << i << "]| = " << abs(roots5[i]) << "\n";

    return 0;
}

API の詳細は Root Finding API リファレンスを参照のこと。

ビルドと実行

cd calx
mkdir build && cd build
cmake .. -G "Visual Studio 17 2022" -A x64
cmake --build . --config Release --target example-roots
examples\Release\example-roots.exe