dde45l1.m

一个时滞系统的工具箱

function [T,X]=dde45l1(t0,tf,x0,initfun,delay,A0,AA,B,delayF,F,tol)
%DDE45L1  Solves systems of functional-differential equations with 
%	  delay of kind
%	  .
%	  x(t)=A0*x(t)+AA*y(-tau)+integral(-tau_F,0,F(t,s)*y(s))ds+B.
%
%	  DDE45L1 integrates a system of functional-differential 
%	  equations using 4th and 5th order Runge-Kutta formulas
%	  as well interpolation by degenerate cubic splines.
%	  [T,X] = dde45l1(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] = dde45l1(t0,tf,x0,initfun,delay,A0,AA,B,delayF,F,tol) 
%	  uses tolerance tol.
%	  [T,X] = dde45l1(t0,tf,x0,initfun,delay,A0,AA,B) and 
%	  [T,X] = dde45l1(t0,tf,x0,initfun,delay,A0,AA,B,tol) are used
%	  for the systems that don't contain integral.
%	  Invoked without left-hand arguments DDE45L1 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   - Constant delay tau. 
%	  A0, AA  - Constant matrices A0, AA.
%	  B	  - Column-vector B.  
%	  delayF  - Constant tau_F.
%	  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, dde1, dde2.



% 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(A0) ~= [row, row] 
   error('Matrix A0 must have as many columns and strings as initial state'); 
end;
if length(delay) ~= 1 
   error('Fifth argument "delay" must be scalar'); 
end;
if delay < 0
   error('Delay must be nonnegative');
end;
if size(AA) ~= [row, row] 
   error('Dimensions of matrices AA and A0 must coincide');
end;
if size(eval(initfun)) ~= [row, 1] 
   error('Dimensions of initial function and initial state must coincide'); 
end;
if size(B) ~= [row, 1]
   error('Dimensions of vector B and initial state must coincide');
end;
if nargin >= 10
   if length(delayF) ~= 1  
     error('Lower limit of integration "delayF" must be scalar'); 
   end;
   if delayF < 0
      error('Lower limit of integration must be nonnegative');
   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.';

% The main loop

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

  % Compute the slopes
  if nargin > 9
    f(:,1) = A0*x+AA*ydelay(-delay,t1,tout(1:k),xout(1:k,:),initfun)...
             +B+int1(-delayF,0,t1,tout(1:k),xout(1:k,:),initfun,F);
    for j = 1:5
      t=t1+alpha(j)*h;
      f(:,j+1) = A0*(x+h*f*beta(:,j))+B+...
                 int1(-delayF,0,t,tout(1:k),xout(1:k,:),initfun,F)+...
                 AA*ydelay(-delay,t,tout(1:k),xout(1:k,:),initfun);
    end
  else
    f(:,1) = A0*x+AA*ydelay(-delay,t1,tout(1:k),xout(1:k,:),initfun)+B;
    for j = 1:5
      t=t1+alpha(j)*h;
      f(:,j+1) = A0*(x+h*f*beta(:,j))+B+...
                 AA*ydelay(-delay,t,tout(1:k),xout(1:k,:),initfun);
    end
  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,:);