viernes, 8 de octubre de 2010

Aproximación de Padé en SciLab

function P = Pade(q,n)
  // Aproximacion de Pade de orden «n» para «exp(-q.s)»
  // Calculado como:
  //      $\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$
  s=poly(0,'s');
  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

Ahora sí, corregido.

No hay comentarios:

Publicar un comentario