octave_packages/m/polynomial/polyaffine.m

   1 ## Copyright (C) 2009-2012 Tony Richardson, Jaroslav Hajek
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn {Function File} {} polyaffine (@var{f}, @var{mu})
  21 ## Return the coefficients of the polynomial vector @var{f} after an affine
  22 ## transformation.  If @var{f} is the vector representing the polynomial f(x),
  23 ## then @code{@var{g} = polyaffine (@var{f}, @var{mu})} is the vector
  24 ## representing:
  25 ##
  26 ## @example
  27 ## g(x) = f( (x - @var{mu}(1)) / @var{mu}(2) )
  28 ## @end example
  29 ##
  30 ## @seealso{polyval, polyfit}
  31 ## @end deftypefn
  32
  33
  34 function g = polyaffine (f, mu)
  35
  36    if (nargin != 2)
  37       print_usage ();
  38    endif
  39
  40    if (! isvector (f))
  41       error ("polyaffine: F must be a vector");
  42    endif
  43
  44    if (! isvector (mu) || length (mu) != 2)
  45       error ("polyaffine: MU must be a two-element vector");
  46    endif
  47
  48    lf = length (f);
  49
  50    ## Ensure that f is a row vector
  51    if (rows (f) > 1)
  52       f = f.';
  53    endif
  54
  55    g = f;
  56
  57    ## Scale.
  58    if (mu(2) != 1)
  59      g = g ./ (mu(2) .^ (lf-1:-1:0));
  60    endif
  61
  62    ## Translate.
  63    if (mu(1) != 0)
  64      w = (-mu(1)) .^ (0:lf-1);
  65      ii = lf:-1:1;
  66      g = g(ii) * (toeplitz (w) .* pascal (lf, -1));
  67      g = g(ii);
  68    endif
  69
  70 endfunction
  71
  72
  73 %!demo
  74 %! f = [1/5 4/5 -7/5 -2];
  75 %! g = polyaffine (f, [1, 1.2]);
  76 %! x = linspace (-4, 4, 100);
  77 %! plot(x, polyval (f, x), x, polyval (g, x));
  78 %! legend ("original", "affine");
  79 %! axis ([-4 4 -3 5]);
  80 %! grid ("on");
  81
  82 %!test
  83 %! f = [1/5 4/5 -7/5 -2];
  84 %! mu = [1, 1.2];
  85 %! g = polyaffine (f, mu);
  86 %! x = linspace (-4, 4, 100);
  87 %! assert (polyval (f, x, [], mu), polyval (g, x), 1e-10);
  88