123.zip_matlab 时滞_时滞_时滞MATLAB_时滞微分_时滞微分方程

function sol = BM_dde23_wtl(ddefun,lags,history,tspan,options,varargin) %DDE23 Solve delay differential equations (DDEs) with constant delays. % SOL = DDE23(DDEFUN,LAGS,HISTORY,TSPAN) integrates a system of DDEs % y'(t) = f(t,y(t),y(t - tau_1),...,y(t - tau_k)). The constant, positive delays % tau_1,...,tau_k are input as the vector LAGS. The function DDEFUN(T,Y,Z) must % return a column vector corresponding to f(t,y(t),y(t - tau_1),...,y(t - tau_k)). % In the call to DDEFUN, T is the current t, the column vector Y approximates y(t), % and Z(:,j) approximates y(t - tau_j) for delay tau_j = LAGS(J). The DDEs % are integrated from T0 to TF where T0 < TF and TSPAN = [T0 TF]. % The solution at t <= T0 is specified by HISTORY in one of three ways: HISTORY % can be a function of t that returns the column vector y(t). If y(t) is constant, % HISTORY can be this column vector. If this call to DDE23 continues a previous % integration to T0, HISTORY is the solution SOL from that call. % % DDE23 produces a solution that is continuous on [T0,TF]. The solution is % evaluated at points TINT using the output SOL of DDE23 and the function % DEVAL: YINT = DEVAL(SOL,TINT). The output SOL is a structure with % SOL.x -- mesh selected by DDE23 % SOL.y -- approximation to y(t) at the mesh points of SOL.x % SOL.yp -- approximation to y'(t) at the mesh points of SOL.x % SOL.solver -- 'dde23' % % SOL = DDE23(DDEFUN,LAGS,HISTORY,TSPAN,OPTIONS) solves as above with default % parameters replaced by values in OPTIONS, a structure created with the % DDESET function. See DDESET for details. Commonly used options are % scalar relative error tolerance 'RelTol' (1e-3 by default) and vector of % absolute error tolerances 'AbsTol' (all components 1e-6 by default). % % SOL = DDE23(DDEFUN,LAGS,HISTORY,TSPAN,OPTIONS,P1,P2,...) passes the % additional parameters P1,P2,... to the DDE function as DDEFUN(T,Y,Z,P1,P2,...), % to the history (if it is a function) as HISTORY(T,P1,P2,...), and to % all functions specified in OPTIONS. Use OPTIONS = [] as a place holder % if no options are set. % % DDE23 can solve problems with discontinuities in the solution prior to T0 % (the history) or discontinuities in coefficients of the equations at known % values of t after T0 if the locations of these discontinuites are % provided in a vector as the value of the 'Jumps' option. % % By default the initial value of the solution is the value returned by % HISTORY at T0. A different initial value can be supplied as the value of % the 'InitialY' property. % % With the 'Events' property in OPTIONS set to a function EVENTS, DDE23 % solves as above while also finding where event functions % g(t,y(t),y(t - tau_1),...,y(t - tau_k)) are zero. For each function % you specify whether the integration is to terminate at a zero and whether % the direction of the zero crossing matters. These are the three vectors % returned by EVENTS: [VALUE,ISTERMINAL,DIRECTION] = EVENTS(T,Y,Z). % For the I-th event function: VALUE(I) is the value of the function, % ISTERMINAL(I) = 1 if the integration is to terminate at a zero of this % event function and 0 otherwise. DIRECTION(I) = 0 if all zeros are to % be computed (the default), +1 if only zeros where the event function is % increasing, and -1 if only zeros where the event function is decreasing. % The field SOL.xe is a column vector of times at which events occur. Rows % of SOL.ye are the corresponding solutions, and indices in vector SOL.ie % specify which event occurred. % % Example % sol = dde23(@ddex1de,[1, 0.2],@ddex1hist,[0, 5]); % solves a DDE on the interval [0, 5] with lags 1 and 0.2 and delay % differential equations computed by the function ddex1de. The history % is evaluated for t <= 0 by the function ddex1hist. The solution is % evaluated at 100 equally spaced points in [0 5] % tint = linspace(0,5); % yint = deval(sol,tint); % and plotted with % plot(tint,yint); % DDEX1 shows how this problem can be coded using subfunctions. For % another example see DDEX2. % % Class support for inputs TSPAN, LAGS, HISTORY, and the result of DDEFUN(T,Y,Z): % float: double, single % % See also DDESET, DDEGET, DEVAL. % DDE23 tracks discontinuities and integrates with the explicit Runge-Kutta % (2,3) pair and interpolant of ODE23. It uses iteration to take steps % longer than the lags. % Details are to be found in Solving DDEs in MATLAB, L.F. Shampine and % S. Thompson, Applied Numerical Mathematics, 37 (2001). % Jacek Kierzenka, Lawrence F. Shampine and Skip Thompson % Copyright 1984-2004 The MathWorks, Inc. % $Revision: 1.5.4.4 $ $Date: 2004/04/16 22:05:21 $ solver_name = 'dde23'; % Check inputs if nargin < 5 options = []; if nargin < 4 error('MATLAB:dde23:NotEnoughInputs',... 'Not enough input arguments. See DDE23.'); end end % Stats nsteps = 0; nfailed = 0; nfevals = 0; t0 = tspan(1); tfinal = tspan(end); % Ignore all entries of tspan except first and last. if tfinal <= t0 error('MATLAB:dde23:TspandEndLTtspan1', 'Must have tspan(1) < tspan(end).') end sol.solver = solver_name; if isnumeric(history) temp = history; sol.history = history; elseif isstruct(history) if history.x(end) ~= t0 error('MATLAB:dde23:NotContinueFromHistoryEnd',... 'Must continue from last point reached.') end temp = history.y(:,end); sol.history = history.history; else temp = feval(history,t0,varargin{:}); sol.history = history; end y0 = temp(:); maxlevel = 4; initialy = ddeget(options,'InitialY',[],'fast'); if ~isempty(initialy) y0 = initialy(:); maxlevel = 5; end neq = length(y0); % If solving a DDE, locate potential discontinuities. We need to step to each of % the points of potential lack of smoothness. Because we start at t0, we can % remove it from discont. The solver always steps to tfinal, so it is % convenient to add it to discont. if isempty(lags) discont = tfinal; minlag = Inf; else lags = lags(:)'; minlag = min(lags); if minlag <= 0 error('MATLAB:dde23:NotPosLags', 'The lags must all be positive.') end vl = t0; maxlag = max(lags); if isstruct(history) indices = find( history.discont < (t0 - maxlag) ); if ~isempty(indices) ndex = indices(end); sol.discont = history.discont(1:ndex); vl = [history.discont(ndex+1:end) t0]; end end jumps = ddeget(options,'Jumps',[],'fast'); if ~isempty(jumps) indices = find( ((t0 - maxlag) <= jumps) & (jumps <= tfinal) ); if ~isempty(indices) jumps = jumps(indices); vl = sort([vl jumps(:)']); maxlevel = 5; end end discont = vl; for level = 2:maxlevel vlp1 = vl(1) + lags; for i = 2:length(vl) vlp1 = [vlp1 (vl(i)+lags)]; end % Restrict to tspan. indices = find(vlp1 <= tfinal); vl = vlp1(indices); if isempty(vl) break; end nvl = length(vl); if nvl > 1 % Purge duplicates in vl. vl = sort(vl); indices = find(abs(diff(vl)) <= 10*eps(vl(1:nvl-1))) + 1; vl(indices) = []; end discont = [discont vl]; end if length(discont) > 1 discont = sort(discont); % Purge duplicates. indices = find(abs(diff(discont)) <= 10*eps(discont(1:end-1))) + 1; discont(indices) = []; end end % if isstruct(history) % sol.discont = [history.discont discont]; % else % sol.discont = discont; % end sol.discont = discont; % Add tfinal to the list of discontinuities if it is not already included. This % is a programming convenience and is not added to sol.discont. if abs(tfinal - discont(end)) <= 10*eps(tfi