dde45lin.m

一个时滞系统的工具箱

function [T,X]=dde45lin(t0,tf,x0,initfun,delay,A0,AA,B,delayF,F,tol)
%DDE45LIN Solves systems of differential equations with 
%	 delay of kind
%
%	 x`(t)=A0(t)*x(t)+summa(i=1,k,A_i(t)*y(-tau_i(t)))+
%	       integral(-tau_{k+1}(t),0,F(t,s)*y(s))ds+B(t).
%
%	 DDE45LIN integrates a system of functional-differential 
%	 equations using 4th and 5th order Runge-Kutta formulas
%	 as well interpolation by degenerate cubic splines.
%	 [T,X] = dde45lin(t0,tf,x0,initfun,delay,A0,AA,B,delayF,F) 
%	 integrates the system of functional-differential equations 
%	 over the segment from t0 to tf, with the initial value 
%	 x0 = x(t0) and the initial function y0(s) = x(t0+s) for s<0.
%	 [T,X] = dde45lin(t0,tf,x0,initfun,delay,A0,AA,B,delayF,F,tol) 
%	 uses tolerance tol.
%	 [T,X] = dde45lin(t0,tf,x0,initfun,delay,A0,AA,B) and 
%	 [T,X] = dde45lin(t0,tf,x0,initfun,delay,A0,AA,B,tol) are used
%	 for the systems that don't contain integral.
%	 Invoked without left-hand arguments DDE45LIN produces the graph.
%
%	 INPUT:
%	 t0      - Initial value of t.
%	 tf      - Final value of t.
%	 x0      - Initial value column-vector.
%	 initfun - String variable with column-vector of initial 
%		        function depending on variable s.
%	 delay   - String variable with string-vector of delays 
%		        (tau_1,...,tau_k) depending on t and corresponding 
%	           to matrices (A_1,...,A_k). 
%	 A0      - String variable with matrix function A0 
%		        depending on t.
%	 AA	   - String variable with string-vector of matrices
%		        (A_1,...,A_k) depending on t.
%	 B	      - String variable with column-vector B
%		        depending on t.
%	 delayF  - String variable with function tau_{k+1} 
%		        depending on variable t.
%	 F       - String variable with matrix F(t,s).
%	 tol     - The desired accuracy. (by default: tol = 1.e-4).
%
%	 OUTPUT:
%	 T  - Returned integration time points (column-vector).
%	 X  - Returned solution, one solution vector per T-value.
%
%	 The result can be displayed by: plot(T,X).
%
%	 See also dde45.



% Check for input arguments
if tf <= t0 
  error('The final point of integration must be greater than initial one'); 
end;
t=t0; s=t0;
[row,col] = size(x0);
if (col ~= 1)
   error('Initial state x0 must be column-vector'); 
end;
if size(eval(A0)) ~= [row, row] 
   error('Matrix A0 must have as many columns and strings as initial state'); 
end;
delayt = eval(delay);
if size(delayt,1) ~= 1 
   error('Fifth argument "delay" must be string-vector'); 
end;
if size(eval(AA)) ~= [row, row*length(eval(delay))] 
   error('Matrix AA must have length(x0) rows and length(x0)*length(delay) columns');
end;
if size(eval(initfun)) ~= [row, 1] 
   error('Dimensions of initial function and initial state must coincide'); 
end;
if size(eval(B)) ~= [row, 1]
   error('Dimensions of vector B and initial state must coincide');
end;
if nargin >= 10
   if length(eval(delayF)) ~= 1  
     error('Lower limit of integration "delayF" must be scalar'); 
   end;
   if size(eval(F)) ~= [row, row] 
     error('The A0 and F matrices must have the same dimensions.'); 
   end;
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)|(nargin == 8), tol = 1.e-4; end;
if nargin == 9  tol = delayF; end;

% Initialization
t1 = t0;
hmax = (tf - t1)/16;
h = hmax/8;
x = x0(:);
f = zeros(row,6);
chunk = 128;
tout = zeros(chunk,1);
xout = zeros(chunk,row);
k = 1;
tout(k) = t1;
xout(k,:) = x.';
len = length(delayt);
ydelayt = zeros(row*len,1);
integral=zeros([row,1]);


% The main loop

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

  % Compute the slopes
  t=t1;
  delayt=eval(delay); 
  A0t=eval(A0); 
  AAt=eval(AA); 
  Bt=eval(B);
  for i=1:len 
      ydelayt((i-1)*row+1:i*row)=ydelay(-delayt(i),t,tout(1:k),xout(1:k,:),initfun);
  end;
  if nargin > 9
      delayFt=eval(delayF);
      integral=int1(-delayFt,0,t,tout(1:k),xout(1:k,:),initfun,F);
  end;
  f(:,1) = A0t*x+AAt*ydelayt+Bt+integral;
  for j = 1:5
      t=t1+alpha(j)*h;
      delayt=eval(delay);
      A0t=eval(A0);
      AAt=eval(AA); 
      Bt=eval(B);
      for i=1:len 
         ydelayt((i-1)*row+1:i*row)=ydelay(-delayt(i),t,tout(1:k),xout(1:k,:),initfun);
      end;
      if nargin > 9
         delayFt=eval(delayF);
         integral=int1(-delayFt,0,t,tout(1:k),xout(1:k,:),initfun,F);
      end;
      f(:,j+1) = A0t*(x+h*f*beta(:,j))+AAt*ydelayt+Bt+integral;
   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
      t1 = t1 + h;
      x = x + h*f*gamma(:,1);
      k = k+1;
      if k > length(tout)
         tout = [tout; zeros(chunk,1)];
         xout = [xout; zeros(chunk,row)];
      end
      tout(k) = t1;
      xout(k,:) = x.';
   end
   % Update the step size
   if delta ~= 0.0
      h = min(hmax, 0.8*h*(tau/delta)^pow);
   end
end

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

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

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