impulsed.m

一个时滞系统的工具箱

function [T,V]=impulsed(tf,a0,a1,b,c,d,delay,F,tol,IU);
%IMPULSED Impulse response of continuous-time linear systems with delays.
%	  IMPULSED 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 impulse 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
%		v(t) = C*x(t)
%
%	  with the initial conditions x(0)=B_IU, x(t)=[0], when t<0.
%	  If parameter IU is given then B_IU is IU-th column of B,
%	  otherwise B_IU=B.
% 
%	  The time vector is automatically determined.  
%
%         When invoked with left hand arguments,
%            [T,V]=impulsed(tf,a0,a1,b,c,d,delay);
%            [T,V]=impulsed(tf,a0,a1,b,c,d,delay,F);
%            [T,V]=impulsed(tf,a0,a1,b,c,d,delay,F,tol);
%            [T,V]=impulsed(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 IMPULSED produces the graph.
%
%	  See also: STEPD, INITIALD, INPUTD.


% Check the input parameters
if tf < 0 
  error('The first argument must be positive.'); 
end;
error(a0a1bcd(a0,a1,b,c,d));
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;
  b=b(:,IU);
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,b,initfunt,delay,a0,a1,...
                       zeros(length(b),1));
end;
if nargin==8
  [tout,xout]=dde45l1(0,tf,b,initfunt,delay,a0,a1,...
                       zeros(length(b),1),delay,F);
end;
if (nargin==9)|(nargin==10)
  [tout,xout]=dde45l1(0,tf,b,initfunt,delay,a0,a1,...
                       zeros(length(b),1),delay,F,tol);
end;


% Compute the output 
vout=xout*c';


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;