RationalFunction — 有理式

概要

有理式 (RationalFunction<T>) は 2 つの多項式の商 $f(x) = P(x)/Q(x)$ を表す型である。 分子 $P$・分母 $Q$ をそれぞれ Polynomial<T> として保持し、 構築時に既定で 最大公約因子 (GCD) を自動約分する。

• 体としての演算 — 有理式は体をなし、四則 ($+ - \times \div$) が閉じている。
• 微分・合成・冪 — 商の微分公式・関数合成・整数冪をサポートする。
• 部分分数分解 — 単純極について Heaviside 法で $\sum_k A_k/(x - r_k)$ に分解する。
• 極零相殺 — 近似 GCD により伝達関数の余分な極零対を相殺し最小実現を得る。

係数型 $T$ は double・float のほか多倍長 Float・有理数 Rational なども使える (部分分数分解・極零など根を扱う関数は浮動小数型が前提)。

ビルド

#include <math/core/RationalFunction.hpp>                  // 有理式 本体
#include <math/core/RationalFunction_partialFractions.hpp> // 部分分数分解・極・零・留数

ヘッダオンリー。リンクするライブラリは不要。RationalFunction.hpp は Polynomial.hpp のみに依存する。

構築

パラメータ:

約分により次数が下がる例 ($\frac{x^2 - 1}{x - 1} = x + 1$):

// (x^2 - 1) / (x - 1)  昇べき係数: {-1,0,1} / {-1,1}
RationalFunction<double> f({-1, 0, 1}, {-1, 1});
std::cout << f.toString();  // 実行結果: x+1   (共通因子 (x-1) を約分)
std::cout << f(3.0);        // 実行結果: 4

アクセサ・約分

四則演算

有理式同士、有理式とスカラー $c$、有理式と多項式の間で四則が定義される (いずれも結果を自動約分)。 複合代入 (+= -= *= /=) と単項 - も用意する。

RationalFunction<double> f({1}, {1, 1});     // 1/(1+x)
RationalFunction<double> g({0, 1}, {1, 1});  // x/(1+x)
std::cout << (f + g).toString();   // 実行結果: 1            (1/(1+x) + x/(1+x) = 1)
std::cout << (f * g).toString();   // 実行結果: x/(x^2+2x+1)
std::cout << (f / g).toString();   // 実行結果: 1/x

評価・合成

RationalFunction<double> f({1}, {1, 1});     // 1/(1+x)
RationalFunction<double> g({0, 1}, {1, 1});  // x/(1+x)
auto h = compose(f, g);            // f(g(x)) = 1 / (1 + x/(1+x)) = (1+x)/(1+2x)
std::cout << h.toString();         // 実行結果: (x+1)/(2x+1)
std::cout << h(1.0);               // 実行結果: 0.666666666666667   (= 2/3)

微分・逆数・冪

RationalFunction<double> f({1}, {1, 1});  // 1/(1+x)
std::cout << f.derivative().toString();   // 実行結果: -1/(x^2+2x+1)   (= -1/(1+x)^2)
std::cout << f.reciprocal().toString();   // 実行結果: x+1
std::cout << f.pow(2).toString();         // 実行結果: 1/(x^2+2x+1)

部分分数分解

分母が単純根 (重根なし) のとき、Heaviside 法で $f(x) = \text{quotient}(x) + \sum_k \dfrac{A_k}{x - r_k}$ に分解する。 留数は $A_k = P(r_k)/Q'(r_k)$。

結果型:

template<typename T>
struct PartialFraction {
    Complex<T> coefficient;  // 留数 A_k
    Complex<T> pole;         // 極 r_k
};

template<typename T>
struct PartialFractionExpansion {
    Polynomial<T>                  quotient;  // deg(P) >= deg(Q) のときの多項式部
    std::vector<PartialFraction<T>> terms;     // 各部分分数項
    bool                          success;   // 分解に成功したか
};

パラメータ:

// 1 / (x^2 - 1) = 1/((x-1)(x+1))  昇べき: {1} / {-1,0,1}
RationalFunction<double> f({1}, {-1, 0, 1});
auto pfe = partialFractions(f);
std::cout << pfe.success;       // 実行結果: true
std::cout << pfe.terms.size();  // 実行結果: 2
std::cout << toString(pfe);     // 実行結果: (0.5)/(x - (1)) + (-0.5)/(x - (-1))

auto ps = poles(f);            // 極 {1, -1}
// ps[0].re = 1,  ps[1].re = -1

極零相殺 (最小実現)

伝達関数を扱う際、分子・分母にほぼ等しい極零対が現れると数値的に不要な次数増加となる。 minimalRealization は近似 GCD でこれを相殺し、最小実現を返す (IIR フィルタの極零相殺など)。

使用例

#include <math/core/RationalFunction.hpp>
#include <math/core/RationalFunction_partialFractions.hpp>
#include <iostream>
using namespace sangi;

int main() {
    // (x^2 - 1)/(x - 1) は (x+1) に自動約分される
    RationalFunction<double> f({-1, 0, 1}, {-1, 1});
    std::cout << "f = " << f.toString() << "\n";      // f = x+1

    // 四則: 1/(1+x) + x/(1+x) = 1
    RationalFunction<double> a({1}, {1, 1}), b({0, 1}, {1, 1});
    std::cout << "a+b = " << (a + b).toString() << "\n"; // a+b = 1

    // 微分: d/dx 1/(1+x) = -1/(1+x)^2
    std::cout << "da  = " << a.derivative().toString() << "\n"; // da = -1/(x^2+2x+1)

    // 部分分数: 1/(x^2-1) = 0.5/(x-1) - 0.5/(x+1)
    RationalFunction<double> r({1}, {-1, 0, 1});
    auto pfe = partialFractions(r);
    if (pfe.success) std::cout << toString(pfe) << "\n";
    // (0.5)/(x - (1)) + (-0.5)/(x - (-1))
}

関連モジュール

有理式は多項式の上に構築され、根の計算に求根モジュールを用いる。

• Polynomial — 多項式 — 分子・分母の表現と多項式演算
• Root Finding — 求根 — 極・零の計算に用いる多項式求根
• Complex — 複素数 — 極・零・留数の値型