stepd.m

一个时滞系统的工具箱

function [T,V]=stepd(tf,a0,a1,b,c,d,delay,F,tol,IU);
%STEPD    Step response of continuous-time linear systems with delays.
%	  STEPD simulates the linear time-delay system
%	.
%	x(t) = A0*x(t)+A1*y(-delay)+int(-delay,0,F(s)*y(s))ds+B*u(t)
%		v(t) = C*x(t) + D*u(t)
%
%	  on the segment [0, tf] with the initial conditions
%         x(0)=[0], x(t)=[0], when t<0, 
%	  where [0] -- column-vector of zeros,
%         to an step applied to the single input.
%	  It's equivalent to the system
%	  .
%	  x(t) = A0*x(t)+A1*y(-delay)+int(-delay,0,F(s)*y(s))ds+B_IU
%		v(t) = C*x(t) + D_IU
%
%	  with the initial conditions x(0)=[0], x(t)=[0], when t<0.
%	  If parameter IU is given then B_IU and D_IU are IU-th 
%	  columns of B and D correspondingly, otherwise B_IU=B, D_IU=D.
%  
%	  The time vector is automatically determined.  
%
%         When invoked with left hand arguments,
%            [T,V]=stepd(tf,a0,a1,b,c,d,delay);
%            [T,V]=stepd(tf,a0,a1,b,c,d,delay,F);
%            [T,V]=stepd(tf,a0,a1,b,c,d,delay,F,tol);
%            [T,V]=stepd(tf,a0,a1,b,c,d,delay,F,tol,IU);
%	  returns the time history in the vector T and the output in
%         the matrix V. No plot is drawn on the screen. V has as many
%         columns as there are outputs and length(T) rows.
%	  Invoked without left hand arguments STEPD produces the graph.
%
%	  See also: IMPULSED, INITIALD, INPUTD.

% Check the input parameters
if tf < 0 
  error('First input argument must be positive.'); 
end;

if (nargin==10)
  if (IU<1)
    error('The number of input ID must be positive integer.'); 
  end;
  if (size(b,2)<IU)
    error('The number of input ID must not exceed the number of inputs.'); 
  end;
end;
error(a0a1bcd(a0,a1,b,c,d));

% Determine the IU-th component
if nargin < 10
   IU = 1;
end;

% Compute the time vector and the corresponding state values
initfunt = mat2str(zeros(size(a0,2),1));
if nargin==7
  [tout,xout]=dde45l1(0,tf,zeros(size(a0,2),1),...
                       initfunt,delay,a0,a1,b(:,IU));
end;
if nargin==8
  [tout,xout]=dde45l1(0,tf,zeros(size(a0,2),1),...
                       initfunt,delay,a0,a1,b(:,IU),delay,F);
end;
if (nargin==9)|(nargin==10)
  [tout,xout]=dde45l1(0,tf,zeros(size(a0,2),1),...
                       initfunt,delay,a0,a1,b(:,IU),delay,F,tol);
end;

% Compute the output 
vout=xout*c'+ones(length(tout),1)*d(:,IU)';

if nargout==0          % If no output arguments, plot graph
  plot(tout,vout)
  xlabel('Time (secs), t'), ylabel('Amplitude, v(t)')
  return    % Suppress output
end

T=tout;   
V=vout;