diff --git a/Project.toml b/Project.toml index ba7dae4c..4bb7917e 100644 --- a/Project.toml +++ b/Project.toml @@ -2,7 +2,7 @@ name = "Polynomials" uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" license = "MIT" author = "JuliaMath" -version = "2.0.8" +version = "2.0.9" [deps] Intervals = "d8418881-c3e1-53bb-8760-2df7ec849ed5" diff --git a/src/Polynomials.jl b/src/Polynomials.jl index ebb5dbcd..62e20acf 100644 --- a/src/Polynomials.jl +++ b/src/Polynomials.jl @@ -29,6 +29,16 @@ include("polynomials/multroot.jl") include("polynomials/ChebyshevT.jl") +# Rational functions +if VERSION >= v"1.2.0" + include("rational-functions/common.jl") + include("rational-functions/rational-function.jl") + include("rational-functions/fit.jl") + #include("rational-transfer-function.jl") + include("rational-functions/plot-recipes.jl") +end + + # compat; opt-in with `using Polynomials.PolyCompat` include("polynomials/Poly.jl") diff --git a/src/common.jl b/src/common.jl index 8e7f1e34..3c0da31c 100644 --- a/src/common.jl +++ b/src/common.jl @@ -1019,6 +1019,7 @@ function Base.isapprox(p1::AbstractPolynomial{T,X}, rtol::Real = (Base.rtoldefault(T, S, 0)), atol::Real = 0,) where {T,X,S,Y} assert_same_variable(p1, p2) + (hasnan(p1) || hasnan(p2)) && return false # NaN poisons comparisons # copy over from abstractarray.jl Δ = norm(p1-p2) if isfinite(Δ) diff --git a/src/polynomials/multroot.jl b/src/polynomials/multroot.jl index cfd72d6f..8b5a4eb9 100644 --- a/src/polynomials/multroot.jl +++ b/src/polynomials/multroot.jl @@ -79,6 +79,9 @@ is misidentified. function multroot(p::Polynomials.StandardBasisPolynomial{T}; verbose=false, kwargs...) where {T} + # degenerate case, constant + degree(p) == 0 && return (values=T[], multiplicities=Int[], κ=NaN, ϵ=NaN) + # degenerate case, all zeros if (nz = findfirst(!iszero, coeffs(p))) == length(coeffs(p)) return (values=zeros(T,1), multiplicities=[nz-1], κ=NaN, ϵ=NaN) @@ -202,7 +205,7 @@ function pejorative_root(p, zs::Vector{S}, ls::Vector{Int}; The multiplicity count may be in error: the initial guess for the roots failed to converge to a pejorative root. """) - return(zₘs) + return(zₖs) end end diff --git a/src/polynomials/ngcd.jl b/src/polynomials/ngcd.jl index 624737f2..ebef5286 100644 --- a/src/polynomials/ngcd.jl +++ b/src/polynomials/ngcd.jl @@ -19,6 +19,7 @@ function ngcd(p::P, q::Q, a,b = divrem(p,q) return ngcd(q,b, args...; λ=100, kwargs...) end + # easy cases degree(p) < 0 && return (u=q, v=p, w=one(q), θ=NaN, κ=NaN) degree(p) == 0 && return (u=one(q), v=p, w=q, θ=NaN, κ=NaN) @@ -28,15 +29,30 @@ function ngcd(p::P, q::Q, Polynomials.assert_same_variable(p,q) R = promote_type(float(T), float(S)) + 𝑷 = Polynomials.constructorof(promote_type(P,Q)){R,X} + ps = R[pᵢ for pᵢ ∈ coeffs(p)] qs = R[qᵢ for qᵢ ∈ coeffs(q)] - p′ = PnPolynomial(ps) - q′ = PnPolynomial(qs) + # cancel zeros + nz = min(findfirst(!iszero, ps), findfirst(!iszero, qs)) + if nz == length(qs) + x = variable(p) + u = x^(nz-1) + v,w = 𝑷(ps[nz:end]), 𝑷(qs[nz:end]) + return (u=u, v=v, w=w, Θ=NaN, κ=NaN) + end + + ## call ngcd + p′ = PnPolynomial{R,X}(ps[nz:end]) + q′ = PnPolynomial{R,X}(qs[nz:end]) out = NGCD.ngcd(p′, q′, args...; kwargs...) - 𝑷 = Polynomials.constructorof(promote_type(P,Q)){R,X} + 𝑷 = Polynomials.constructorof(promote_type(P,Q)){R,X} u,v,w = convert.(𝑷, (out.u,out.v,out.w)) + if nz > 1 + u *= variable(u)^(nz-1) + end (u=u,v=v,w=w, Θ=out.Θ, κ = out.κ) end diff --git a/src/rational-functions/common.jl b/src/rational-functions/common.jl new file mode 100644 index 00000000..6b8ec250 --- /dev/null +++ b/src/rational-functions/common.jl @@ -0,0 +1,573 @@ +export RationalFunction +export poles, residues +export lowest_terms + +## ---- +""" + AbstractRationalFunction{T,X,P} + +Abstract type for holding ratios of polynomials of type `P{T,X}`. + +Default methods for basic arithmetic operations are provided. + +Numeric methods to cancel common factors, compute the poles, and return the residues are provided. + +!!! Note: + Requires `VERSION >= v"1.2.0"` +""" +abstract type AbstractRationalFunction{T,X,P} end + + +function Base.show(io::IO, pq::AbstractRationalFunction) + p,q = pq + print(io,"(") + print(io, p) + print(io, ") // (") + print(io, q) + print(io, ")") +end + +## helper to make a rational function of given type +function rational_function(::Type{R}, p::P, q::Q) where {R<:AbstractRationalFunction, + T,X, P<:AbstractPolynomial{T,X}, + S, Q<:AbstractPolynomial{S,X}} + constructorof(R)(promote(p,q)...) +end + + +## ---- conversion +function Base.convert(::Type{PQ}, pq′::PQ′) where {T,X,P,PQ <: AbstractRationalFunction{T,X,P}, + T′,X′,P′,PQ′<:AbstractRationalFunction{T′,X′,P′} } + !isconstant(pq′) && assert_same_variable(X,X′) + p′,q′=pqs(pq′) + p,q = convert(P, p′), convert(P, q′) + rational_function(PQ, p, q) +end + + +# R ⊂ R[x] ⊂ R(x) +function Base.convert(::Type{PQ}, p::Number) where {PQ <: AbstractRationalFunction} + P = eltype(PQ) + rational_function(PQ, p * one(P), one(P)) +end + +function Base.convert(::Type{PQ}, p::P) where {PQ <: AbstractRationalFunction, P<:AbstractPolynomial} + Q = eltype(PQ) + q = convert(Q, p) + rational_function(PQ, q, one(q)) +end + +function Base.convert(::Type{P}, pq::PQ) where {P<:AbstractPolynomial, PQ<:AbstractRationalFunction} + p,q = pqs(pq) + isconstant(q) || throw(ArgumentError("Can't convert rational function with non-constant denominator to a polynomial.")) + convert(P, p) / constantterm(q) +end + +function Base.convert(::Type{S}, pq::PQ) where {S<:Number, T,X,P,PQ<:AbstractRationalFunction} + !isconstant(pq) && throw(ArgumentError("Can't convert non-constant rational function to a number")) + S(pq(0)) +end + + + +# promotion rule to promote things upwards +function Base.promote_rule(::Type{PQ}, ::Type{PQ′}) where {T,X,P,PQ <: AbstractRationalFunction{T,X,P}, + T′,X′,P′,PQ′<:AbstractRationalFunction{T′,X′,P′} } + assert_same_variable(X,X′) + PQ_, PQ′_ = constructorof(PQ), constructorof(PQ′) + 𝑷𝑸 = PQ_ == PQ′ ? PQ_ : RationalFunction + 𝑷 = constructorof(typeof(variable(P)+variable(P′))) + 𝑻 = promote_type(T,T′) + 𝑷𝑸{𝑻,X,𝑷{𝑻,X}} +end +Base.promote_rule(::Type{PQ}, ::Type{P}) where {PQ <: AbstractRationalFunction, P<:AbstractPolynomial} = PQ +Base.promote_rule(::Type{PQ}, ::Type{P}) where {PQ <: AbstractRationalFunction, P<:Number} = PQ + + +## Look like rational numbers +# The p//q constructor is reserved for the `RationalFunction` type, but these should all be defined +# for other possible types. +function Base.://(p::PQ,q::PQ′) where {PQ <: AbstractRationalFunction, PQ′ <: AbstractRationalFunction} + p0,p1 = p + q0,q1 = q + rational_function(promote_type(PQ, PQ′), p0*q1, p1*q0) +end + +function Base.://(p::Union{Number,AbstractPolynomial},q::PQ) where {PQ <: AbstractRationalFunction} + q0,q1 = q + rational_function(PQ, p*q1, q0) +end +function Base.://(p::PQ, q::Union{Number,AbstractPolynomial}) where {PQ <: AbstractRationalFunction} + p0, p1 = p + rational_function(PQ, p0, p1*q) +end + +function Base.copy(pq::PQ) where {PQ <: AbstractRationalFunction} + p,q = pq + rational_function(PQ, p, q) +end + + +## ---- + +# requires these field names +Base.numerator(pq::AbstractRationalFunction) = pq.num +Base.denominator(pq::AbstractRationalFunction) = pq.den + +# Treat a RationalFunction as a tuple (num=p, den=q) +Base.length(pq::AbstractRationalFunction) = 2 +function Base.iterate(pq, state=nothing) + state == nothing && return (numerator(pq), 1) + state == 1 && return (denominator(pq), 2) + nothing +end +Base.collect(pq::AbstractRationalFunction{T,X,P}) where {T,X,P} = collect(P, pq) + + +Base.eltype(pq::Type{<:AbstractRationalFunction{T,X,P}}) where {T,X,P} = P +Base.eltype(pq::Type{<:AbstractRationalFunction{T,X}}) where {T,X} = Polynomial{T,X} +Base.eltype(pq::Type{<:AbstractRationalFunction{T}}) where {T} = Polynomial{T,:x} +Base.eltype(pq::Type{<:AbstractRationalFunction}) = Polynomial{Float64,:x} + + +""" + pqs(pq) + +Return `(p,q)`, where `pq=p/q`, as polynomials. +Alternative to simply `p,q=pq` in case `pq` is not stored as two polynomials. +""" +pqs(pq::AbstractRationalFunction) = (numerator(pq), denominator(pq)) + +## ---- + +Base.size(F::AbstractRationalFunction) = () +function Base.isinf(pq::AbstractRationalFunction) + p,q=pqs(pq) + iszero(q) && !iszero(p) +end +function Base.isnan(pq::AbstractRationalFunction) + p,q= pqs(pq) + iszero(p) && iszero(q) +end +function Base.iszero(pq::AbstractRationalFunction) + p,q= pqs(pq) + iszero(p) && !iszero(q) +end +function Base.isone(pq::AbstractRationalFunction) + p,q = pqs(pq) + isconstant(p) && isconstant(q) && p == q +end + + + +## ----- + +_indeterminate(::Type{<:AbstractRationalFunction}) = nothing +_indeterminate(::Type{PQ}) where {T,X,PQ<:AbstractRationalFunction{T,X}} = X +indeterminate(pq::PQ) where {PQ<:AbstractRationalFunction} = indeterminate(PQ) +function indeterminate(::Type{PQ}, var=:x) where {PQ<:AbstractRationalFunction} + X′ = _indeterminate(PQ) + X = X′ == nothing ? Symbol(var) : X′ + X +end + +isconstant(pq::AbstractRationalFunction; kwargs...) = all(isconstant.(lowest_terms(pq;kwargs...))) +isconstant(::Number) = true + +function constantterm(pq::AbstractRationalFunction; kwargs...) + p,q = pqs(pq) + isconstant(pq) && return constantterm(p)/constantterm(q) + throw(ArgumentError("No constant term defined for a non-constant polynomial")) +end + +function Base.zero(pq::R) where {R <: AbstractRationalFunction} + p,q = pqs(pq) + rational_function(R, zero(p), one(q)) +end +function Base.one(pq::R) where {R <: AbstractRationalFunction} + p,q = pqs(pq) + rational_function(R, one(p), one(q)) +end +function variable(pq::R) where {R <: AbstractRationalFunction} + p,q = pqs(pq) + rational_function(R, variable(p), one(q)) +end + +function variable(::Type{PQ}) where {PQ <: AbstractRationalFunction} + x = variable(eltype(PQ)) + rational_function(PQ, x, one(x)) +end + + +# use degree as largest degree of p,q after reduction +function degree(pq::AbstractRationalFunction) + pq′ = lowest_terms(pq) + maximum(degree.(pq′)) +end + +# Evaluation +function eval_rationalfunction(x, pq::AbstractRationalFunction) + md = minimum(degree.(pq)) + num, den = pqs(pq) + result = num(x)/den(x) + while md >= 0 + !isnan(result) && return result + num,den = derivative(num), derivative(den) + result = num(x)/den(x) + md -= 1 + end + + x*NaN +end + +# equality +import Base: == +function ==(p::AbstractRationalFunction{T,X,P}, q::AbstractRationalFunction{S,Y,Q}) where {T,X,P,S,Y,Q} + isconstant(p) && isconstant(q) && p(0) == q(0) && return true + X == Y || return false + p₀, p₁ = pqs(p) + q₀, q₁ = pqs(q) + p₀ * q₁ == q₀ * p₁ || return false +end + +function ==(p::AbstractRationalFunction{T,X,P}, q::Union{AbstractPolynomial, Number}) where {T,X,P} + ==(promote(p,q)...) +end + +function ==(p::Union{AbstractPolynomial, Number}, q::AbstractRationalFunction{T,X,P}) where {T,X,P} + ==(promote(p,q)...) +end + + +function Base.isapprox(pq₁::PQ₁, pq₂::PQ₂, + rtol::Real = sqrt(eps(real(promote_type(T,S)))), + atol::Real = zero(real(promote_type(T,S)))) where {T,X,P,PQ₁<:AbstractRationalFunction{T,X,P}, + S,Y,Q,PQ₂<:AbstractRationalFunction{S,Y,Q}} + + p₁,q₁ = pqs(pq₁) + p₂,q₂ = pqs(pq₂) + + isapprox(p₁*q₂, q₁*p₂; rtol=rtol, atol=atol) +end + + +# Arithmetic +function Base.:-(pq::PQ) where {PQ <: AbstractRationalFunction} + p, q = copy.(pq) + rational_function(PQ, -p, q) +end + +Base.:+(p::Number, q::AbstractRationalFunction) = q + p +Base.:+(p::AbstractRationalFunction, q::Number) = p + q*one(p) +Base.:+(p::AbstractPolynomial, q::AbstractRationalFunction) = q + p +Base.:+(p::AbstractRationalFunction, q::AbstractPolynomial) = p + (q//one(q)) +Base.:+(p::AbstractRationalFunction, q::AbstractRationalFunction) = sum(promote(p,q)) +# type should implement this +function Base.:+(p::R, q::R) where {T,X,P,R <: AbstractRationalFunction{T,X,P}} + p0,p1 = pqs(p) + q0,q1 = pqs(q) + rational_function(R, p0*q1 + p1*q0, p1*q1) +end + +Base.:-(p::Number, q::AbstractRationalFunction) = -q + p +Base.:-(p::AbstractRationalFunction, q::Number) = p - q*one(p) +Base.:-(p::AbstractPolynomial, q::AbstractRationalFunction) = -q + p +Base.:-(p::PQ, q::AbstractPolynomial) where {PQ <: AbstractRationalFunction} = p - rational_function(PQ,q, one(q)) +function Base.:-(p::AbstractRationalFunction, q::AbstractRationalFunction) + p′, q′ = promote(p,q) + p′ - q′ +end +# type should implement this +function Base.:-(p::R, q::R) where {T,X,P,R <: AbstractRationalFunction{T,X,P}} + p0,p1 = pqs(p) + q0,q1 = pqs(q) + rational_function(R, p0*q1 - p1*q0, p1*q1) +end + +function Base.:*(p::Number, q::R) where {T, X, R <: AbstractRationalFunction{T,X}} + q0,q1 = pqs(q) + rational_function(R, (p*q0), q1) +end +function Base.:*(p::R, q::Number) where {R <:AbstractRationalFunction} + p0,p1 = pqs(p) + rational_function(R, p0*q, p1) +end +function Base.:*(p::AbstractPolynomial, q::R) where {R <: AbstractRationalFunction} + rational_function(R, p, one(p)) * q +end +Base.:*(p::R, q::AbstractPolynomial) where {R <: AbstractRationalFunction} = p * rational_function(R,q, one(q)) +# type should implement this +Base.:*(p::AbstractRationalFunction, q::AbstractRationalFunction) = prod(promote(p,q)) +function Base.:*(p::R, q::R) where {T,X,P,R <: AbstractRationalFunction{T,X,P}} + p0,p1 = pqs(p) + q0,q1 = pqs(q) + rational_function(R, p0*q0, p1*q1) +end + + + +function Base.:/(p::Number, q::R) where {R <: AbstractRationalFunction} + q0,q1 = pqs(q) + rational_function(R,p*q1, q0) +end +function Base.:/(p::R, q::Number) where {R <: AbstractRationalFunction} + p0,p1 = pqs(p) + rational_function(R, p0, (p1*q)) +end +Base.:/(p::AbstractPolynomial, q::PQ) where {PQ <: AbstractRationalFunction} = rational_function(PQ, p,one(p)) / q +function Base.:/(p::PQ, q::AbstractPolynomial) where {PQ <: AbstractRationalFunction} + p0,p1 = pqs(p) + rational_function(PQ,p0, p1*q) +end +function Base.:/(p::AbstractRationalFunction, q::AbstractRationalFunction) + p′,q′ = promote(p,q) + p′ / q′ +end +# type should implement this +function Base.:/(p::PQ, q::PQ) where {T,X,P,PQ <: AbstractRationalFunction{T,X,P}} + p0,p1 = pqs(p) + q0,q1 = pqs(q) + rational_function(PQ, p0*q1, p1*q0) +end + +function Base.:^(pq::P, n::Int) where {P <: AbstractRationalFunction} + p,q = pqs(pq) + rational_function(P, p^n, q^n) +end + +function Base.inv(pq::P) where {P <: AbstractRationalFunction} + p,q = pqs(pq) + rational_function(P, q, p) +end + +# conj, transpose... TODO + +## derivative and integrals +function derivative(pq::P, n::Int=1) where {P <: AbstractRationalFunction} + n <= 0 && return pq + while n >= 1 + p,q = pqs(pq) + pq = rational_function(P, (derivative(p)*q - p * derivative(q)), q^2) + n -= 1 + end + pq +end + +function integrate(pq::P) where {P <: AbstractRationalFunction} + p,q = pqs(pq) + isconstant(q) && return rational_function(P, integrate(p), q) + # XXX could work here, e.g.: + # d,r = partial_fraction + # ∫d + Σ∫r for each residue (issue with logs) + throw(ArgumentError("Can only integrate rational functions with constant denominators")) +end + +## ---- + +""" + divrem(pq::AbstractRationalFunction; method=:numerical, kargs...) + +Return `d,r` with `p/q = d + r/q` where `degree(numerator(r)) < degree(denominator(q))`, `d` a Polynomial, `r` a `AbstractRationalFunction`. + +* `method`: passed to `gcd` +* `kwargs...`: passed to `gcd` +""" +function Base.divrem(pq::PQ; method=:numerical, kwargs...) where {PQ <: AbstractRationalFunction} + p,q = pqs(pq) + degree(p) < degree(q) && return (zero(p), pq) + + d,r = divrem(p,q) + (d, rational_function(PQ, r, q)) + +end + + +# like Base.divgcd in rational.jl +# divide p,q by u +function _divgcd(V::Val{:euclidean}, pq; kwargs...) + p, q = pqs(pq) + u = gcd(V,p,q; kwargs...) + p÷u, q÷u +end +function _divgcd(V::Val{:noda_sasaki}, pq; kwargs...) + p, q = pqs(pq) + u = gcd(V,p,q; kwargs...) + p÷u, q÷u +end +function _divgcd(v::Val{:numerical}, pq; kwargs...) + u,v,w,θ,κ = ngcd(pqs(pq)...; kwargs...) # u⋅v=p, u⋅w=q + v, w +end + + +""" + lowest_terms(pq::AbstractRationalFunction, method=:numerical) + +Find GCD of `(p,q)`, `u`, and return `(p÷u)//(q÷u)`. Commonly referred to as lowest terms. + +* `method`: passed to `gcd(p,q)` +* `kwargs`: passed to `gcd(p,q)` + +By default, `AbstractRationalFunction` types do not cancel common factors. This method will numerically cancel common factors, returning the normal form, canonicalized here by `q[end]=1`. The result and original may be considered equivalent as rational expressions, but different when seen as functions of the indeterminate. + +""" +function lowest_terms(pq::PQ; method=:numerical, kwargs...) where {T,X, + P<:StandardBasisPolynomial{T,X}, + PQ<:AbstractRationalFunction{T,X,P}} + v,w = _divgcd(Val(method), pq; kwargs...) + rational_function(PQ, v/w[end], w/w[end]) +end + +## ---- zeros, poles, ... +""" + poles(pq::AbstractRationalFunction; method=:numerical, kwargs...) + +For a rational function `p/q`, first reduces to normal form, then finds the roots and multiplicities of the resulting denominator. + +""" +function poles(pq::AbstractRationalFunction; method=:numerical, kwargs...) + pq′ = lowest_terms(pq; method=method, kwargs...) + den = denominator(pq′) + mr = Multroot.multroot(den) + (zs=mr.values, multiplicities = mr.multiplicities) +end + +""" + roots(pq::AbstractRationalFunction; kwargs...) + +Return the `zeros` of the rational function (after cancelling commong factors, the `zeros` are the roots of the numerator. + +""" +function roots(pq::AbstractRationalFunction; method=:numerical, kwargs...) + pq′ = lowest_terms(pq; method=method, kwargs...) + den = numerator(pq′) + mr = Multroot.multroot(den) + (zs=mr.values, multiplicities = mr.multiplicities) +end + +""" + residues(pq::AbstractRationalFunction; method=:numerical, kwargs...) + +If `p/q =d + r/q`, returns `d` and the residues of a rational fraction `r/q`. + +First expresses `p/q =d + r/q` with `r` of lower degree than `q` through `divrem`. +Then finds the poles of `r/q`. +For a pole, `λj` of multiplicity `k` there are `k` residues, +`rⱼ[k]/(z-λⱼ)^k`, `rⱼ[k-1]/(z-λⱼ)^(k-1)`, `rⱼ[k-2]/(z-λⱼ)^(k-2)`, …, `rⱼ[1]/(z-λⱼ)`. +The residues are found using this formula: +`1/j! * dʲ/dsʲ (F(s)(s - λⱼ)^k` evaluated at `λⱼ` ([5-28](https://stanford.edu/~boyd/ee102/rational.pdf)). + + +## Example + +(From page 5-33 of above pdf) + +```jldoctest +julia> s = variable(Polynomial, :s) +Polynomial(1.0*s) + +julia> pq = (-s^2 + s + 1) // (s * (s+1)^2) +(1.0 + 1.0*s - 1.0*s^2) // (1.0*s + 2.0*s^2 + 1.0*s^3) + +julia> d,r = residues(pq) +(Polynomial(0.0), Dict(-1.0 => [-1.9999999999999978, 1.0000000000000029], 2.0153919257182735e-18 => [1.0])) + +julia> iszero(d) +true + +julia> z = variable(pq) +(1.0*s) // (1.0) + +julia> for (λ, rs) ∈ r # reconstruct p/q from output of `residues` + for (i,rᵢ) ∈ enumerate(rs) + d += rᵢ/(z-λ)^i + end + end + +julia> p′, q′ = lowest_terms(d) +(1.0 + 1.0*s - 1.0*s^2) // (6.64475e-16 + 1.0*s + 2.0*s^2 + 1.0*s^3) + +julia> q′ ≈ (s * (s+1)^2) # works, as q is monic +true + +julia> p′ ≈ (-s^2 + s + 1) +true +``` + + + +!!! Note: + There are several areas where numerical issues can arise. The `divrem`, the identification of multiple roots (`multroot`), the evaluation of the derivatives, ... + +""" +function residues(pq::AbstractRationalFunction; method=:numerical, kwargs...) + + + d,r′ = divrem(pq) + r = lowest_terms(r′; method=method, kwargs...) + b,a = pqs(r) + a′ = derivative(a) + + residues = Any[] + mr = Multroot.multroot(a) + + for (λₖ, mₖ) ∈ zip(mr.values, mr.multiplicities) + + if mₖ == 1 + push!(residues, λₖ => [b(λₖ)/a′(λₖ)]) + else + # returns rₖ₁, rₖ₂,...rₖₘ where rₖ,i/(s-λₖ)ⁱ is part of the decomposition + s = variable(a) + F = lowest_terms(r*(s-λₖ)^mₖ) + rs = [F(λₖ)] + j! = 1 + for j ∈ 1:mₖ-1 + dF = lowest_terms(derivative(F)) + pushfirst!(rs, 1/j! * dF(λₖ)) + j! *= (j+1) + end + push!(residues, λₖ => rs) + end + end + d, Dict(residues...) +end + +""" + partial_fraction(pq::AbstractRationalFunction; method=:numerical, kwargs...) + +For a rational function `p/q = d + r/q`, with `degree(r) < degree(q)` returns `d` and +the terms that comprise `r/q`: For each pole with multiplicity returns +`rⱼ[k]/(z-λⱼ)^k`, `rⱼ[k-1]/(z-λⱼ)^(k-1)`, `rⱼ[k-2]/(z-λⱼ)^(k-2)`, …, `rⱼ[1]/(z-λⱼ)`. + +Should be if `p/q` is in normal form and `d,r=partial_fraction(p//q)` that +`d + sum(r) - p//q ≈ 0` + +""" +function partial_fraction(pq::AbstractRationalFunction; method=:numerical, kwargs...) + d,r = residues(pq; method=method, kwargs...) + s = variable(pq) + d, partial_fraction(Val(:residue), r, s) +end + +function partial_fraction(::Val{:residue}, r, s::PQ) where {PQ} + terms = [] + for (λₖ,rsₖ) ∈ r + for (rⱼ,j) ∈ zip(rsₖ, length(rsₖ):-1:1) + push!(terms, rⱼ/(s-λₖ)^j) + end + end + terms + +end + + + + + + + + + + + + diff --git a/src/rational-functions/fit.jl b/src/rational-functions/fit.jl new file mode 100644 index 00000000..c1f7557e --- /dev/null +++ b/src/rational-functions/fit.jl @@ -0,0 +1,238 @@ +module RationalFunctionFit +using ..Polynomials +import ..Polynomials: RationalFunction, indeterminate, constructorof +import ..Polynomials: evalpoly +using LinearAlgebra + +""" + fit(::Type{RationalFunction}, xs::AbstractVector{S}, ys::AbstractVector{T}, m, n; var=:x) + +Fit a rational function of the form `pq = (a₀ + a₁x¹ + … + aₘxᵐ) / (1 + b₁x¹ + … + bₙxⁿ)` to the data `(x,y)`. + + + +!!! Note: + + This uses a simple implementation of the Gauss-Newton method + to solve the non-linear least squares problem: + `minᵦ Σ(yᵢ - pq(xᵢ,β)²`, where `β=(a₀,a₁,…,aₘ,b₁,…,bₙ)`. + + A more rapidly convergent method is used in the `LsqFit.jl` + package, and if performance is important, re-expressing the + problem for use with that package is suggested. + + Further, if an accurate rational function fit of adaptive degrees + is of interest, the `BaryRational.jl` package provides an + implementation of the AAA algorithm ("which offers speed, + flexibility, and robustness we have not seen in other algorithms" + [Nakatsukasa, Sète, + Trefethen](https://arxiv.org/pdf/1612.00337.pdf)) and one using + Floater-Hormann weights [Floater, + Hormann](https://doi.org/10.1007/s00211-007-0093-y) ("that have no + real poles and arbitrarily high approximation orders on any real + interval, regardless of the distribution of the points") + + The [RationalApproximations](https://github.com/billmclean/RationalApproximations) package also has implementations of the AAA algorithm. + + A python libary, [polyrat](https://github.com/jeffrey-hokanson/polyrat), has implementations of other algorithms. + +## Example +``` +julia> x = variable(Polynomial{Float64}) +Polynomial(1.0*x) + +julia> pq = (1+x)//(1-x) +(1.0 + 1.0*x) // (1.0 - 1.0*x) + +julia> xs = 2.0:.1:3; + +julia> ys = pq.(xs); + +julia> v = fit(RationalFunction, xs, ys, 2, 2) +(1.0 + 1.0*x - 6.82121e-13*x^2) // (1.0 - 1.0*x + 2.84217e-13*x^2) + +julia> maximum(abs, v(x)-pq(x) for x ∈ 2.1:0.1:3.0) +1.06314956838105e-12 + +julia> using BaryRational + +julia> u = aaa(xs,ys) +(::BaryRational.AAAapprox{Vector{Float64}}) (generic function with 1 method) + +julia> maximum(abs, u(x)-pq(x) for x ∈ 2.1:0.1:3.0) +4.440892098500626e-16 + +julia> u(variable(pq)) # to see which polynomial is used +(2.68328 + 0.447214*x - 1.78885*x^2 + 0.447214*x^3) // (2.68328 - 4.91935*x + 2.68328*x^2 - 0.447214*x^3) +``` + +""" +function Polynomials.fit(::Type{PQ}, xs::AbstractVector{S}, ys::AbstractVector{T}, m, n; var=:x) where {T,S, PQ<:RationalFunction} + + + β₁,β₂ = gauss_newton(collect(xs), convert(Vector{float(T)}, ys), m, n) + P = eltype(PQ) + T′ = Polynomials._eltype(P) == nothing ? eltype(β₁) : eltype(P) + X = indeterminate(PQ, var) + P′ = constructorof(P){T′,X} + p = P′(β₁) + q = P′(vcat(1, β₂)) + + p // q +end + + +""" + fit(::Type{RationalFunction}, r::Polynomial, m, n; var=:x) + +Fit a Pade approximant ([`pade_fit`](@ref)) to `r`. + +Examples: + +```jldoctext +julia> using Polynomials, PolynomialRatios + +julia> x = variable() +Polynomial(x) + +julia> ex = 1 + x + x^2/2 + x^3/6 + x^4/24 + x^5/120 # Taylor polynomial for e^x +Polynomial(1.0 + 1.0*x + 0.5*x^2 + 0.16666666666666666*x^3 + 0.041666666666666664*x^4 + 0.008333333333333333*x^5) + +julia> maximum(abs, exp(x) - fit(RationalFunction, ex, 1,1)(x) for x ∈ 0:.05:0.5) +0.017945395966538547 + +julia> maximum(abs, exp(x) - fit(RationalFunction, ex, 1,2)(x) for x ∈ 0:.05:0.5) +0.0016624471707165078 + +julia> maximum(abs, exp(x) - fit(RationalFunction, ex, 2,1)(x) for x ∈ 0:.05:0.5) +0.001278729299871717 + +julia> maximum(abs, exp(x) - fit(RationalFunction, ex, 2,2)(x) for x ∈ 0:.05:0.5) +7.262205147950951e-5 +``` +""" +function Polynomials.fit(::Type{RationalFunction},r::Polynomial, m::Integer, n::Integer;var=:x) + p,q = pade_fit(r, m,n, var=var) + p // q +end + + + +## ---- Pade +## https://mathworld.wolfram.com/PadeApproximant.html +""" + pade_fit(r::Polynomial, m,n) + +For a polynomial `r` of degree `d ≥ m + n`, find a rational function `p/q` with +`degree(p) ≤ m`, `degree(q) ≤ n` and `q*r - p = x^{m+n+1}*s(x)` for some polynomial `s`. + +This implementation sets up a system of equations to identify `p` and `q`. +""" +function pade_fit(p::Polynomial{T}, m::Integer, n::Integer; var=:x) where {T} + d = degree(p) + @assert (0 <= m) && (1 <= n) && (m + n <= d) + + # could be much more perfomant + c = convert(LaurentPolynomial, p) # for better indexing + cs = [c[m+j-i] for j ∈ 1:n, i ∈ 0:n] + + qs′ = cs[:, 2:end] \ cs[:,1] + qs = vcat(1, -qs′) + + cs = [c[0 + j - i] for j ∈ 0:m, i∈0:n] + ps = cs * qs + + Polynomial(ps, var), Polynomial(qs,var) +end + + +## ---- Least Squares +## avoiding dependency on another package, LsqFit + +using LinearAlgebra + +# return pair ((a₀,...,aₙ), (b₁,...,bₘ)) +function initial_guess(xs::Vector{T}, ys::Vector{S}, n, m) where {T, S} + # p(x) = a₀ + ... + aₙx^n + # q(x) = 1 + b₁x + ... + bₘx^m = 1 + r(x) + # yᵢ + yᵢ * r(xᵢ) = p(xᵢ) + # yᵢ = p(xᵢ) - r(xᵢ) * yᵢ + k = n+1+m + A = zeros(T, k, k) + A[:,1] .= one(T) + xs′ = xs[1:k] + xs′′ = copy(xs′) + for i ∈ 1:n + A[:,1+i] = xs′ + xs′ .*= xs′′ + end + xs′ = -copy(xs′′) .* ys[1:k] + for i ∈ 1:m + A[:, 1+n+i] = xs′ + xs′ .*= xs′′ + end + β = pinv(A) * ys[1:k] + +end + +function make_model(n) + (x, β) -> begin + β₁, β₂ = β[1:n+1], β[n+2:end] + evalpoly.(x,(β₁,)) ./ evalpoly.(x, (vcat(1, β₂),)) + end +end + +function J!(Jᵣ, xs::Vector{T}, β, n) where {T} + β₁, β₂ = β[1:n+1],β[n+2:end] + ps = evalpoly.(xs, (β₁,)) + qs = evalpoly.(xs, (vcat(1, β₂),)) + λ = one(T) ./ qs + for i ∈ eachindex(β₁) + Jᵣ[:,1] = λ + λ .*= xs + end + λ = xs .* ps ./ (qs .* qs) + for i ∈ eachindex(β₂) + Jᵣ[:, n+1+i]=λ + λ .*= xs + end + nothing +end + + +function gauss_newton(xs, ys::Vector{T}, n, m, tol=sqrt(eps(T))) where {T} + + β = initial_guess(xs, ys, n, m) + model = make_model(n) + + Jᵣ = zeros(T, length(xs), 1 + n + m) + + Δ = norm(ys, Inf) * tol + + ϵₘ = norm(ys - model(xs, β), Inf) + βₘ = copy(β) + + no_steps = 0 + + while no_steps < 25 + no_steps += 1 + + r = ys - model(xs, β) + ϵ = norm(r, Inf) + ϵ < Δ && return (β[1:n+1], β[n+2:end]) + if ϵ < ϵₘ + ϵₘ = ϵ + βₘ .= β + end + J!(Jᵣ, xs, β, n) + Δᵦ = pinv(Jᵣ' * Jᵣ) * (Jᵣ' * r) + β .-= Δᵦ + + end + + @warn "no convergence; returning best fit of many steps" + return (βₘ[1:n+1], βₘ[n+2:end]) +end + + +end diff --git a/src/rational-functions/plot-recipes.jl b/src/rational-functions/plot-recipes.jl new file mode 100644 index 00000000..a5d9391c --- /dev/null +++ b/src/rational-functions/plot-recipes.jl @@ -0,0 +1,66 @@ +## Plot recipe +## define a hueristic to work around asymptotes +## just sort of succesful +@recipe function f(pq::AbstractRationalFunction{T}, a=nothing, b=nothing) where {T} + + xlims = get(plotattributes, :xlims, (nothing,nothing)) + ylims = get(plotattributes, :ylims, (nothing, nothing)) + rational_function_trim(pq, a, b, xlims, ylims) + +end + +isapproxreal(x::Real) = true +isapproxreal(x::Complex{T}) where {T} = imag(x) <= sqrt(eps(real(T))) +function toobig(pq) + x -> begin + y = pq(x) + isinf(y) && return true + isnan(y) && return true + abs(y) > 1e8 && return true + return false + end +end + +function rational_function_trim(pq, a, b, xlims, ylims) + + p,q = lowest_terms(//(pq...), method=:numerical) + dpq = derivative(p//q) + p′,q′ = lowest_terms(dpq) + + λs = Multroot.multroot(q).values + λs = isempty(λs) ? λs : real.(filter(isapproxreal, λs)) + + cps = Multroot.multroot(p′).values + cps = isempty(cps) ? cps : real.(filter(isapproxreal, cps)) + cps = isempty(cps) ? cps : filter(!toobig(pq), cps) + + a = a == nothing ? xlims[1] : a + b = b == nothing ? xlims[2] : b + + if a==nothing && b==nothing + u= poly_interval(p) + v= poly_interval(q) + a,b = min(first(u), first(v)), max(last(u), last(v)) + + if !isempty(λs) + a,b = min(a, real(minimum(λs))), max(b, real(maximum(λs))) + end + if !isempty(cps) + a,b = min(a, real(minimum(cps))), max(b, real(maximum(cps))) + end + a *= (a > 0 ? 1/1.5 : 1.25) + b *= (b < 0 ? 1/1.5 : 1.25) + end + + n = 601 + xs = range(a,stop=b, length=n) + ys = pq.(xs) + M = max(5, 3*maximum(abs, pq.(cps)), 1.25*maximum(abs, pq.((a,b)))) + + lo = ylims[1] == nothing ? -M : ylims[1] + hi = ylims[2] == nothing ? M : ylims[2] + ys′ = [lo <= yᵢ <= hi ? yᵢ : NaN for yᵢ ∈ ys] + xs, ys′ + +end + diff --git a/src/rational-functions/rational-function.jl b/src/rational-functions/rational-function.jl new file mode 100644 index 00000000..ffa570b6 --- /dev/null +++ b/src/rational-functions/rational-function.jl @@ -0,0 +1,67 @@ +""" + RationalFunction(p::AbstractPolynomial, q::AbstractPolynomial) + p // q + +Create a rational expression (`p/q`) from the two polynomials. + +There is no attempt to cancel common factors. The [`lowest_terms`](@ref) function attempts to do that. + +For purposes of iteration, a rational function is treated like a tuple. + +## Examples +``` +julia> using Polynomials + +julia> p,q = fromroots(Polynomial, [1,2,3]), fromroots(Polynomial, [2,3,4]) +(Polynomial(-6 + 11*x - 6*x^2 + x^3), Polynomial(-24 + 26*x - 9*x^2 + x^3)) + +julia> pq = p // q +(-6 + 11*x - 6*x^2 + x^3) // (-24 + 26*x - 9*x^2 + x^3) + +julia> lowest_terms(pq) +(-0.333333 + 0.333333*x) // (-1.33333 + 0.333333*x) + +julia> pq(2.5) +-1.0 + +julia> pq(2) # uses first non-`0/0` ratio of `p⁽ᵏ⁾/q⁽ᵏ⁾` +-0.5 + +julia> pq^2 +(36 - 132*x + 193*x^2 - 144*x^3 + 58*x^4 - 12*x^5 + x^6) // (576 - 1248*x + 1108*x^2 - 516*x^3 + 133*x^4 - 18*x^5 + x^6) + +julia> derivative(pq) +(-108 + 180*x - 111*x^2 + 30*x^3 - 3*x^4) // (576 - 1248*x + 1108*x^2 - 516*x^3 + 133*x^4 - 18*x^5 + x^6) +``` + +!!! Note: + The [RationalFunctions.jl](https://github.com/aytekinar/RationalFunctions.jl) was a helpful source of ideas. + + +""" +struct RationalFunction{T, X, P<:AbstractPolynomial{T,X}} <: AbstractRationalFunction{T,X,P} + num::P + den::P + function RationalFunction(p::P, q::P) where {T,X, P<:AbstractPolynomial{T,X}} + new{T,X,P}(p, q) + end + function RationalFunction(p::P, q::T) where {T,X, P<:AbstractPolynomial{T,X}} + new{T,X,P}(p, q*one(P)) + end + function RationalFunction(p::T, q::Q) where {T,X, Q<:AbstractPolynomial{T,X}} + new{T,X,Q}(p*one(Q), q) + end +end + +RationalFunction(p,q) = RationalFunction(promote(p,q)...) +RationalFunction(p::ImmutablePolynomial,q::ImmutablePolynomial) = throw(ArgumentError("Sorry, immutable #polynomials are not a valid polynomial type for RationalFunction")) + +# evaluation +(pq::RationalFunction)(x) = eval_rationalfunction(x, pq) + +# Look like rational numbers +function Base.://(p::AbstractPolynomial,q::AbstractPolynomial) + RationalFunction(p,q) +end + + diff --git a/src/rational-functions/rational-transfer-function.jl b/src/rational-functions/rational-transfer-function.jl new file mode 100644 index 00000000..ad6372a2 --- /dev/null +++ b/src/rational-functions/rational-transfer-function.jl @@ -0,0 +1,169 @@ +## An example subtype of AbstractRationalFunction +## that might prove useful as a guide +module TransferFunction +# https://github.com/andreasvarga/DescriptorSystems.jl/blob/main/src/types/RationalFunction.jl + +using Polynomials +import Polynomials: AbstractPolynomial, AbstractRationalFunction, RationalFunction +import Polynomials: pqs + +export RationalTransferFunction +export sampling_time, gain + +struct RationalTransferFunction{T,X,P<:AbstractPolynomial{T,X},Ts} <: AbstractRationalFunction{T,X,P} + num::P + den::P + function RationalTransferFunction{T,X,P,Ts}(num::P, den::P) where{T,X,P<:AbstractPolynomial{T,X}, Ts} + check_den(den) + new{T,X,P,Ts}(num,den) + end + function RationalTransferFunction{T,X,P,Ts}(num::P, den::P,ts::Union{Real,Nothing}) where{T,X,P<:AbstractPolynomial{T,X}, Ts} + check_den(den) + check_Ts(Ts,ts) + new{T,X,P,Ts}(num,den) + end + # can promote constants to polynomials too + function RationalTransferFunction{T,X,P,Ts}(num::S, den::P, ts::Union{Real,Nothing}) where{S, T,X,P<:AbstractPolynomial{T,X}, Ts} + check_den(den) + check_Ts(Ts,ts) + new{T,X,P,Ts}(convert(P,num),den) + end + function RationalTransferFunction{T,X,P,Ts}(num::P,den::S, ts::Union{Real,Nothing}) where{S, T,X,P<:AbstractPolynomial{T,X}, Ts} + check_den(den) + check_Ts(Ts,ts) + new{T,X,P,Ts}(num, convert(P,den)) + end + function RationalTransferFunction{T,X,P}(num::P, den::P, Ts::Union{Real,Nothing}) where{T,X,P<:AbstractPolynomial{T,X}} + + check_den(den) + Ts′ = standardize_Ts(Ts) + new{T,X,P,Val(Ts′)}(num,den) + end +end + +(pq::RationalTransferFunction)(x) = Polynomials.eval_rationalfunction(x, pq) + +function Base.convert(::Type{PQ}, pq::RationalTransferFunction) where {PQ <:RationalFunction} + p,q = pq + p//q +end + +# alternate constructor +function RationalTransferFunction(p′::P, q′::Q, Ts::Union{Real,Nothing}) where {T,X,P<:AbstractPolynomial{T,X}, + S, Q<:AbstractPolynomial{S,X}} + + p,q = promote(p′, q′) + R = eltype(p) + RationalTransferFunction{R,X,typeof(p)}(p,q,Ts) +end + + +function Polynomials.rational_function(::Type{PQ}, p::P, q::Q) where {PQ <:RationalTransferFunction, + T,X, P<:AbstractPolynomial{T,X}, + S, Q<:AbstractPolynomial{S,X}} + RationalTransferFunction(promote(p,q)..., sampling_time(PQ)) +end + + + +## helpers for constructors +# standardize Ts or throw error +function standardize_Ts(Ts) + isnothing(Ts) || Ts >= 0 || Ts == -1 || + throw(ArgumentError("Ts must be either a positive number, 0 (continuous system), or -1 (unspecified)")) + Ts′ = isnothing(Ts) ? Ts : Float64(Ts) +end +function check_Ts(Ts, ts) + ValT(Ts) == promote_Ts(ValT(Ts), ts) || throw(ArgumentError("sampling times have mismatch")) +end +function check_den(den) + iszero(den) && throw(ArgumentError("Cannot create a rational function with zero denominator")) +end + +ValT(::Val{T}) where {T} = T +sampling_time(pq::RationalTransferFunction{T,X,P,Ts}) where {T,X,P,Ts} = ValT(Ts) +sampling_time(::Type{𝑷}) where {T,X,P,Ts, 𝑷<:RationalTransferFunction{T,X,P,Ts}} = ValT(Ts) + + + +## ---- + +function Base.convert(PQ′::Type{PQ}, p::P) where {PQ <: RationalTransferFunction, P<:AbstractPolynomial} + PQ(p, one(p), sampling_time(PQ)) +end +function Base.convert(PQ′::Type{PQ}, p::Number) where {PQ <: RationalTransferFunction} + PQ(p, one(eltype(PQ)), sampling_time(PQ)) +end + +function promote_Ts(p,q) + Ts1,Ts2 = sampling_time.((p,q)) + promote_Ts(Ts1, Ts2) +end + +function promote_Ts(Ts1::Union{Float64,Nothing}, Ts2::Union{Float64,Nothing}) + isnothing(Ts1) && (return Ts2) + isnothing(Ts2) && (return Ts1) + Ts1 == Ts2 && (return Ts1) + Ts1 == -1 && (Ts2 > 0 ? (return Ts2) : error("Sampling time mismatch")) + Ts2 == -1 && (Ts1 > 0 ? (return Ts1) : error("Sampling time mismatch")) + error("Sampling time mismatch") +end + +function Base.promote_rule(::Type{PQ}, ::Type{PQ′}) where {T,X,P,Ts,PQ <: RationalTransferFunction{T,X,P,Ts}, + T′,X′,P′,Ts′,PQ′ <: RationalTransferFunction{T′,X′,P′,Ts′}} + + S = promote_type(T,T′) + Polynomials.assert_same_variable(X,X′) + Y = X + Q = promote_type(P, P′) + ts = promote_Ts(PQ, PQ′) + RationalTransferFunction{S,Y,Q,Val(ts)} +end + + + + + +# zeros (roots) and poles are defined in common. The gain is specific to this type +""" + gain(pq::RationalTransferFunction) + +The ratio of the leading coefficients +""" +function gain(pq::PQ) where {PQ <: RationalTransferFunction} + p,q = pqs(pq) + p[end]/q[end] +end + + +# need to check here +# +""" + rt = adjoint(r) +Compute the adjoint `rt(λ)` of the rational transfer function `r(λ)` such that for +`r(λ) = num(λ)/den(λ)` we have: + (1) `rt(λ) = conj(num(-λ))/conj(num(-λ))`, if `r.Ts = 0`; + (2) `rt(λ) = conj(num(1/λ))/conj(num(1/λ))`, if `r.Ts = -1` or `r.Ts > 0`. +""" +function Base.adjoint(pq::PQ) where {PQ <: RationalTransferFunction} + p,q = pqs(pq) + Ts = sampling_time(pq) + if Ts != nothing && iszero(Ts) + # p(-λ)/q(-λ) + p′ = poly_scale(p, -1) + q′ = poly_scale(q, -1) + return RationalTransferFunction(p′, q′, Ts) + else + # p(1/λ) / q(1/λ) = poly_inversion(p) / poly_inversion(q) + # maps oo -> 0 + p′ = poly_inversion(p) + q′ = poly_inversion(q) + return RationalTransferFunction(p′, q′, Ts) + end +end + +## XXX confmap ... + + + +end diff --git a/test/StandardBasis.jl b/test/StandardBasis.jl index 84076766..d1a16b13 100644 --- a/test/StandardBasis.jl +++ b/test/StandardBasis.jl @@ -313,6 +313,8 @@ end @test (P([1,eps(), 1]) ≈ P([1,0,1])) == ([1,eps(), 1] ≈ [1,0,1]) # true @test (P([1,2]) ≈ P([1,2,eps()])) == ([1,2,0] ≈ [1,2,eps()]) + # NaN poisons comparison + @test !(P([NaN, 1.0, 2.0]) ≈ P([NaN, 1.0, 2.0])) # check how ==, ===, isapprox ignore variable mismatch when constants are involved, issue #217, issue #219 @test zero(P, :x) == zero(P, :y) @@ -582,6 +584,10 @@ end out = Polynomials.Multroot.multroot(x^3) @test out.values == zeros(T,1) @test out.multiplicities == [3] + # bug with constant + out = Polynomials.Multroot.multroot(P(1)) + @test isempty(out.values) + @test isempty(out.multiplicities) end end end @@ -1044,6 +1050,16 @@ end q = 2.0 + 5.0*x + 8.0*x^2 + 7.0*x^3 + 4.0*x^4 + 1.0*x^5 out = Polynomials.ngcd(p,q) @test out.u * out.v ≈ p + + # check for canceling of x^k terms + x = variable(P{Float64}) + p,q = x^2 + 1, x^2 - 1 + for j ∈ 0:2 + for k ∈ 0:j + out = Polynomials.ngcd(x^j*p, x^k*q) + @test out.u == x^k + end + end end end diff --git a/test/rational-functions.jl b/test/rational-functions.jl new file mode 100644 index 00000000..86624143 --- /dev/null +++ b/test/rational-functions.jl @@ -0,0 +1,169 @@ +using Polynomials +using Test +using LinearAlgebra + +@testset "Basic constructor, arithmetic" begin + + p,q = fromroots(Polynomial, [1,2,3]), fromroots(Polynomial, [2,3,4]) + r,s = SparsePolynomial([1,2,3], :x), SparsePolynomial([1,2,3], :y) + t,u = Polynomial([1,2,pi]), SparsePolynomial([1.1, 2.2, 3.4], :y) + + # constructor + @test p // q isa RationalFunction + @test p // r isa RationalFunction + @test_throws ArgumentError r // s + + # We expect p::P // q::P (same type polynomial). + # As Immutable Polynomials have N as type parameter, we disallow + @test_throws ArgumentError variable(ImmutablePolynomial) // variable(ImmutablePolynomial) + + pq = p // t # promotes to type of t + @test pq isa RationalFunction{Float64, :x} + + # iterations, broadcast + pp,qq = p // q + @test (pp,qq) == (p,q) + @test eltype(p//q) == typeof(pp) + u = gcd(p//q...) + @test u/u[end] ≈ fromroots(Polynomial, [2,3]) + @test degree.(p // q) == [degree(p), degree(q)] + + # evaluation + pq = p//q + @test pq(2.5) ≈ p(2.5) / q(2.5) + @test pq(2) ≈ fromroots([1,3])(2) / fromroots([3,4])(2) + + # arithmetic + rs = r // (r-1) + x = 1.2 + @test (pq + rs)(x) ≈ (pq(x) + rs(x)) + @test (pq * rs)(x) ≈ (pq(x) * rs(x)) + @test (-pq)(x) ≈ -p(x)/q(x) + + # derivative + pq = p // one(p) + x = 1.23 + @test derivative(pq)(x) ≈ derivative(p)(x) + pq = p // q + dp,dq = derivative.((p,q)) + @test derivative(pq)(x) ≈ (dp(x)*q(x) - p(x)*dq(x))/q(x)^2 + + # integral + pq = p // (2*one(p)) + @test iszero(derivative(integrate(pq)) - pq) + pq = p // q + @test_throws ArgumentError integrate(pq) # need logs terms in general + + # lowest terms + pq = p // q + pp, qq = lowest_terms(pq) + @test all(abs.(pp.(roots(qq))) .> 1/2) + + +end + +@testset "zeros, poles, residues" begin + x = variable(Polynomial{Float64}) + p,q = x^7,(x^5 + 4*x^4 + 7*x^3 + 8*x^2 + 5*x + 2) + pq = p // q + + + ## zeros + zs = roots(pq) + @test length(zs.zs) == 1 + @test zs.multiplicities == [7] + + ## poles + ps = poles(pq) + @test length(ps.zs) == 3 + + ## residues + d,r = residues(pq) + @test d ≈ 9 - 4x + x^2 + + z = variable(pq) + for (λ, rs) ∈ r # reconstruct p/q from output of `residues` + for (i,rᵢ) ∈ enumerate(rs) + d += rᵢ/(z-λ)^i + end + end + @test norm(numerator(lowest_terms(d - pq)), Inf) <= sqrt(eps()) + +end + +@testset "As matrix elements" begin + + p, q = Polynomial([1,2], :x), Polynomial([1,2],:y) + pp = p // (p-1) + PP = typeof(pp) + + r, s = SparsePolynomial([1,2], :x), SparsePolynomial([1,2],:y) + rr = r // (r-1) + + ## T, Polynomial{T} promotes + @test eltype([1, p, pp]) == PP + + ## test mixed types promote polynomial type + @test eltype([pp rr p r]) == PP + + ## test non-constant polys of different symbols throw error + @test_throws ArgumentError [pp, q] + @test_throws ArgumentError [pp, q//(q-1)] + + ## constant polys will work -- if specified + @test eltype(PP[pp one(q)]) == PP + @test eltype(PP[pp one(q//(q-1))]) == PP + + ## if not specified, the rational function will not promote + @test eltype([pp one(q)]) == PP + @test_throws ArgumentError [pp one(q//(q-1))] + +end + + +@testset "Rational function fit" begin + + p,q = Polynomial([1,1,1]), Polynomial([1,2]) + xs = range(-0.25,stop=1, length=15) + ys = (p//q).(xs) + pq = fit(RationalFunction, xs, ys, 2, 2) + @test numerator(pq) ≈ p + @test denominator(pq) ≈ q + + x = variable(Polynomial{Float64}) + pq = (1+x)//(1-x) + xs = 2.0:.1:3; + ys = pq.(xs); + v = fit(RationalFunction, xs, ys, 2, 2) + @test maximum(abs, v(x)-pq(x) for x ∈ 2.1:0.1:3.0) <= sqrt(eps()) + +end + +@testset "Pade fit" begin + + #Tests for Pade approximants from Pade.jl + + # exponential + a = Polynomial( 1 .// factorial.(big(0):17) ) + pq = fit(RationalFunction,a,8,8) + @test pq(1.0) ≈ exp(1.0) + @test pq(-1.0) ≈ exp(-1.0) + + # sine + b = Polynomial(Int.(sinpi.((0:16)/2)) .// factorial.(big(0):16)) + pq = fit(RationalFunction, b, 7, 8) + @test pq(1.0) ≈ sin(1.0) + @test pq(-1.0) ≈ sin(-1.0) + + # cosine + c = Polynomial(Int.(sinpi.((1:17)/2)) .// factorial.(big(0):16)) + pq = fit(RationalFunction, c, 8, 8) + @test pq(1.0) ≈ cos(1.0) + @test pq(-1.0) ≈ cos(-1.0) + + # summation of a factorially divergent series + d = Polynomial( (-big(1)).^(0:60) .* factorial.(big(0):60) ) + pq = fit(RationalFunction, d, 30, 30) + @test pq(1.0) ≈ exp(1)*(-Base.MathConstants.γ-sum([(-1).^k/k./factorial(k) for k=1:20])) + +end diff --git a/test/runtests.jl b/test/runtests.jl index 8df9fe87..d6adbf00 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,4 +10,7 @@ using OffsetArrays @testset "Standard basis" begin include("StandardBasis.jl") end @testset "ChebyshevT" begin include("ChebyshevT.jl") end +if VERSION >= v"1.2.0" + @testset "Rational functions" begin include("rational-functions.jl") end +end @testset "Poly, Pade (compatability)" begin include("Poly.jl") end