RationalFunction

Overview

A rational function (RationalFunction<T>) represents the quotient $f(x) = P(x)/Q(x)$ of two polynomials. The numerator $P$ and denominator $Q$ are each stored as a Polynomial<T>, and common factors are cancelled automatically via GCD at construction time (by default).

• Field operations — rational functions form a field; $+ - \times \div$ are closed.
• Derivative, composition, power — quotient-rule differentiation, function composition, integer powers.
• Partial fractions — Heaviside decomposition into $\sum_k A_k/(x - r_k)$ for simple poles.
• Pole-zero cancellation — approximate-GCD minimal realization of a transfer function.

The coefficient type $T$ may be double/float as well as multiprecision Float or exact Rational (root-based functions such as partial fractions and poles require a floating-point type).

Build

#include <math/core/RationalFunction.hpp>                  // rational function
#include <math/core/RationalFunction_partialFractions.hpp> // partial fractions, poles, zeros, residues

Header-only; no library to link. RationalFunction.hpp depends only on Polynomial.hpp.

Construction

Parameters:

Reduction can lower the degree ($\frac{x^2 - 1}{x - 1} = x + 1$):

// (x^2 - 1) / (x - 1)  ascending coeffs: {-1,0,1} / {-1,1}
RationalFunction<double> f({-1, 0, 1}, {-1, 1});
std::cout << f.toString();  // Run output: x+1   (common factor (x-1) cancelled)
std::cout << f(3.0);        // Run output: 4

Accessors & Reduction

Arithmetic

Four operations are defined between two rational functions, a rational function and a scalar $c$, and a rational function and a polynomial (all auto-reduce the result). Compound assignment (+= -= *= /=) and unary - are also provided.

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

Evaluation & Composition

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();         // Run output: (x+1)/(2x+1)
std::cout << h(1.0);               // Run output: 0.666666666666667   (= 2/3)

Derivative, Reciprocal, Power

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

Partial Fractions

When the denominator has simple roots (no repeated roots), the Heaviside method decomposes $f(x) = \text{quotient}(x) + \sum_k \dfrac{A_k}{x - r_k}$, with residues $A_k = P(r_k)/Q'(r_k)$.

Result types:

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

template<typename T>
struct PartialFractionExpansion {
    Polynomial<T>                  quotient;  // polynomial part when deg(P) >= deg(Q)
    std::vector<PartialFraction<T>> terms;     // each partial-fraction term
    bool                          success;   // whether the decomposition succeeded
};

Parameters:

// 1 / (x^2 - 1) = 1/((x-1)(x+1))  ascending: {1} / {-1,0,1}
RationalFunction<double> f({1}, {-1, 0, 1});
auto pfe = partialFractions(f);
std::cout << pfe.success;       // Run output: true
std::cout << pfe.terms.size();  // Run output: 2
std::cout << toString(pfe);     // Run output: (0.5)/(x - (1)) + (-0.5)/(x - (-1))

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

Pole-Zero Cancellation

When working with transfer functions, near-equal pole-zero pairs in numerator and denominator produce a numerically unnecessary degree increase. minimalRealization cancels them via an approximate GCD and returns the minimal realization (e.g. IIR-filter pole-zero cancellation).

Example

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

int main() {
    // (x^2 - 1)/(x - 1) auto-reduces to (x+1)
    RationalFunction<double> f({-1, 0, 1}, {-1, 1});
    std::cout << "f = " << f.toString() << "\n";      // f = x+1

    // Arithmetic: 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

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

    // Partial fractions: 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))
}

Related Modules

Rational functions are built on polynomials and use the root-finding module for poles and zeros.

• Polynomial — numerator/denominator representation and polynomial algebra
• Root Finding — polynomial root finding used for poles and zeros
• Complex — value type for poles, zeros, and residues