Mercurial > octave-libtiff
view scripts/polynomial/padecoef.m @ 23219:3ac9f9ecfae5 stable
maint: Update copyright dates.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Wed, 22 Feb 2017 12:39:29 -0500 |
parents | e9a0469dedd9 |
children | 092078913d54 |
line wrap: on
line source
## Copyright (C) 2014-2017 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 @nospell{@code{exp (-sT)}} is defined by the ## following equation ## ## @example ## @group ## Pn(s) ## exp (-sT) ~ ------- ## Qn(s) ## @end group ## @end example ## ## Where both @nospell{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 <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)