powerin.m

计算动力学系统的分岔图

function ans = powerin(x,y)
% X.^Y for adtayl objects.
% The case where X is a general adtayl object and Y is an integer.
switch y
  case 1, ans = x;
  case 2, x.*x;
  case 3, x.*x.*x;
  case 4, x.*x;  ans.*ans; 
  otherwise
    if ~(isnumeric(y) && fix(y)==y)
      error('X^Y for adtayl objects currently requires Y integer')
    end
    xtc = x.tc;
    if y>0
      pow = power0(xtc,y);
    else
      pow = zeros(size(xtc)); pow(1)=1; % vector [1,0,...,0], representing constant 1.
      if y<0
        pow = filter(1,power0(xtc,-y), pow);
      end
    end
    class(struct('tc',pow),'adtayl');
end

%---------------------------
function ans = power0(xtc,y)
% Assumes xtc is a vector and y>=1.
if y==1
  xtc;
else
 %yover2 = fix(y/2); yrem2 = rem(y,2);
  pyover2 = power0(xtc,fix(y/2));
  filter(pyover2,1, pyover2);
  if rem(y,2)~=0
    filter(xtc,1, ans);
  end
end