[Scilab-users] Simplify with threshold

[Federico Miyara]

Dear All,

Just in case anyone finds it useful, I'm attaching a function to 
simplify rationals which allows, unlike the native function simp, to set 
a threshold below which a zero and a pole can be cancelled.

Regards,

Federico Miyara


-- 
El software de antivirus Avast ha analizado este correo electrónico en busca de virus.
https://www.avast.com/antivirus
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/users/attachments/20200827/2602ca94/attachment.htm>
-------------- next part --------------
function [Hsimp, p, z] = simplify(H, th)
    // Simplify rational function when poles and zeros are 
    // found whose values differ less than a given numerical
    // threshold
    // 
    // Usage:
    //       [Hsimp, p, z] = simplify(H, th)
    // where
    //       H:     Rational function
    //       th:    Threshold; if omitted its default value is 1e-6 
    //       Hsimp: Version of H with coincident poles and zeros
    //              removed
    //       p:     Poles after simplification 
    //       z:     Zeros after simplification
    //
    // NOTE: The simplification is imperfect if poles and zeros
    //       are not identical, but the rational function computed 
    //       at specific values of the variable will differ very 
    //       little if the threshold is carefully selected. The
    //       threshold should be small, otherwise poles and zeros
    //       which differ considerably might pe cancelled out. 
    //
    // ----------------------------
    // Author: Federico Miyara
    // Date:   2020-08-27
    
    // Input argument error handling
    if argn(2)<1 | argn(2)>2
        error(msprintf(gettext(..
        "%s: Wrong number of input arguments: %d or %d expected.\n"),..
        "simplify",1,2)); 
    end
    if argn(2)<2
        th = 0.000001;
    end
    if type(th)<>1 | (type(th)==1 & ~isreal(th))
        error(msprintf(gettext(..
        "%s: Argument #%d: Decimal number expected.\n"), ..
        "simplify", 2))
    end
    if type(H)<>16 
        error(msprintf(gettext(..
        "%s: Argument #%d: Rational fraction expected.\n"), ..
        "simplify", 1))
    end
    if argn(1)>3
        error(msprintf(gettext(..
        "%s: Wrong number of output arguments: %d to %d expected.\n"),..
        "simplify",0,3));     
    end
   
    // Zeros and poles of H
    z = roots(H.num);
    p = roots(H.den);
    
    // Compute the differences between all poles and all 
    // zeros. Find the pole and zero whose difference is 
    // minimum. If that difference is smaller than the 
    // threshold, remove that pole anz zero. Iterate 
    // until no pole and zero differing less than the
    // threshold are left.
    terminate = 0;
    while terminate==0
        // Number of zeros and poles
        nz = length(z);
        np = length(p);
        // Initialize matrix of differences 
        delta = zeros(np,nz);
        // Compute matrix of differences
        for k=1:nz
            for n=1:np
                delta(n,k) = abs(z(k) - p(n));
            end 
        end
        // Find the minimum differenc
        [mini, I] = min(delta);
        // If the minimum is smaller than the threshold
        // remove the corresponding pole and zero.
        // Otherwise, terminate the iteration.
        if mini<th
            p(I(1)) = []; 
            z(I(2)) = [];
        else
            terminate = 1;
        end
    end 
    
    // Recreate simplified numerator and denominator
    Nsimp = poly(z, varn(N));
    Dsimp = poly(p, varn(D));
    
    // Simplified version of rational
    Hsimp = Nsimp/Dsimp;
     
endfunction