clsim.m

一个时滞系统的工具箱

function [T, X]=clsim(tf,x0,initfun,delay,A0,A1,B,C,D,tol)
%CLSIM   solves closed-loop system of differential equations with delay 
%
%        x`(t)=(A0+B*C)*x(t)+A1*x(t-delay)+B*int(-delay,0,D(s)*x(t+s))ds
%
%			corresponding to the linear system with delay
%
%			x`(t)=A0*x(t)+A1*x(t-delay)+B*u
%
%			and the feedback control
%
%			u=C*x(t)+int(-delay,0,D(s)*x(t+s))ds.
%
%      	CLSIM uses 4th and 5th order Runge-Kutta formulas
%			as well interpolation by degenerate cubic splines.
%
%	 [T,X] = clsim(tf,x0,initfun,delay,A0,A1,B,C,D) 
%	 integrates the system of functional-differential equations 
%	 over the segment from 0 to tf, with the initial value 
%	 x0 = x(0) and the initial function y0(s) = x(s) for s<0.
%	 [T,X] = clsim(tf,x0,initfun,delay,A0,A1,B,C,D,tol) 
%	 uses tolerance tol.
%	 Invoked without left-hand arguments CLSIM produces the graph.
%
%	INPUT:
%	tf        - Final value of t.
%	x0        - Initial value column-vector.
%  initfun   - String variable with column-vector of initial function 
%					depending on variable s.
%	delay     - A constant positive delay. 
%  A0,A1,B,C - Matrices of constant coefficients.
%  D         - String variable with matrix D depending on s.
%	tol       - The desired accuracy. (by default: tol = 1.e-4).
%
%	OUTPUT:
%	T      - returned integration time points (column-vector).
%	X      - returned solution, one solution column-vector per T-value.
%
%	The result can be displayed by: plot(T, X).
%
%	See also: clso.



% Check the integration segment
if tf <= 0 
  error('The final point of integration must be greater than zero'); 
end;

% The Fehlberg coefficients:
alpha = [1/4  3/8  12/13  1  1/2]';
beta  = [ [    1      0      0     0      0    0]/4
          [    3      9      0     0      0    0]/32
          [ 1932  -7200   7296     0      0    0]/2197
          [ 8341 -32832  29440  -845      0    0]/4104
          [-6080  41040 -28352  9295  -5643    0]/20520 ]';
gamma = [ [902880  0  3953664  3855735  -1371249  277020]/7618050
          [ -2090  0    22528    21970    -15048  -27360]/752400 ]';
pow = 1/5;
if nargin < 10, tol = 1.e-4; end

% Initialization
t = 0;
hmax = (tf - t)/16;
h = hmax/8;
x = x0(:);
f = zeros(length(x),6);
chunk = 128;
tout = zeros(chunk,1);
xout = zeros(chunk,length(x));
k = 1;
tout(k) = t;
xout(k,:) = x.';
A = A0+B*C;

% The main loop

while (t < tf) & (t + h > t)
   if t + h > tf, h = tf - t; end

	% Compute the slopes
	temp1 = ydelay(-delay,t,tout(1:k),xout(1:k,:),initfun);
	temp2 = int1(-delay,0,t,tout(1:k),xout(1:k,:),initfun,D);
   f(:,1) = A*x + A1*temp1 + B*temp2;
   for j = 1:5
      temp3 = t+alpha(j)*h;
      temp1 = ydelay(-delay,temp3,tout(1:k),xout(1:k,:),initfun);
		temp2 = int1(-delay,0,temp3,tout(1:k),xout(1:k,:),initfun,D);
      f(:,j+1) = A*(x+h*f*beta(:,j)) + A1*temp1 + B*temp2;
   end
   
   % Estimate the error and the acceptable error
   delta = norm(h*f*gamma(:,2),'inf');
   tau = tol*max(norm(x,'inf'),1.0);
   
   % Update the solution only if the error is acceptable
   if delta <= tau
      t = t + h;
      x = x + h*f*gamma(:,1);
      k = k+1;
      if k > length(tout)
         tout = [tout; zeros(chunk,1)];
         xout = [xout; zeros(chunk,length(x))];
      end
      tout(k) = t;
      xout(k,:) = x.';

   end
   % Update the step size
   if delta ~= 0.0
      h = min(hmax, 0.8*h*(tau/delta)^pow);
   end
end

if (t < tf)
   disp('Singularity likely.')
   t
end

% Plot Graph
if nargout == 0
   plot(tout(1:k),xout(1:k,:));
   return;
end;

T = tout(1:k);
X = xout(1:k,:);