Aproximación de Padé de tiempos muertos en Octave

Padé approximation of delays for GNU Octave

--- pade.m ---

## -*- texinfo -*-
## @deftypefn  {Function File} {} pade (@var{q})
## @deftypefnx {Function File} {} pade (@var{q}, @var{n})
##
## Calcula la aproximacion de Pade de orden @var{n} para exp(-q.s) (s = variable compleja),
## y devuelve una funcion de transferencia continua.
##
## @var{q} es un numero entero.
##
## @var{n} es el orden de la aproximacion de Pade. Si se omite se asume @samp{"1"}.
##
## La aproximacion se calcula segun (OZBAY, Hitay. Introduction to Feedback Control Theory. USA: CRC Press, 1999. 232 p. ISBN 0-8493-1867-X):
##
##    $\textrm{e}^{-\Theta\,t}\approxeq\dfrac{\sum_{k=0}^{n}\,(-1)^{k}\, c_{k}\,(\theta s)^{k}}{\sum_{k=0}^{n}\, c_{k}\,(\theta s)^{k}}$
##
##    donde: $c_{k}=\dfrac{(2n-k)!\, n!}{(2n)!\, k!\,(n-k)!}$
##
##    para $k=1,2,\ldots,n$
##
## @end deftypefn

## Author: Alejandro Regodesebes -
## v.1.0 - Copyright (C) 2012 Alejandro Regodesebes

#    This program is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program.  If not, see .
 
function p = pade(q,n=1)
  if (not(isnumeric(q)))
    error("Los argumentos deben ser constantes enteras.");
  endif
  if (nargin == 0 || nargin > 2)
    print_usage ();
  endif
  s = tf("s"); % s es la variable de una función de transferencia
  Num = 0; % Inicializo el numerador
  Den = 0; % Inicializo el denominador
  nn = prod(1:n); % nn = n!
  n2 = prod(1:(2*n)); % n2 = (2n)!
  for k=0:n
    c = (prod(1:(2*n-k))*nn)/(n2*prod(1:k)*prod(1:(n-k)));
    Num = Num + ((-1)**k)*c*(q*s)**k;
    Den = Den + c*(q*s)**k;
  end;
  p = Num/Den;
endfunction

--- fin pade.m ---