1 %# Copyright (C) 2008-2012, Thomas Treichl <treichl@users.sourceforge.net>
2 %# OdePkg - A package for solving ordinary differential equations and more
4 %# This program is free software; you can redistribute it and/or modify
5 %# it under the terms of the GNU General Public License as published by
6 %# the Free Software Foundation; either version 2 of the License, or
7 %# (at your option) any later version.
9 %# This program is distributed in the hope that it will be useful,
10 %# but WITHOUT ANY WARRANTY; without even the implied warranty of
11 %# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 %# GNU General Public License for more details.
14 %# You should have received a copy of the GNU General Public License
15 %# along with this program; If not, see <http://www.gnu.org/licenses/>.
18 %# @deftypefn {Function File} {[@var{}] =} ode78d (@var{@@fun}, @var{slot}, @var{init}, @var{lags}, @var{hist}, [@var{opt}], [@var{par1}, @var{par2}, @dots{}])
19 %# @deftypefnx {Command} {[@var{sol}] =} ode78d (@var{@@fun}, @var{slot}, @var{init}, @var{lags}, @var{hist}, [@var{opt}], [@var{par1}, @var{par2}, @dots{}])
20 %# @deftypefnx {Command} {[@var{t}, @var{y}, [@var{xe}, @var{ye}, @var{ie}]] =} ode78d (@var{@@fun}, @var{slot}, @var{init}, @var{lags}, @var{hist}, [@var{opt}], [@var{par1}, @var{par2}, @dots{}])
22 %# This function file can be used to solve a set of non--stiff delay differential equations (non--stiff DDEs) with a modified version of the well known explicit Runge--Kutta method of order (7,8).
24 %# If this function is called with no return argument then plot the solution over time in a figure window while solving the set of DDEs that are defined in a function and specified by the function handle @var{@@fun}. The second input argument @var{slot} is a double vector that defines the time slot, @var{init} is a double vector that defines the initial values of the states, @var{lags} is a double vector that describes the lags of time, @var{hist} is a double matrix and describes the history of the DDEs, @var{opt} can optionally be a structure array that keeps the options created with the command @command{odeset} and @var{par1}, @var{par2}, @dots{} can optionally be other input arguments of any type that have to be passed to the function defined by @var{@@fun}.
26 %# In other words, this function will solve a problem of the form
28 %# dy/dt = fun (t, y(t), y(t-lags(1), y(t-lags(2), @dots{})))
30 %# y(slot(1)-lags(1)) = hist(1), y(slot(1)-lags(2)) = hist(2), @dots{}
33 %# If this function is called with one return argument then return the solution @var{sol} of type structure array after solving the set of DDEs. The solution @var{sol} has the fields @var{x} of type double column vector for the steps chosen by the solver, @var{y} of type double column vector for the solutions at each time step of @var{x}, @var{solver} of type string for the solver name and optionally the extended time stamp information @var{xe}, the extended solution information @var{ye} and the extended index information @var{ie} all of type double column vector that keep the informations of the event function if an event function handle is set in the option argument @var{opt}.
35 %# If this function is called with more than one return argument then return the time stamps @var{t}, the solution values @var{y} and optionally the extended time stamp information @var{xe}, the extended solution information @var{ye} and the extended index information @var{ie} all of type double column vector.
40 %# the following code solves an anonymous implementation of a chaotic behavior
43 %# fcao = @@(vt, vy, vz) [2 * vz / (1 + vz^9.65) - vy];
45 %# vopt = odeset ("NormControl", "on", "RelTol", 1e-3);
46 %# vsol = ode78d (fcao, [0, 100], 0.5, 2, 0.5, vopt);
48 %# vlag = interp1 (vsol.x, vsol.y, vsol.x - 2);
49 %# plot (vsol.y, vlag); legend ("fcao (t,y,z)");
53 %# to solve the following problem with two delayed state variables
56 %# d y1(t)/dt = -y1(t)
57 %# d y2(t)/dt = -y2(t) + y1(t-5)
58 %# d y3(t)/dt = -y3(t) + y2(t-10)*y1(t-10)
61 %# one might do the following
64 %# function f = fun (t, y, yd)
65 %# f(1) = -y(1); %% y1' = -y1(t)
66 %# f(2) = -y(2) + yd(1,1); %% y2' = -y2(t) + y1(t-lags(1))
67 %# f(3) = -y(3) + yd(2,2)*yd(1,2); %% y3' = -y3(t) + y2(t-lags(2))*y1(t-lags(2))
70 %# res = ode78d (@@fun, T, [1;1;1], [5, 10], ones (3,2));
78 function [varargout] = ode78d (vfun, vslot, vinit, vlags, vhist, varargin)
80 if (nargin == 0) %# Check number and types of all input arguments
82 error ('OdePkg:InvalidArgument', ...
83 'Number of input arguments must be greater than zero');
88 elseif (~isa (vfun, 'function_handle'))
89 error ('OdePkg:InvalidArgument', ...
90 'First input argument must be a valid function handle');
92 elseif (~isvector (vslot) || length (vslot) < 2)
93 error ('OdePkg:InvalidArgument', ...
94 'Second input argument must be a valid vector');
96 elseif (~isvector (vinit) || ~isnumeric (vinit))
97 error ('OdePkg:InvalidArgument', ...
98 'Third input argument must be a valid numerical value');
100 elseif (~isvector (vlags) || ~isnumeric (vlags))
101 error ('OdePkg:InvalidArgument', ...
102 'Fourth input argument must be a valid numerical value');
104 elseif ~(isnumeric (vhist) || isa (vhist, 'function_handle'))
105 error ('OdePkg:InvalidArgument', ...
106 'Fifth input argument must either be numeric or a function handle');
110 if (~isstruct (varargin{1}))
111 %# varargin{1:len} are parameters for vfun
112 vodeoptions = odeset;
113 vfunarguments = varargin;
115 elseif (length (varargin) > 1)
116 %# varargin{1} is an OdePkg options structure vopt
117 vodeoptions = odepkg_structure_check (varargin{1}, 'ode78d');
118 vfunarguments = {varargin{2:length(varargin)}};
120 else %# if (isstruct (varargin{1}))
121 vodeoptions = odepkg_structure_check (varargin{1}, 'ode78d');
126 else %# if (nargin == 5)
127 vodeoptions = odeset;
131 %# Start preprocessing, have a look which options have been set in
132 %# vodeoptions. Check if an invalid or unused option has been set and
134 vslot = vslot(:)'; %# Create a row vector
135 vinit = vinit(:)'; %# Create a row vector
136 vlags = vlags(:)'; %# Create a row vector
138 %# Check if the user has given fixed points of time
139 if (length (vslot) > 2), vstepsizegiven = true; %# Step size checking
140 else vstepsizegiven = false; end
142 %# Get the default options that can be set with 'odeset' temporarily
145 %# Implementation of the option RelTol has been finished. This option
146 %# can be set by the user to another value than default value.
147 if (isempty (vodeoptions.RelTol) && ~vstepsizegiven)
148 vodeoptions.RelTol = 1e-6;
149 warning ('OdePkg:InvalidOption', ...
150 'Option "RelTol" not set, new value %f is used', vodeoptions.RelTol);
151 elseif (~isempty (vodeoptions.RelTol) && vstepsizegiven)
152 warning ('OdePkg:InvalidOption', ...
153 'Option "RelTol" will be ignored if fixed time stamps are given');
154 %# This implementation has been added to odepkg_structure_check.m
155 %# elseif (~isscalar (vodeoptions.RelTol) && ~vstepsizegiven)
156 %# error ('OdePkg:InvalidOption', ...
157 %# 'Option "RelTol" must be set to a scalar value for this solver');
160 %# Implementation of the option AbsTol has been finished. This option
161 %# can be set by the user to another value than default value.
162 if (isempty (vodeoptions.AbsTol) && ~vstepsizegiven)
163 vodeoptions.AbsTol = 1e-6;
164 warning ('OdePkg:InvalidOption', ...
165 'Option "AbsTol" not set, new value %f is used', vodeoptions.AbsTol);
166 elseif (~isempty (vodeoptions.AbsTol) && vstepsizegiven)
167 warning ('OdePkg:InvalidOption', ...
168 'Option "AbsTol" will be ignored if fixed time stamps are given');
169 else %# create column vector
170 vodeoptions.AbsTol = vodeoptions.AbsTol(:);
173 %# Implementation of the option NormControl has been finished. This
174 %# option can be set by the user to another value than default value.
175 if (strcmp (vodeoptions.NormControl, 'on')), vnormcontrol = true;
176 else vnormcontrol = false;
179 %# Implementation of the option NonNegative has been finished. This
180 %# option can be set by the user to another value than default value.
181 if (~isempty (vodeoptions.NonNegative))
182 if (isempty (vodeoptions.Mass)), vhavenonnegative = true;
184 vhavenonnegative = false;
185 warning ('OdePkg:InvalidOption', ...
186 'Option "NonNegative" will be ignored if mass matrix is set');
188 else vhavenonnegative = false;
191 %# Implementation of the option OutputFcn has been finished. This
192 %# option can be set by the user to another value than default value.
193 if (isempty (vodeoptions.OutputFcn) && nargout == 0)
194 vodeoptions.OutputFcn = @odeplot;
195 vhaveoutputfunction = true;
196 elseif (isempty (vodeoptions.OutputFcn)), vhaveoutputfunction = false;
197 else vhaveoutputfunction = true;
200 %# Implementation of the option OutputSel has been finished. This
201 %# option can be set by the user to another value than default value.
202 if (~isempty (vodeoptions.OutputSel)), vhaveoutputselection = true;
203 else vhaveoutputselection = false; end
205 %# Implementation of the option Refine has been finished. This option
206 %# can be set by the user to another value than default value.
207 if (isequal (vodeoptions.Refine, vodetemp.Refine)), vhaverefine = true;
208 else vhaverefine = false; end
210 %# Implementation of the option Stats has been finished. This option
211 %# can be set by the user to another value than default value.
213 %# Implementation of the option InitialStep has been finished. This
214 %# option can be set by the user to another value than default value.
215 if (isempty (vodeoptions.InitialStep) && ~vstepsizegiven)
216 vodeoptions.InitialStep = abs (vslot(1,1) - vslot(1,2)) / 10;
217 vodeoptions.InitialStep = vodeoptions.InitialStep / 10^vodeoptions.Refine;
218 warning ('OdePkg:InvalidOption', ...
219 'Option "InitialStep" not set, new value %f is used', vodeoptions.InitialStep);
222 %# Implementation of the option MaxStep has been finished. This option
223 %# can be set by the user to another value than default value.
224 if (isempty (vodeoptions.MaxStep) && ~vstepsizegiven)
225 vodeoptions.MaxStep = abs (vslot(1,1) - vslot(1,length (vslot))) / 10;
226 %# vodeoptions.MaxStep = vodeoptions.MaxStep / 10^vodeoptions.Refine;
227 warning ('OdePkg:InvalidOption', ...
228 'Option "MaxStep" not set, new value %f is used', vodeoptions.MaxStep);
231 %# Implementation of the option Events has been finished. This option
232 %# can be set by the user to another value than default value.
233 if (~isempty (vodeoptions.Events)), vhaveeventfunction = true;
234 else vhaveeventfunction = false; end
236 %# The options 'Jacobian', 'JPattern' and 'Vectorized' will be ignored
237 %# by this solver because this solver uses an explicit Runge-Kutta
238 %# method and therefore no Jacobian calculation is necessary
239 if (~isequal (vodeoptions.Jacobian, vodetemp.Jacobian))
240 warning ('OdePkg:InvalidOption', ...
241 'Option "Jacobian" will be ignored by this solver');
243 if (~isequal (vodeoptions.JPattern, vodetemp.JPattern))
244 warning ('OdePkg:InvalidOption', ...
245 'Option "JPattern" will be ignored by this solver');
247 if (~isequal (vodeoptions.Vectorized, vodetemp.Vectorized))
248 warning ('OdePkg:InvalidOption', ...
249 'Option "Vectorized" will be ignored by this solver');
251 if (~isequal (vodeoptions.NewtonTol, vodetemp.NewtonTol))
252 warning ('OdePkg:InvalidArgument', ...
253 'Option "NewtonTol" will be ignored by this solver');
255 if (~isequal (vodeoptions.MaxNewtonIterations,...
256 vodetemp.MaxNewtonIterations))
257 warning ('OdePkg:InvalidArgument', ...
258 'Option "MaxNewtonIterations" will be ignored by this solver');
261 %# Implementation of the option Mass has been finished. This option
262 %# can be set by the user to another value than default value.
263 if (~isempty (vodeoptions.Mass) && isnumeric (vodeoptions.Mass))
264 vhavemasshandle = false; vmass = vodeoptions.Mass; %# constant mass
265 elseif (isa (vodeoptions.Mass, 'function_handle'))
266 vhavemasshandle = true; %# mass defined by a function handle
267 else %# no mass matrix - creating a diag-matrix of ones for mass
268 vhavemasshandle = false; %# vmass = diag (ones (length (vinit), 1), 0);
271 %# Implementation of the option MStateDependence has been finished.
272 %# This option can be set by the user to another value than default
274 if (strcmp (vodeoptions.MStateDependence, 'none'))
275 vmassdependence = false;
276 else vmassdependence = true;
279 %# Other options that are not used by this solver. Print a warning
280 %# message to tell the user that the option(s) is/are ignored.
281 if (~isequal (vodeoptions.MvPattern, vodetemp.MvPattern))
282 warning ('OdePkg:InvalidOption', ...
283 'Option "MvPattern" will be ignored by this solver');
285 if (~isequal (vodeoptions.MassSingular, vodetemp.MassSingular))
286 warning ('OdePkg:InvalidOption', ...
287 'Option "MassSingular" will be ignored by this solver');
289 if (~isequal (vodeoptions.InitialSlope, vodetemp.InitialSlope))
290 warning ('OdePkg:InvalidOption', ...
291 'Option "InitialSlope" will be ignored by this solver');
293 if (~isequal (vodeoptions.MaxOrder, vodetemp.MaxOrder))
294 warning ('OdePkg:InvalidOption', ...
295 'Option "MaxOrder" will be ignored by this solver');
297 if (~isequal (vodeoptions.BDF, vodetemp.BDF))
298 warning ('OdePkg:InvalidOption', ...
299 'Option "BDF" will be ignored by this solver');
302 %# Starting the initialisation of the core solver ode78d
303 vtimestamp = vslot(1,1); %# timestamp = start time
304 vtimelength = length (vslot); %# length needed if fixed steps
305 vtimestop = vslot(1,vtimelength); %# stop time = last value
308 vstepsize = vodeoptions.InitialStep;
309 vminstepsize = (vtimestop - vtimestamp) / (1/eps);
310 else %# If step size is given then use the fixed time steps
311 vstepsize = abs (vslot(1,1) - vslot(1,2));
312 vminstepsize = eps; %# vslot(1,2) - vslot(1,1) - eps;
315 vretvaltime = vtimestamp; %# first timestamp output
316 if (vhaveoutputselection) %# first solution output
317 vretvalresult = vinit(vodeoptions.OutputSel);
318 else vretvalresult = vinit;
321 %# Initialize the OutputFcn
322 if (vhaveoutputfunction)
323 feval (vodeoptions.OutputFcn, vslot', ...
324 vretvalresult', 'init', vfunarguments{:});
327 %# Initialize the History
328 if (isnumeric (vhist))
330 vhavehistnumeric = true;
331 else %# it must be a function handle
332 for vcnt = 1:length (vlags);
333 vhmat(:,vcnt) = feval (vhist, (vslot(1)-vlags(vcnt)), vfunarguments{:});
335 vhavehistnumeric = false;
338 %# Initialize DDE variables for history calculation
339 vsaveddetime = [vtimestamp - vlags, vtimestamp]';
340 vsaveddeinput = [vhmat, vinit']';
341 vsavedderesult = [vhmat, vinit']';
343 %# Initialize the EventFcn
344 if (vhaveeventfunction)
345 odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
346 {vretvalresult', vhmat}, 'init', vfunarguments{:});
349 vpow = 1/8; %# MC2001: see p.91 in Ascher & Petzold
350 va = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; %# The 7(8) coefficients
351 1/18, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; %# Coefficients proved, tt 20060827
352 1/48, 1/16, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
353 1/32, 0, 3/32, 0, 0, 0, 0, 0, 0, 0, 0, 0;
354 5/16, 0, -75/64, 75/64, 0, 0, 0, 0, 0, 0, 0, 0;
355 3/80, 0, 0, 3/16, 3/20, 0, 0, 0, 0, 0, 0, 0;
356 29443841/614563906, 0, 0, 77736538/692538347, -28693883/1125000000, ...
357 23124283/1800000000, 0, 0, 0, 0, 0, 0;
358 16016141/946692911, 0, 0, 61564180/158732637, 22789713/633445777, ...
359 545815736/2771057229, -180193667/1043307555, 0, 0, 0, 0, 0;
360 39632708/573591083, 0, 0, -433636366/683701615, -421739975/2616292301, ...
361 100302831/723423059, 790204164/839813087, 800635310/3783071287, 0, 0, 0, 0;
362 246121993/1340847787, 0, 0, -37695042795/15268766246, -309121744/1061227803, ...
363 -12992083/490766935, 6005943493/2108947869, 393006217/1396673457, ...
364 123872331/1001029789, 0, 0, 0;
365 -1028468189/846180014, 0, 0, 8478235783/508512852, 1311729495/1432422823, ...
366 -10304129995/1701304382, -48777925059/3047939560, 15336726248/1032824649, ...
367 -45442868181/3398467696, 3065993473/597172653, 0, 0;
368 185892177/718116043, 0, 0, -3185094517/667107341, -477755414/1098053517, ...
369 -703635378/230739211, 5731566787/1027545527, 5232866602/850066563, ...
370 -4093664535/808688257, 3962137247/1805957418, 65686358/487910083, 0;
371 403863854/491063109, 0, 0, -5068492393/434740067, -411421997/543043805, ...
372 652783627/914296604, 11173962825/925320556, -13158990841/6184727034, ...
373 3936647629/1978049680, -160528059/685178525, 248638103/1413531060, 0];
374 vb7 = [13451932/455176623; 0; 0; 0; 0; -808719846/976000145; ...
375 1757004468/5645159321; 656045339/265891186; -3867574721/1518517206; ...
376 465885868/322736535; 53011238/667516719; 2/45; 0];
377 vb8 = [14005451/335480064; 0; 0; 0; 0; -59238493/1068277825; 181606767/758867731; ...
378 561292985/797845732; -1041891430/1371343529; 760417239/1151165299; ...
379 118820643/751138087; -528747749/2220607170; 1/4];
382 %# The solver main loop - stop if the endpoint has been reached
383 vcntloop = 2; vcntcycles = 1; vu = vinit; vk = vu' * zeros(1,13);
384 vcntiter = 0; vunhandledtermination = true;
385 while ((vtimestamp < vtimestop && vstepsize >= vminstepsize))
387 %# Hit the endpoint of the time slot exactely
388 if ((vtimestamp + vstepsize) > vtimestop)
389 vstepsize = vtimestop - vtimestamp; end
391 %# Estimate the thirteen results when using this solver
393 vthetime = vtimestamp + vc(j,1) * vstepsize;
394 vtheinput = vu' + vstepsize * vk(:,1:j-1) * va(j,1:j-1)';
395 %# Claculate the history values (or get them from an external
396 %# function) that are needed for the next step of solving
397 if (vhavehistnumeric)
398 for vcnt = 1:length (vlags)
399 %# Direct implementation of a 'quadrature cubic Hermite interpolation'
400 %# found at the Faculty for Mathematics of the University of Stuttgart
401 %# http://mo.mathematik.uni-stuttgart.de/inhalt/aussage/aussage1269
402 vnumb = find (vthetime - vlags(vcnt) >= vsaveddetime);
403 velem = min (vnumb(end), length (vsaveddetime) - 1);
404 vstep = vsaveddetime(velem+1) - vsaveddetime(velem);
405 vdiff = (vthetime - vlags(vcnt) - vsaveddetime(velem)) / vstep;
407 %# Calculation of the coefficients for the interpolation algorithm
408 vua = (1 + 2 * vdiff) * vsubs^2;
409 vub = (3 - 2 * vdiff) * vdiff^2;
410 vva = vstep * vdiff * vsubs^2;
411 vvb = -vstep * vsubs * vdiff^2;
412 vhmat(:,vcnt) = vua * vsaveddeinput(velem,:)' + ...
413 vub * vsaveddeinput(velem+1,:)' + ...
414 vva * vsavedderesult(velem,:)' + ...
415 vvb * vsavedderesult(velem+1,:)';
417 else %# the history must be a function handle
418 for vcnt = 1:length (vlags)
419 vhmat(:,vcnt) = feval ...
420 (vhist, vthetime - vlags(vcnt), vfunarguments{:});
424 if (vhavemasshandle) %# Handle only the dynamic mass matrix,
425 if (vmassdependence) %# constant mass matrices have already
426 vmass = feval ... %# been set before (if any)
427 (vodeoptions.Mass, vthetime, vtheinput, vfunarguments{:});
428 else %# if (vmassdependence == false)
429 vmass = feval ... %# then we only have the time argument
430 (vodeoptions.Mass, vthetime, vfunarguments{:});
432 vk(:,j) = vmass \ feval ...
433 (vfun, vthetime, vtheinput, vhmat, vfunarguments{:});
436 (vfun, vthetime, vtheinput, vhmat, vfunarguments{:});
440 %# Compute the 7th and the 8th order estimation
441 y7 = vu' + vstepsize * (vk * vb7);
442 y8 = vu' + vstepsize * (vk * vb8);
443 if (vhavenonnegative)
444 vu(vodeoptions.NonNegative) = abs (vu(vodeoptions.NonNegative));
445 y7(vodeoptions.NonNegative) = abs (y7(vodeoptions.NonNegative));
446 y8(vodeoptions.NonNegative) = abs (y8(vodeoptions.NonNegative));
448 vSaveVUForRefine = vu;
450 %# Calculate the absolute local truncation error and the acceptable error
454 vtau = max (vodeoptions.RelTol * vu', vodeoptions.AbsTol);
456 vdelta = norm (y8 - y7, Inf);
457 vtau = max (vodeoptions.RelTol * max (norm (vu', Inf), 1.0), ...
460 else %# if (vstepsizegiven == true)
461 vdelta = 1; vtau = 2;
464 %# If the error is acceptable then update the vretval variables
465 if (all (vdelta <= vtau))
466 vtimestamp = vtimestamp + vstepsize;
467 vu = y8'; %# MC2001: the higher order estimation as "local extrapolation"
468 vretvaltime(vcntloop,:) = vtimestamp;
469 if (vhaveoutputselection)
470 vretvalresult(vcntloop,:) = vu(vodeoptions.OutputSel);
472 vretvalresult(vcntloop,:) = vu;
474 vcntloop = vcntloop + 1; vcntiter = 0;
476 %# Update DDE values for next history calculation
477 vsaveddetime(end+1) = vtimestamp;
478 vsaveddeinput(end+1,:) = vtheinput';
479 vsavedderesult(end+1,:) = vu;
481 %# Call plot only if a valid result has been found, therefore this
482 %# code fragment has moved here. Stop integration if plot function
484 if (vhaveoutputfunction)
485 if (vhaverefine) %# Do interpolation
486 for vcnt = 0:vodeoptions.Refine %# Approximation between told and t
487 vapproxtime = (vcnt + 1) * vstepsize / (vodeoptions.Refine + 2);
488 vapproxvals = vSaveVUForRefine' + vapproxtime * (vk * vb8);
489 if (vhaveoutputselection)
490 vapproxvals = vapproxvals(vodeoptions.OutputSel);
492 feval (vodeoptions.OutputFcn, (vtimestamp - vstepsize) + vapproxtime, ...
493 vapproxvals, [], vfunarguments{:});
496 vpltret = feval (vodeoptions.OutputFcn, vtimestamp, ...
497 vretvalresult(vcntloop-1,:)', [], vfunarguments{:});
498 if (vpltret), vunhandledtermination = false; break; end
501 %# Call event only if a valid result has been found, therefore this
502 %# code fragment has moved here. Stop integration if veventbreak is
504 if (vhaveeventfunction)
506 odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
507 {vu(:), vhmat}, [], vfunarguments{:});
508 if (~isempty (vevent{1}) && vevent{1} == 1)
509 vretvaltime(vcntloop-1,:) = vevent{3}(end,:);
510 vretvalresult(vcntloop-1,:) = vevent{4}(end,:);
511 vunhandledtermination = false; break;
514 end %# If the error is acceptable ...
516 %# Update the step size for the next integration step
518 %# vdelta may be 0 or even negative - could be an iteration problem
519 vdelta = max (vdelta, eps);
520 vstepsize = min (vodeoptions.MaxStep, ...
521 min (0.8 * vstepsize * (vtau ./ vdelta) .^ vpow));
522 elseif (vstepsizegiven)
523 if (vcntloop < vtimelength)
524 vstepsize = vslot(1,vcntloop-1) - vslot(1,vcntloop-2);
528 %# Update counters that count the number of iteration cycles
529 vcntcycles = vcntcycles + 1; %# Needed for postprocessing
530 vcntiter = vcntiter + 1; %# Needed to find iteration problems
532 %# Stop solving because the last 1000 steps no successful valid
533 %# value has been found
534 if (vcntiter >= 5000)
535 error (['Solving has not been successful. The iterative', ...
536 ' integration loop exited at time t = %f before endpoint at', ...
537 ' tend = %f was reached. This happened because the iterative', ...
538 ' integration loop does not find a valid solution at this time', ...
539 ' stamp. Try to reduce the value of "InitialStep" and/or', ...
540 ' "MaxStep" with the command "odeset".\n'], vtimestamp, vtimestop);
545 %# Check if integration of the ode has been successful
546 if (vtimestamp < vtimestop)
547 if (vunhandledtermination == true)
548 error (['Solving has not been successful. The iterative', ...
549 ' integration loop exited at time t = %f', ...
550 ' before endpoint at tend = %f was reached. This may', ...
551 ' happen if the stepsize grows smaller than defined in', ...
552 ' vminstepsize. Try to reduce the value of "InitialStep" and/or', ...
553 ' "MaxStep" with the command "odeset".\n'], vtimestamp, vtimestop);
555 warning ('OdePkg:HideWarning', ...
556 ['Solver has been stopped by a call of "break" in', ...
557 ' the main iteration loop at time t = %f before endpoint at', ...
558 ' tend = %f was reached. This may happen because the @odeplot', ...
559 ' function returned "true" or the @event function returned "true".'], ...
560 vtimestamp, vtimestop);
564 %# Postprocessing, do whatever when terminating integration algorithm
565 if (vhaveoutputfunction) %# Cleanup plotter
566 feval (vodeoptions.OutputFcn, vtimestamp, ...
567 vretvalresult(vcntloop-1,:)', 'done', vfunarguments{:});
569 if (vhaveeventfunction) %# Cleanup event function handling
570 odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
571 {vretvalresult(vcntloop-1,:), vhmat}, 'done', vfunarguments{:});
574 %# Print additional information if option Stats is set
575 if (strcmp (vodeoptions.Stats, 'on'))
577 vnsteps = vcntloop-2; %# vcntloop from 2..end
578 vnfailed = (vcntcycles-1)-(vcntloop-2)+1; %# vcntcycl from 1..end
579 vnfevals = 13*(vcntcycles-1); %# number of ode evaluations
580 vndecomps = 0; %# number of LU decompositions
581 vnpds = 0; %# number of partial derivatives
582 vnlinsols = 0; %# no. of solutions of linear systems
583 %# Print cost statistics if no output argument is given
585 vmsg = fprintf (1, 'Number of successful steps: %d', vnsteps);
586 vmsg = fprintf (1, 'Number of failed attempts: %d', vnfailed);
587 vmsg = fprintf (1, 'Number of function calls: %d', vnfevals);
589 else vhavestats = false;
592 if (nargout == 1) %# Sort output variables, depends on nargout
593 varargout{1}.x = vretvaltime; %# Time stamps are saved in field x
594 varargout{1}.y = vretvalresult; %# Results are saved in field y
595 varargout{1}.solver = 'ode78d'; %# Solver name is saved in field solver
596 if (vhaveeventfunction)
597 varargout{1}.ie = vevent{2}; %# Index info which event occured
598 varargout{1}.xe = vevent{3}; %# Time info when an event occured
599 varargout{1}.ye = vevent{4}; %# Results when an event occured
602 varargout{1}.stats = struct;
603 varargout{1}.stats.nsteps = vnsteps;
604 varargout{1}.stats.nfailed = vnfailed;
605 varargout{1}.stats.nfevals = vnfevals;
606 varargout{1}.stats.npds = vnpds;
607 varargout{1}.stats.ndecomps = vndecomps;
608 varargout{1}.stats.nlinsols = vnlinsols;
610 elseif (nargout == 2)
611 varargout{1} = vretvaltime; %# Time stamps are first output argument
612 varargout{2} = vretvalresult; %# Results are second output argument
613 elseif (nargout == 5)
614 varargout{1} = vretvaltime; %# Same as (nargout == 2)
615 varargout{2} = vretvalresult; %# Same as (nargout == 2)
616 varargout{3} = []; %# LabMat doesn't accept lines like
617 varargout{4} = []; %# varargout{3} = varargout{4} = [];
619 if (vhaveeventfunction)
620 varargout{3} = vevent{3}; %# Time info when an event occured
621 varargout{4} = vevent{4}; %# Results when an event occured
622 varargout{5} = vevent{2}; %# Index info which event occured
624 %# else nothing will be returned, varargout{1} undefined
627 %! # We are using a "pseudo-DDE" implementation for all tests that
628 %! # are done for this function. We also define an Events and a
629 %! # pseudo-Mass implementation. For further tests we also define a
630 %! # reference solution (computed at high accuracy) and an OutputFcn.
631 %!function [vyd] = fexp (vt, vy, vz, varargin)
632 %! vyd(1,1) = exp (- vt) - vz(1); %# The DDEs that are
633 %! vyd(2,1) = vy(1) - vz(2); %# used for all examples
634 %!function [vval, vtrm, vdir] = feve (vt, vy, vz, varargin)
635 %! vval = fexp (vt, vy, vz); %# We use the derivatives
636 %! vtrm = zeros (2,1); %# don't stop solving here
637 %! vdir = ones (2,1); %# in positive direction
638 %!function [vval, vtrm, vdir] = fevn (vt, vy, vz, varargin)
639 %! vval = fexp (vt, vy, vz); %# We use the derivatives
640 %! vtrm = ones (2,1); %# stop solving here
641 %! vdir = ones (2,1); %# in positive direction
642 %!function [vmas] = fmas (vt, vy, vz, varargin)
643 %! vmas = [1, 0; 0, 1]; %# Dummy mass matrix for tests
644 %!function [vmas] = fmsa (vt, vy, vz, varargin)
645 %! vmas = sparse ([1, 0; 0, 1]); %# A dummy sparse matrix
646 %!function [vref] = fref () %# The reference solution
647 %! vref = [0.12194462133618, 0.01652432423938];
648 %!function [vout] = fout (vt, vy, vflag, varargin)
649 %! if (regexp (char (vflag), 'init') == 1)
650 %! if (any (size (vt) ~= [2, 1])) error ('"fout" step "init"'); end
651 %! elseif (isempty (vflag))
652 %! if (any (size (vt) ~= [1, 1])) error ('"fout" step "calc"'); end
654 %! elseif (regexp (char (vflag), 'done') == 1)
655 %! if (any (size (vt) ~= [1, 1])) error ('"fout" step "done"'); end
656 %! else error ('"fout" invalid vflag');
659 %! %# Turn off output of warning messages for all tests, turn them on
660 %! %# again if the last test is called
661 %!error %# input argument number one
662 %! warning ('off', 'OdePkg:InvalidOption');
663 %! B = ode78d (1, [0 5], [1; 0], 1, [1; 0]);
664 %!error %# input argument number two
665 %! B = ode78d (@fexp, 1, [1; 0], 1, [1; 0]);
666 %!error %# input argument number three
667 %! B = ode78d (@fexp, [0 5], 1, 1, [1; 0]);
668 %!error %# input argument number four
669 %! B = ode78d (@fexp, [0 5], [1; 0], [1; 1], [1; 0]);
670 %!error %# input argument number five
671 %! B = ode78d (@fexp, [0 5], [1; 0], 1, 1);
672 %!test %# one output argument
673 %! vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0]);
674 %! assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 0.2);
675 %! assert (isfield (vsol, 'solver'));
676 %! assert (vsol.solver, 'ode78d');
677 %!test %# two output arguments
678 %! [vt, vy] = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0]);
679 %! assert ([vt(end), vy(end,:)], [5, fref], 0.2);
680 %!test %# five output arguments and no Events
681 %! [vt, vy, vxe, vye, vie] = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0]);
682 %! assert ([vt(end), vy(end,:)], [5, fref], 0.2);
683 %! assert ([vie, vxe, vye], []);
684 %!test %# anonymous function instead of real function
685 %! faym = @(vt, vy, vz) [exp(-vt) - vz(1); vy(1) - vz(2)];
686 %! vsol = ode78d (faym, [0 5], [1; 0], 1, [1; 0]);
687 %! assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 0.2);
688 %!test %# extra input arguments passed trhough
689 %! vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], 'KL');
690 %! assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 0.2);
691 %!test %# empty OdePkg structure *but* extra input arguments
693 %! vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt, 12, 13, 'KL');
694 %! assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 0.2);
695 %!error %# strange OdePkg structure
696 %! vopt = struct ('foo', 1);
697 %! vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
698 %!test %# AbsTol option
699 %! vopt = odeset ('AbsTol', 1e-5);
700 %! vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
701 %! assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 0.2);
702 %!test %# AbsTol and RelTol option
703 %! vopt = odeset ('AbsTol', 1e-7, 'RelTol', 1e-7);
704 %! vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
705 %! assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 0.2);
706 %!test %# RelTol and NormControl option
707 %! vopt = odeset ('AbsTol', 1e-7, 'NormControl', 'on');
708 %! vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
709 %! assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 0.2);
710 %!test %# NonNegative for second component
711 %! vopt = odeset ('NonNegative', 1);
712 %! vsol = ode78d (@fexp, [0 2.5], [1; 0], 1, [1; 0], vopt);
713 %! assert ([vsol.x(end), vsol.y(end,:)], [2.5, 0.001, 0.237], 0.2);
714 %!test %# Details of OutputSel and Refine can't be tested
715 %! vopt = odeset ('OutputFcn', @fout, 'OutputSel', 1, 'Refine', 5);
716 %! vsol = ode78d (@fexp, [0 2.5], [1; 0], 1, [1; 0], vopt);
717 %!test %# Stats must add further elements in vsol
718 %! vopt = odeset ('Stats', 'on');
719 %! vsol = ode78d (@fexp, [0 2.5], [1; 0], 1, [1; 0], vopt);
720 %! assert (isfield (vsol, 'stats'));
721 %! assert (isfield (vsol.stats, 'nsteps'));
722 %!test %# InitialStep option
723 %! vopt = odeset ('InitialStep', 1e-8);
724 %! vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
725 %! assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 0.2);
726 %!test %# MaxStep option
727 %! vopt = odeset ('MaxStep', 1e-2);
728 %! vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
729 %! assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 0.2);
730 %!test %# Events option add further elements in vsol
731 %! vopt = odeset ('Events', @feve);
732 %! vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
733 %! assert (isfield (vsol, 'ie'));
734 %! assert (vsol.ie, [1; 1]);
735 %! assert (isfield (vsol, 'xe'));
736 %! assert (isfield (vsol, 'ye'));
737 %!test %# Events option, now stop integration
738 %! warning ('off', 'OdePkg:HideWarning');
739 %! vopt = odeset ('Events', @fevn, 'NormControl', 'on');
740 %! vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
741 %! assert ([vsol.ie, vsol.xe, vsol.ye], ...
742 %! [1.0000, 2.9219, -0.2127, -0.2671], 0.2);
743 %!test %# Events option, five output arguments
744 %! vopt = odeset ('Events', @fevn, 'NormControl', 'on');
745 %! [vt, vy, vxe, vye, vie] = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
746 %! assert ([vie, vxe, vye], ...
747 %! [1.0000, 2.9219, -0.2127, -0.2671], 0.2);
749 %! %# test for Jacobian option is missing
750 %! %# test for Jacobian (being a sparse matrix) is missing
751 %! %# test for JPattern option is missing
752 %! %# test for Vectorized option is missing
753 %! %# test for NewtonTol option is missing
754 %! %# test for MaxNewtonIterations option is missing
756 %!test %# Mass option as function
757 %! vopt = odeset ('Mass', eye (2,2));
758 %! vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
759 %! assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 0.2);
760 %!test %# Mass option as matrix
761 %! vopt = odeset ('Mass', eye (2,2));
762 %! vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
763 %! assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 0.2);
764 %!test %# Mass option as sparse matrix
765 %! vopt = odeset ('Mass', sparse (eye (2,2)));
766 %! vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
767 %! assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 0.2);
768 %!test %# Mass option as function and sparse matrix
769 %! vopt = odeset ('Mass', @fmsa);
770 %! vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
771 %! assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 0.2);
772 %!test %# Mass option as function and MStateDependence
773 %! vopt = odeset ('Mass', @fmas, 'MStateDependence', 'strong');
774 %! vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
775 %! assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 0.2);
776 %!test %# Set BDF option to something else than default
777 %! vopt = odeset ('BDF', 'on');
778 %! [vt, vy] = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
779 %! assert ([vt(end), vy(end,:)], [5, fref], 0.5);
781 %! %# test for MvPattern option is missing
782 %! %# test for InitialSlope option is missing
783 %! %# test for MaxOrder option is missing
785 %! warning ('on', 'OdePkg:InvalidOption');
787 %# Local Variables: ***