inputd.m

一个时滞系统的工具箱

function [T,V]=inputd(tf,a0,a1,b,c,d,u,delay,F,tol,IU);
%INPUTD   Simulation of continuous linear systems with delays to
%         arbitrary continuous input.
%	  INPUTD 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 the continuous perturbation function u(t) applied 
%	  to the single input.
%
%	  The time vector is automatically determined.  
%
%         When invoked with left hand arguments,
%            [T,V]=inputd(tf,a0,a1,b,c,d,u,delay);
%            [T,V]=inputd(tf,a0,a1,b,c,d,u,delay,F);
%            [T,V]=inputd(tf,a0,a1,b,c,d,u,delay,F,tol);
%            [T,V]=inputd(tf,a0,a1,b,c,d,u,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 INPUTD produces the graph.
%
%	  See also: STEPD, INITIALD, IMPULSED.

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

t = 0;
u1=eval(u);

if (nargin==11)
  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;
if size(u1,1)~=size(b,2)
  error('The length of vector u must be equal to the number of columns of matrix B');
end;
if size(u1,2)>1
  error('Control u must be column-vector'); end;
error(a0a1bcd(a0,a1,b,c,d));


% Separate the IU-th component of control
if nargin ~= 11
  IU = 1;
else
  punct = find(u == ';');
  temp1 = find(u == '[');
  temp2 = find(u == ']');
  punct = [temp1(1), punct, temp2(length(temp2))];
  u = u(punct(IU)+1:punct(IU+1)-1);
end;


% Compute the time vector and the corresponding state values
initfunt = mat2str(zeros(size(a0,2),1));
delayt = mat2str(delay);
a0t = mat2str(a0);
a1t = mat2str(a1);
bt = ['(' mat2str(b(:,IU)) ')' '*' '(' u ')'];
if nargin==8
  [tout,xout]=dde45lin(0,tf,zeros(size(a0,2),1),...
                       initfunt,delayt,a0t,a1t,bt);
end;
if nargin==9
  Ft = mat2str(F);
  [tout,xout]=dde45lin(0,tf,zeros(size(a0,2),1),...
                       initfunt,delayt,a0t,a1t,bt,delayt,Ft);
end;
if (nargin==10)|(nargin==11)
  Ft = mat2str(F);
  [tout,xout]=dde45lin(0,tf,zeros(size(a0,2),1),...
                       initfunt,delayt,a0t,a1t,bt,delayt,Ft,tol);
end;

% Compute the output 
u2=zeros(length(tout),1);
for i=1:length(tout)
  t=tout(i);
  u2(i)=eval(u);
end;
vout=xout*c'+u2*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;