Mercurial > octave
changeset 21314:a4235a7eeb2c
Add new padecoef.m function.
* padecoef.m: New function.
* NEWS: Announce new function in 4.2 release.
* scripts/polynomial/module.mk: Add to build system.
* __unimplemented__.m: Remove from unimplemented lst.
* poly.txi: Add padecoef to manual.
| author | Endre Kozma <endre.kozma@gmx.com> |
|---|---|
| date | Thu, 11 Sep 2014 19:18:21 +0200 |
| parents | 4dda463ca576 |
| children | ea2c2e08ceda |
| files | NEWS doc/interpreter/poly.txi scripts/help/__unimplemented__.m scripts/polynomial/module.mk scripts/polynomial/padecoef.m |
| diffstat | 5 files changed, 185 insertions(+), 1 deletions(-) [+] |
line wrap: on
line diff
--- a/NEWS Sat Feb 20 12:08:24 2016 -0800 +++ b/NEWS Thu Sep 11 19:18:21 2014 +0200 @@ -73,6 +73,7 @@ odeset odeget ode45 + padecoef rad2deg ** Deprecated functions.
--- a/doc/interpreter/poly.txi Sat Feb 20 12:08:24 2016 -0800 +++ b/doc/interpreter/poly.txi Thu Sep 11 19:18:21 2014 +0200 @@ -349,6 +349,12 @@ @end float @end ifnotinfo +A very specific form of polynomial interpretation is the Pade' approximant. +For control systems, a continuous-time delay can be modeled very simply with +the approximant. + +@DOCSTRING(padecoef) + The function, @code{ppval}, evaluates the piecewise polynomials, created by @code{mkpp} or other means, and @code{unmkpp} returns detailed information about the piecewise polynomial.
--- a/scripts/help/__unimplemented__.m Sat Feb 20 12:08:24 2016 -0800 +++ b/scripts/help/__unimplemented__.m Thu Sep 11 19:18:21 2014 +0200 @@ -767,7 +767,6 @@ "ordeig", "ordqz", "outerjoin", - "padecoef", "parseSoapResponse", "partition", "pathtool",
--- a/scripts/polynomial/module.mk Sat Feb 20 12:08:24 2016 -0800 +++ b/scripts/polynomial/module.mk Thu Sep 11 19:18:21 2014 +0200 @@ -11,6 +11,7 @@ scripts/polynomial/deconv.m \ scripts/polynomial/mkpp.m \ scripts/polynomial/mpoles.m \ + scripts/polynomial/padecoef.m \ scripts/polynomial/pchip.m \ scripts/polynomial/poly.m \ scripts/polynomial/polyaffine.m \
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/polynomial/padecoef.m Thu Sep 11 19:18:21 2014 +0200 @@ -0,0 +1,177 @@ +## Copyright (C) 2014 Endre Kozma +## +## This file is part of Octave. +## +## Octave 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. +## +## Octave 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 Octave; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {} {[@var{num}, @var{den}] =} padecoef (@var{T}) +## @deftypefnx {} {[@var{num}, @var{den}] =} padecoef (@var{T}, @var{N}) +## Compute the @var{N}th-order Pad@'e approximant of the continuous-time +## delay @var{T} in transfer function form. +## +## @tex +## The Pad\'e approximant of $e^{-sT}$ is defined by the following equation +## $$ e^{-sT} \approx {P_n(s) \over Q_n(s)} $$ +## where both $P_n(s)$ and $Q_n(s)$ are $N^{th}$-order rational functions +## defined by the following expressions +## $$ P_n(s)=\sum_{k=0}^N {(2N - k)!N!\over (2N)!k!(N - k)!}(-sT)^k $$ +## $$ Q_n(s) = P_n(-s) $$ +## @end tex +## @ifnottex +## The Pad@'e approximant of exp (-sT) is defined by the following equation +## +## @example +## @group +## Pn(s) +## exp (-sT) ~ ------- +## Qn(s) +## @end group +## @end example +## +## Where both Pn(s) and Qn(s) are @var{N}th-order rational functions defined by +## the following expressions +## +## @example +## @group +## N (2N - k)!N! k +## Pn(s) = SUM --------------- (-sT) +## k=0 (2N)!k!(N - k)! +## +## Qn(s) = Pn(-s) +## @end group +## @end example +## +## @end ifnottex +## +## The inputs @var{T} and @var{N} must be non-negative numeric scalars. If +## @var{N} is unspecified it defaults to 1. +## +## The output row vectors @var{num} and @var{den} contain the numerator and +## denominator coefficients in descending powers of s. Both are +## @var{N}th-order polynomials. +## +## For example: +## +## @smallexample +## @group +## t = 0.1; +## n = 4; +## [num, den] = padecoef (t, n) +## @result{} num = +## +## 1.0000e-04 -2.0000e-02 1.8000e+00 -8.4000e+01 1.6800e+03 +## +## @result{} den = +## +## 1.0000e-04 2.0000e-02 1.8000e+00 8.4000e+01 1.6800e+03 +## @end group +## @end smallexample +## @end deftypefn + +function [num, den] = padecoef (T, N = 1) + + if (nargin < 1 || nargin > 2) + print_usage (); + endif + + if (! (isscalar (T) && isnumeric (T) && T >= 0)) + error ("padecoef: T must be a non-negative scalar"); + elseif (! (isscalar (N) && isnumeric (N) && N >= 0)) + error ("padecoef: N must be a non-negative scalar"); + endif + + N = round (N); + k = N : -1 : 0; + num = prod (linspace ((N - k + 1), (2 * N - k), N)', ones (1, N)) ... + / prod (N + 1 : 2 * N) ./ factorial (k); + num /= num(1); + den = num .* (T .^ k); + num .*= ((-T) .^ k); + +endfunction + + +%!test +%! T = 1; +%! [n_obs, d_obs] = padecoef (T); +%! n_exp = [1, 2] .* [-T, 1]; +%! d_exp = [1, 2] .* [T, 1]; +%! assert ([n_obs, d_obs], [n_exp, d_exp], eps); + +%!test +%! T = 0.1; +%! [n_obs, d_obs] = padecoef (T); +%! n_exp = [1, 2] .* [-T, 1]; +%! d_exp = [1, 2] .* [T, 1]; +%! assert ([n_obs, d_obs], [n_exp, d_exp], eps); + +%!test +%! T = 1; +%! N = 2; +%! k = N : -1 : 0; +%! [n_obs, d_obs] = padecoef (T, N); +%! n_exp = [1, 6, 12] .* ((-T) .^ k); +%! d_exp = [1, 6, 12] .* (T .^ k); +%! assert ([n_obs, d_obs], [n_exp, d_exp], eps); + +%!test +%! T = 0.25; +%! N = 2; +%! k = N : -1 : 0; +%! [n_obs, d_obs] = padecoef (T, 2); +%! n_exp = [1, 6, 12] .* ((-T) .^ k); +%! d_exp = [1, 6, 12] .* (T .^ k); +%! assert ([n_obs, d_obs], [n_exp, d_exp], eps); + +%!test +%! T = 0.47; +%! N = 3; +%! k = N : -1 : 0; +%! [n_obs, d_obs] = padecoef (T, N); +%! n_exp = [1, 12, 60, 120] .* ((-T) .^ k); +%! d_exp = [1, 12, 60, 120] .* (T .^ k); +%! assert ([n_obs, d_obs], [n_exp, d_exp], eps); + +%!test +%! T = 1; +%! N = 7; +%! i = 0 : 2 * N; +%! b = ((-T) .^ i) ./ factorial (i); +%! A = [[eye(N + 1); zeros(N, N + 1)], ... +%! [zeros(1, N); toeplitz(-b(1 : 2 * N), [-b(1), zeros(1, N-1)])]]; +%! x = A \ b'; +%! k = N : -1 : 0; +%! d_exp = [flipud(x(N + 2 : 2 * N + 1)); 1]'; +%! n_exp = flipud(x(1 : N + 1))'; +%! n_exp ./= d_exp(1); +%! d_exp ./= d_exp(1); +%! [n_obs, d_obs] = padecoef (T, N); +%! assert ([n_obs, d_obs], [n_exp, d_exp], 1e-2); + +## For checking in Wolfram Alpha (look at Alternate forms -> more): +## PadeApproximant[Exp[-x * T], {x, 0, {n, n}}] + +%% Test input validation +%!error padecoef () +%!error padecoef (1,2,3) + error ("padecoef: T must be a non-negative scalar"); +%!error <T must be a non-negative scalar> padecoef ([1,2]) +%!error <T must be a non-negative scalar> padecoef ({1}) +%!error <T must be a non-negative scalar> padecoef (-1) +%!error <N must be a non-negative scalar> padecoef (1, [1,2]) +%!error <N must be a non-negative scalar> padecoef (1, {1}) +%!error <N must be a non-negative scalar> padecoef (1, -1) +