JuliaMath / Polynomials.jl

Polynomial manipulations in Julia
http://juliamath.github.io/Polynomials.jl/
Other
303 stars 74 forks source link

Rational function related issues #338

Closed andreasvarga closed 3 years ago

andreasvarga commented 3 years ago

I upgraded the DescriptorSystems to use the new rational function type, by defining according to your suggestion the rational transfer function type as a subtype of the rational function. With this occasion I eliminated all stuff related to concatenation of rational matrices, relying entirely on the existing concatenation techniques. I am convinced my adaptation to the new type is not optimal and some code may even be redundant (duplicating your codes), but in this moment I have no better solution at hand, so I will keep this implementation (at least for a while).

Here are listed some issues I observed:

  1. I implemented for my purposes the conversion from polynomial to rational transfer function. I think this would be useful also for rational functions, since now:
    julia> p = Polynomial([1,2],:s)
    Polynomial(1 + 2*s) 
    RationalFunction(p)
    ERROR: MethodError: no method matching RationalFunction(::Polynomial{Int64, :s})
  2. isapproxfails for rational functions with integer coefficients:
    julia> p = Polynomial([1,2],:s)
    Polynomial(1 + 2*s) 
    isapprox(p/p,p/p)
    ERROR: MethodError: no method matching eps(::Type{Int64})
  3. The length of a rational function is defined as 2 (to allow iteration). However, this imped to perform operations with rational matrices involving broadcasting:
    julia> p/(p+2).*[p/(p+1) p/(p+3)]
    2×2 Matrix{RationalFunction{Int64, :s, Polynomial{Int64, :s}}}:
    (1 + 4*s + 4*s^2) // (2 + 2*s)  (1 + 4*s + 4*s^2) // (4 + 2*s)
    (3 + 8*s + 4*s^2) // (2 + 2*s)  (3 + 8*s + 4*s^2) // (4 + 2*s)

    which is different from the correct result

    julia> [p/(p+1)*p/(p+2) p/(p+3)*p/(p+2)]
    1×2 Matrix{RationalFunction{Int64, :s, Polynomial{Int64, :s}}}:
    (1 + 4*s + 4*s^2) // (6 + 10*s + 4*s^2)  (1 + 4*s + 4*s^2) // (12 + 14*s + 4*s^2)
  4. In the prototype implementation of the RationalTransferFunction type, the operations == and \isapprox must compare also the sampling times.
  5. The folowing is bug
    julia> r=Polynomial(0)/Polynomial(1)
    (0) // (1)
    julia> r(10)
    NaN
  6. The following should probably work:
    julia> p = Polynomial([1,2],:s)
    Polynomial(1 + 2*s)
    julia> 1/p
    ERROR: MethodError: no method matching /(::Int64, ::Polynomial{Int64, :s})
  7. The following constructions should probably work:
    julia> s = Polynomial([0, 1],'s');
    julia> G = [s^2 s/(s+1); 0 1/s]
    ERROR: MethodError: no method matching /(::Int64, ::Polynomial{Int64, :s})
    julia> G = [s^2 s/(s+1); im one(s)/s]
    ERROR: InexactError: Int64(0 + 1im)

    The type should be Float64 I guess:

    julia> G = [s^2 s/(s+1); 0. one(s)/s]
    2×2 Matrix{RationalFunction{Int64, :s, Polynomial{Int64, :s}}}:
    (s^2) // (1)  (s) // (1 + s)
    (0) // (1)    (1) // (s)

Finally, I would like to thank you for the very valuable work implementing the rational function type.

jverzani commented 3 years ago

Awesome feedback. I might not get to this until next week, but very much appreciated.

Let me do a "hot take":

  1. yes, that is a good idea
  2. oops, no reason that shouldn't work
  3. hmm, that one is annoying. I really like that destructuring idea, but maybe it will have to go
  4. Yes, that effort was only partial. It would likely make more sense for that special case to live in your package where you can tweak things
  5. Huh, definitely a bug
  6. Yeah 1//p`` (maybe not1/p`) should work (When I first wrote the type I had thought of a different package for it in which case that would have been type piracy. )
  7. Oops, promotion of the element type isn't right. Hopefully that can be adjusted easily enough.
jverzani commented 3 years ago

Thanks again! #339 should address all of these except for 4 (as I'm not sure where that should live). One difference is that // is needed to construct rational functions, not just /.