dde45.m

一个时滞系统的工具箱

function [tout, xout]=dde45(xpfun,t0,tf,x0,initfun,tol)
%DDE45	Solves functional differential equations (systems with delays).
%	DDE45 integrates a system of functional-differential equations using
%	4th and 5th order Runge-Kutta formulas.
%	[tout,xout] = dde45(xpfun, t0, tf, x0, initfun) integrates the system
%  of functional-differential equations described by the M-file XPRIME.M
%	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.
%	[tout,xout] = dde45(xpfun, t0, tf, x0, initial, tol) uses tolerance tol.
%
%	INPUT:
%	xpfun   - String containing name of user-supplied problem description.
%	          Call: xprime = fun(t,x,tt,xt,initfun) where xpfun = 'fun'.
%	          t       - Time (scalar).
%	          x       - Solution column-vector.
%            tt      - Time history array.
%            xt      - Values x corresponding to array tt.
%            initfun - Initial function.
%	          xprime  - Returned derivative column-vector. 
%                      xprime(i) = dx(i)/dt.
%            (For more details see descriptions.)
%	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.
%	tol     - The desired accuracy. (by default: tol = 1.e-4).
%
%	OUTPUT:
%	tout  - Returned integration time points (column-vector).
%	xout  - Returned solution, one solution vector per tout-value.
%
%	The result can be displayed by: plot(tout, xout).
%
%	See also dde45lin.

% Check the integration segment
if tf <= t0 
  error('The final point of integration must be greater than initial one'); 
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 < 6, tol = 1.e-4; end

% Initialization
t = t0;
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.';

% The main loop

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

   % Compute the slopes
   temp = feval(xpfun,t,x,tout(1:k),xout(1:k,:),initfun);
   f(:,1) = temp(:);
   for j = 1:5
      temp = feval(xpfun,t+alpha(j)*h,x+h*f*beta(:,j),...
                           tout(1:k),xout(1:k,:),initfun);
      f(:,j+1) = temp(:);
   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

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