lv45.m

一个时滞系统的工具箱

function [T,X]=lv45(t0,tf,x0,initfun,delay,eps1,gamma1,gamma2,F1,F2,tol)
%LV45  Solves the Lotke-Volterra systems with delay of kind  
%
%      x1`(t)=[ eps1-gamma1*x2(t)-integral(-tau,0,F1(s)*y2(s))ds]*x1(t)
%      x2`(t)=[-eps1+gamma2*x1(t)+integral(-tau,0,F2(s)*y1(s))ds]*x2(t)
%
%      LV45 integrates a Lotke-Volterra system of functional-differential 
%      equations using 4th and 5th order Runge-Kutta formulas
%      as well interpolation by degenerate cubic splines.
%      [T,X] = lv45(t0,tf,x0,initfun,delay,eps1,gamma1,gamma2,F1,F2) 
%      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] = lv45(t0,tf,x0,initfun,delay,eps1,gamma1,gamma2,F1,F2,tol) 
%      uses tolerance tol.
%      Invoked without left-hand arguments LV45 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. 
%      eps1    - Scalar parameter.
%      gamma1  - Scalar parameter.
%      gamma2  - Scalar parameter.
%      F1,F2   - String variables with scalar functions F1(s), F2(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 input parameters
if tf <= t0 
  error('The final point of integration must be greater than initial one'); 
end;
if delay < 0
  error('Delay must be nonnegative constant');
end;
if (size(x0)~=[2, 1])
  error('Initial state x0 must be two-element column-vector'); 
end;
s=t0;
if (size(eval(initfun))~=[2, 1]) 
  error('Initial function must be two-element column-vector'); 
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 <11, tol = 1.e-4; end


% Initialization
t1 = t0;
hmax = (tf - t1)/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) = t1;
xout(k,:) = x.';

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

  % Compute the slopes
  f(:,1)=[(eps1-gamma1*x(2)-...
          int1(-delay,0,t1,tout(1:k),xout(1:k,:),initfun,F1,2))*x(1);
          (-eps1+gamma2*x(1)+...
          int1(-delay,0,t1,tout(1:k),xout(1:k,:),initfun,F2,1))*x(2)];
  for j = 1:5
    t=t1+alpha(j)*h;
    x1=x+h*f*beta(:,j);
    f(:,j+1)=[(eps1-gamma1*x1(2)-...
              int1(-delay,0,t,tout(1:k),xout(1:k,:),initfun,F1,2))*x1(1);
              (-eps1+gamma2*x1(1)+...
              int1(-delay,0,t,tout(1:k),xout(1:k,:),initfun,F2,1))*x1(2)];
  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,length(x))];
      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,:);