1 %# Copyright (C) 2006-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{}] =} ode45 (@var{@@fun}, @var{slot}, @var{init}, [@var{opt}], [@var{par1}, @var{par2}, @dots{}])
19 %# @deftypefnx {Command} {[@var{sol}] =} ode45 (@var{@@fun}, @var{slot}, @var{init}, [@var{opt}], [@var{par1}, @var{par2}, @dots{}])
20 %# @deftypefnx {Command} {[@var{t}, @var{y}, [@var{xe}, @var{ye}, @var{ie}]] =} ode45 (@var{@@fun}, @var{slot}, @var{init}, [@var{opt}], [@var{par1}, @var{par2}, @dots{}])
22 %# This function file can be used to solve a set of non--stiff ordinary differential equations (non--stiff ODEs) or non--stiff differential algebraic equations (non--stiff DAEs) with the well known explicit Runge--Kutta method of order (4,5).
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 ODEs 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{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 %# If this function is called with one return argument then return the solution @var{sol} of type structure array after solving the set of ODEs. 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}.
28 %# 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.
30 %# For example, solve an anonymous implementation of the Van der Pol equation
33 %# fvdb = @@(vt,vy) [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)];
35 %# vopt = odeset ("RelTol", 1e-3, "AbsTol", 1e-3, \
36 %# "NormControl", "on", "OutputFcn", @@odeplot);
37 %# ode45 (fvdb, [0 20], [2 0], vopt);
44 %# 20010703 the function file "ode45.m" was written by Marc Compere
45 %# under the GPL for the use with this software. This function has been
46 %# taken as a base for the following implementation.
47 %# 20060810, Thomas Treichl
48 %# This function was adapted to the new syntax that is used by the
49 %# new OdePkg for Octave and is compatible to Matlab's ode45.
51 function [varargout] = ode45 (vfun, vslot, vinit, varargin)
53 if (nargin == 0) %# Check number and types of all input arguments
55 error ('OdePkg:InvalidArgument', ...
56 'Number of input arguments must be greater than zero');
61 elseif ~(isa (vfun, 'function_handle') || isa (vfun, 'inline'))
62 error ('OdePkg:InvalidArgument', ...
63 'First input argument must be a valid function handle');
65 elseif (~isvector (vslot) || length (vslot) < 2)
66 error ('OdePkg:InvalidArgument', ...
67 'Second input argument must be a valid vector');
69 elseif (~isvector (vinit) || ~isnumeric (vinit))
70 error ('OdePkg:InvalidArgument', ...
71 'Third input argument must be a valid numerical value');
75 if (~isstruct (varargin{1}))
76 %# varargin{1:len} are parameters for vfun
78 vfunarguments = varargin;
80 elseif (length (varargin) > 1)
81 %# varargin{1} is an OdePkg options structure vopt
82 vodeoptions = odepkg_structure_check (varargin{1}, 'ode45');
83 vfunarguments = {varargin{2:length(varargin)}};
85 else %# if (isstruct (varargin{1}))
86 vodeoptions = odepkg_structure_check (varargin{1}, 'ode45');
91 else %# if (nargin == 3)
96 %# Start preprocessing, have a look which options are set in
97 %# vodeoptions, check if an invalid or unused option is set
98 vslot = vslot(:).'; %# Create a row vector
99 vinit = vinit(:).'; %# Create a row vector
100 if (length (vslot) > 2) %# Step size checking
101 vstepsizefixed = true;
103 vstepsizefixed = false;
106 %# Get the default options that can be set with 'odeset' temporarily
109 %# Implementation of the option RelTol has been finished. This option
110 %# can be set by the user to another value than default value.
111 if (isempty (vodeoptions.RelTol) && ~vstepsizefixed)
112 vodeoptions.RelTol = 1e-6;
113 warning ('OdePkg:InvalidArgument', ...
114 'Option "RelTol" not set, new value %f is used', vodeoptions.RelTol);
115 elseif (~isempty (vodeoptions.RelTol) && vstepsizefixed)
116 warning ('OdePkg:InvalidArgument', ...
117 'Option "RelTol" will be ignored if fixed time stamps are given');
120 %# Implementation of the option AbsTol has been finished. This option
121 %# can be set by the user to another value than default value.
122 if (isempty (vodeoptions.AbsTol) && ~vstepsizefixed)
123 vodeoptions.AbsTol = 1e-6;
124 warning ('OdePkg:InvalidArgument', ...
125 'Option "AbsTol" not set, new value %f is used', vodeoptions.AbsTol);
126 elseif (~isempty (vodeoptions.AbsTol) && vstepsizefixed)
127 warning ('OdePkg:InvalidArgument', ...
128 'Option "AbsTol" will be ignored if fixed time stamps are given');
130 vodeoptions.AbsTol = vodeoptions.AbsTol(:); %# Create column vector
133 %# Implementation of the option NormControl has been finished. This
134 %# option can be set by the user to another value than default value.
135 if (strcmp (vodeoptions.NormControl, 'on')) vnormcontrol = true;
136 else vnormcontrol = false; end
138 %# Implementation of the option NonNegative has been finished. This
139 %# option can be set by the user to another value than default value.
140 if (~isempty (vodeoptions.NonNegative))
141 if (isempty (vodeoptions.Mass)), vhavenonnegative = true;
143 vhavenonnegative = false;
144 warning ('OdePkg:InvalidArgument', ...
145 'Option "NonNegative" will be ignored if mass matrix is set');
147 else vhavenonnegative = false;
150 %# Implementation of the option OutputFcn has been finished. This
151 %# option can be set by the user to another value than default value.
152 if (isempty (vodeoptions.OutputFcn) && nargout == 0)
153 vodeoptions.OutputFcn = @odeplot;
154 vhaveoutputfunction = true;
155 elseif (isempty (vodeoptions.OutputFcn)), vhaveoutputfunction = false;
156 else vhaveoutputfunction = true;
159 %# Implementation of the option OutputSel has been finished. This
160 %# option can be set by the user to another value than default value.
161 if (~isempty (vodeoptions.OutputSel)), vhaveoutputselection = true;
162 else vhaveoutputselection = false; end
164 %# Implementation of the option OutputSave has been finished. This
165 %# option can be set by the user to another value than default value.
166 if (isempty (vodeoptions.OutputSave)), vodeoptions.OutputSave = 1;
169 %# Implementation of the option Refine has been finished. This option
170 %# can be set by the user to another value than default value.
171 if (vodeoptions.Refine > 0), vhaverefine = true;
172 else vhaverefine = false; end
174 %# Implementation of the option Stats has been finished. This option
175 %# can be set by the user to another value than default value.
177 %# Implementation of the option InitialStep has been finished. This
178 %# option can be set by the user to another value than default value.
179 if (isempty (vodeoptions.InitialStep) && ~vstepsizefixed)
180 vodeoptions.InitialStep = (vslot(1,2) - vslot(1,1)) / 10;
181 vodeoptions.InitialStep = vodeoptions.InitialStep / 10^vodeoptions.Refine;
182 warning ('OdePkg:InvalidArgument', ...
183 'Option "InitialStep" not set, new value %f is used', vodeoptions.InitialStep);
186 %# Implementation of the option MaxStep has been finished. This option
187 %# can be set by the user to another value than default value.
188 if (isempty (vodeoptions.MaxStep) && ~vstepsizefixed)
189 vodeoptions.MaxStep = abs (vslot(1,2) - vslot(1,1)) / 10;
190 warning ('OdePkg:InvalidArgument', ...
191 'Option "MaxStep" not set, new value %f is used', vodeoptions.MaxStep);
194 %# Implementation of the option Events has been finished. This option
195 %# can be set by the user to another value than default value.
196 if (~isempty (vodeoptions.Events)), vhaveeventfunction = true;
197 else vhaveeventfunction = false; end
199 %# The options 'Jacobian', 'JPattern' and 'Vectorized' will be ignored
200 %# by this solver because this solver uses an explicit Runge-Kutta
201 %# method and therefore no Jacobian calculation is necessary
202 if (~isequal (vodeoptions.Jacobian, vodetemp.Jacobian))
203 warning ('OdePkg:InvalidArgument', ...
204 'Option "Jacobian" will be ignored by this solver');
206 if (~isequal (vodeoptions.JPattern, vodetemp.JPattern))
207 warning ('OdePkg:InvalidArgument', ...
208 'Option "JPattern" will be ignored by this solver');
210 if (~isequal (vodeoptions.Vectorized, vodetemp.Vectorized))
211 warning ('OdePkg:InvalidArgument', ...
212 'Option "Vectorized" will be ignored by this solver');
214 if (~isequal (vodeoptions.NewtonTol, vodetemp.NewtonTol))
215 warning ('OdePkg:InvalidArgument', ...
216 'Option "NewtonTol" will be ignored by this solver');
218 if (~isequal (vodeoptions.MaxNewtonIterations,...
219 vodetemp.MaxNewtonIterations))
220 warning ('OdePkg:InvalidArgument', ...
221 'Option "MaxNewtonIterations" will be ignored by this solver');
224 %# Implementation of the option Mass has been finished. This option
225 %# can be set by the user to another value than default value.
226 if (~isempty (vodeoptions.Mass) && isnumeric (vodeoptions.Mass))
227 vhavemasshandle = false; vmass = vodeoptions.Mass; %# constant mass
228 elseif (isa (vodeoptions.Mass, 'function_handle'))
229 vhavemasshandle = true; %# mass defined by a function handle
230 else %# no mass matrix - creating a diag-matrix of ones for mass
231 vhavemasshandle = false; %# vmass = diag (ones (length (vinit), 1), 0);
234 %# Implementation of the option MStateDependence has been finished.
235 %# This option can be set by the user to another value than default
237 if (strcmp (vodeoptions.MStateDependence, 'none'))
238 vmassdependence = false;
239 else vmassdependence = true;
242 %# Other options that are not used by this solver. Print a warning
243 %# message to tell the user that the option(s) is/are ignored.
244 if (~isequal (vodeoptions.MvPattern, vodetemp.MvPattern))
245 warning ('OdePkg:InvalidArgument', ...
246 'Option "MvPattern" will be ignored by this solver');
248 if (~isequal (vodeoptions.MassSingular, vodetemp.MassSingular))
249 warning ('OdePkg:InvalidArgument', ...
250 'Option "MassSingular" will be ignored by this solver');
252 if (~isequal (vodeoptions.InitialSlope, vodetemp.InitialSlope))
253 warning ('OdePkg:InvalidArgument', ...
254 'Option "InitialSlope" will be ignored by this solver');
256 if (~isequal (vodeoptions.MaxOrder, vodetemp.MaxOrder))
257 warning ('OdePkg:InvalidArgument', ...
258 'Option "MaxOrder" will be ignored by this solver');
260 if (~isequal (vodeoptions.BDF, vodetemp.BDF))
261 warning ('OdePkg:InvalidArgument', ...
262 'Option "BDF" will be ignored by this solver');
265 %# Starting the initialisation of the core solver ode45
266 vtimestamp = vslot(1,1); %# timestamp = start time
267 vtimelength = length (vslot); %# length needed if fixed steps
268 vtimestop = vslot(1,vtimelength); %# stop time = last value
269 %# 20110611, reported by Nils Strunk
270 %# Make it possible to solve equations from negativ to zero,
271 %# eg. vres = ode45 (@(t,y) y, [-2 0], 2);
272 vdirection = sign (vtimestop - vtimestamp); %# Direction flag
275 if (sign (vodeoptions.InitialStep) == vdirection)
276 vstepsize = vodeoptions.InitialStep;
277 else %# Fix wrong direction of InitialStep.
278 vstepsize = - vodeoptions.InitialStep;
280 vminstepsize = (vtimestop - vtimestamp) / (1/eps);
281 else %# If step size is given then use the fixed time steps
282 vstepsize = vslot(1,2) - vslot(1,1);
283 vminstepsize = sign (vstepsize) * eps;
286 vretvaltime = vtimestamp; %# first timestamp output
287 vretvalresult = vinit; %# first solution output
289 %# Initialize the OutputFcn
290 if (vhaveoutputfunction)
291 if (vhaveoutputselection) vretout = vretvalresult(vodeoptions.OutputSel);
292 else vretout = vretvalresult; end
293 feval (vodeoptions.OutputFcn, vslot.', ...
294 vretout.', 'init', vfunarguments{:});
297 %# Initialize the EventFcn
298 if (vhaveeventfunction)
299 odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
300 vretvalresult.', 'init', vfunarguments{:});
303 vpow = 1/5; %# 20071016, reported by Luis Randez
304 va = [0, 0, 0, 0, 0; %# The Runge-Kutta-Fehlberg 4(5) coefficients
305 1/4, 0, 0, 0, 0; %# Coefficients proved on 20060827
306 3/32, 9/32, 0, 0, 0; %# See p.91 in Ascher & Petzold
307 1932/2197, -7200/2197, 7296/2197, 0, 0;
308 439/216, -8, 3680/513, -845/4104, 0;
309 -8/27, 2, -3544/2565, 1859/4104, -11/40];
310 %# 4th and 5th order b-coefficients
311 vb4 = [25/216; 0; 1408/2565; 2197/4104; -1/5; 0];
312 vb5 = [16/135; 0; 6656/12825; 28561/56430; -9/50; 2/55];
315 %# The solver main loop - stop if the endpoint has been reached
316 vcntloop = 2; vcntcycles = 1; vu = vinit; vk = vu.' * zeros(1,6);
317 vcntiter = 0; vunhandledtermination = true; vcntsave = 2;
318 while ((vdirection * (vtimestamp) < vdirection * (vtimestop)) && ...
319 (vdirection * (vstepsize) >= vdirection * (vminstepsize)))
321 %# Hit the endpoint of the time slot exactely
322 if (vdirection * (vtimestamp + vstepsize) > vdirection * vtimestop)
323 %# vstepsize = vtimestop - vdirection * vtimestamp;
324 %# 20110611, reported by Nils Strunk
325 %# The endpoint of the time slot must be hit exactly,
326 %# eg. vsol = ode45 (@(t,y) y, [0 -1], 1);
327 vstepsize = vdirection * abs (abs (vtimestop) - abs (vtimestamp));
330 %# Estimate the six results when using this solver
332 vthetime = vtimestamp + vc(j,1) * vstepsize;
333 vtheinput = vu.' + vstepsize * vk(:,1:j-1) * va(j,1:j-1).';
334 if (vhavemasshandle) %# Handle only the dynamic mass matrix,
335 if (vmassdependence) %# constant mass matrices have already
336 vmass = feval ... %# been set before (if any)
337 (vodeoptions.Mass, vthetime, vtheinput, vfunarguments{:});
338 else %# if (vmassdependence == false)
339 vmass = feval ... %# then we only have the time argument
340 (vodeoptions.Mass, vthetime, vfunarguments{:});
342 vk(:,j) = vmass \ feval ...
343 (vfun, vthetime, vtheinput, vfunarguments{:});
346 (vfun, vthetime, vtheinput, vfunarguments{:});
350 %# Compute the 4th and the 5th order estimation
351 y4 = vu.' + vstepsize * (vk * vb4);
352 y5 = vu.' + vstepsize * (vk * vb5);
353 if (vhavenonnegative)
354 vu(vodeoptions.NonNegative) = abs (vu(vodeoptions.NonNegative));
355 y4(vodeoptions.NonNegative) = abs (y4(vodeoptions.NonNegative));
356 y5(vodeoptions.NonNegative) = abs (y5(vodeoptions.NonNegative));
358 if (vhaveoutputfunction && vhaverefine)
359 vSaveVUForRefine = vu;
362 %# Calculate the absolute local truncation error and the acceptable error
365 vdelta = abs (y5 - y4);
366 vtau = max (vodeoptions.RelTol * abs (vu.'), vodeoptions.AbsTol);
368 vdelta = norm (y5 - y4, Inf);
369 vtau = max (vodeoptions.RelTol * max (norm (vu.', Inf), 1.0), ...
372 else %# if (vstepsizefixed == true)
373 vdelta = 1; vtau = 2;
376 %# If the error is acceptable then update the vretval variables
377 if (all (vdelta <= vtau))
378 vtimestamp = vtimestamp + vstepsize;
379 vu = y5.'; %# MC2001: the higher order estimation as "local extrapolation"
380 %# Save the solution every vodeoptions.OutputSave steps
381 if (mod (vcntloop-1,vodeoptions.OutputSave) == 0)
382 vretvaltime(vcntsave,:) = vtimestamp;
383 vretvalresult(vcntsave,:) = vu;
384 vcntsave = vcntsave + 1;
386 vcntloop = vcntloop + 1; vcntiter = 0;
388 %# Call plot only if a valid result has been found, therefore this
389 %# code fragment has moved here. Stop integration if plot function
391 if (vhaveoutputfunction)
392 for vcnt = 0:vodeoptions.Refine %# Approximation between told and t
393 if (vhaverefine) %# Do interpolation
394 vapproxtime = (vcnt + 1) * vstepsize / (vodeoptions.Refine + 2);
395 vapproxvals = vSaveVUForRefine.' + vapproxtime * (vk * vb5);
396 vapproxtime = (vtimestamp - vstepsize) + vapproxtime;
399 vapproxtime = vtimestamp;
401 if (vhaveoutputselection)
402 vapproxvals = vapproxvals(vodeoptions.OutputSel);
404 vpltret = feval (vodeoptions.OutputFcn, vapproxtime, ...
405 vapproxvals, [], vfunarguments{:});
406 if vpltret %# Leave refinement loop
410 if (vpltret) %# Leave main loop
411 vunhandledtermination = false;
416 %# Call event only if a valid result has been found, therefore this
417 %# code fragment has moved here. Stop integration if veventbreak is
419 if (vhaveeventfunction)
421 odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
422 vu(:), [], vfunarguments{:});
423 if (~isempty (vevent{1}) && vevent{1} == 1)
424 vretvaltime(vcntloop-1,:) = vevent{3}(end,:);
425 vretvalresult(vcntloop-1,:) = vevent{4}(end,:);
426 vunhandledtermination = false; break;
429 end %# If the error is acceptable ...
431 %# Update the step size for the next integration step
433 %# 20080425, reported by Marco Caliari
434 %# vdelta cannot be negative (because of the absolute value that
435 %# has been introduced) but it could be 0, then replace the zeros
436 %# with the maximum value of vdelta
437 vdelta(find (vdelta == 0)) = max (vdelta);
438 %# It could happen that max (vdelta) == 0 (ie. that the original
439 %# vdelta was 0), in that case we double the previous vstepsize
440 vdelta(find (vdelta == 0)) = max (vtau) .* (0.4 ^ (1 / vpow));
443 vstepsize = min (vodeoptions.MaxStep, ...
444 min (0.8 * vstepsize * (vtau ./ vdelta) .^ vpow));
446 vstepsize = max (- vodeoptions.MaxStep, ...
447 max (0.8 * vstepsize * (vtau ./ vdelta) .^ vpow));
450 else %# if (vstepsizefixed)
451 if (vcntloop <= vtimelength)
452 vstepsize = vslot(vcntloop) - vslot(vcntloop-1);
453 else %# Get out of the main integration loop
458 %# Update counters that count the number of iteration cycles
459 vcntcycles = vcntcycles + 1; %# Needed for cost statistics
460 vcntiter = vcntiter + 1; %# Needed to find iteration problems
462 %# Stop solving because the last 1000 steps no successful valid
463 %# value has been found
464 if (vcntiter >= 5000)
465 error (['Solving has not been successful. The iterative', ...
466 ' integration loop exited at time t = %f before endpoint at', ...
467 ' tend = %f was reached. This happened because the iterative', ...
468 ' integration loop does not find a valid solution at this time', ...
469 ' stamp. Try to reduce the value of "InitialStep" and/or', ...
470 ' "MaxStep" with the command "odeset".\n'], vtimestamp, vtimestop);
475 %# Check if integration of the ode has been successful
476 if (vdirection * vtimestamp < vdirection * vtimestop)
477 if (vunhandledtermination == true)
478 error ('OdePkg:InvalidArgument', ...
479 ['Solving has not been successful. The iterative', ...
480 ' integration loop exited at time t = %f', ...
481 ' before endpoint at tend = %f was reached. This may', ...
482 ' happen if the stepsize grows smaller than defined in', ...
483 ' vminstepsize. Try to reduce the value of "InitialStep" and/or', ...
484 ' "MaxStep" with the command "odeset".\n'], vtimestamp, vtimestop);
486 warning ('OdePkg:InvalidArgument', ...
487 ['Solver has been stopped by a call of "break" in', ...
488 ' the main iteration loop at time t = %f before endpoint at', ...
489 ' tend = %f was reached. This may happen because the @odeplot', ...
490 ' function returned "true" or the @event function returned "true".'], ...
491 vtimestamp, vtimestop);
495 %# Postprocessing, do whatever when terminating integration algorithm
496 if (vhaveoutputfunction) %# Cleanup plotter
497 feval (vodeoptions.OutputFcn, vtimestamp, ...
498 vu.', 'done', vfunarguments{:});
500 if (vhaveeventfunction) %# Cleanup event function handling
501 odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
502 vu.', 'done', vfunarguments{:});
504 %# Save the last step, if not already saved
505 if (mod (vcntloop-2,vodeoptions.OutputSave) ~= 0)
506 vretvaltime(vcntsave,:) = vtimestamp;
507 vretvalresult(vcntsave,:) = vu;
510 %# Print additional information if option Stats is set
511 if (strcmp (vodeoptions.Stats, 'on'))
513 vnsteps = vcntloop-2; %# vcntloop from 2..end
514 vnfailed = (vcntcycles-1)-(vcntloop-2)+1; %# vcntcycl from 1..end
515 vnfevals = 6*(vcntcycles-1); %# number of ode evaluations
516 vndecomps = 0; %# number of LU decompositions
517 vnpds = 0; %# number of partial derivatives
518 vnlinsols = 0; %# no. of solutions of linear systems
519 %# Print cost statistics if no output argument is given
521 vmsg = fprintf (1, 'Number of successful steps: %d\n', vnsteps);
522 vmsg = fprintf (1, 'Number of failed attempts: %d\n', vnfailed);
523 vmsg = fprintf (1, 'Number of function calls: %d\n', vnfevals);
529 if (nargout == 1) %# Sort output variables, depends on nargout
530 varargout{1}.x = vretvaltime; %# Time stamps are saved in field x
531 varargout{1}.y = vretvalresult; %# Results are saved in field y
532 varargout{1}.solver = 'ode45'; %# Solver name is saved in field solver
533 if (vhaveeventfunction)
534 varargout{1}.ie = vevent{2}; %# Index info which event occured
535 varargout{1}.xe = vevent{3}; %# Time info when an event occured
536 varargout{1}.ye = vevent{4}; %# Results when an event occured
539 varargout{1}.stats = struct;
540 varargout{1}.stats.nsteps = vnsteps;
541 varargout{1}.stats.nfailed = vnfailed;
542 varargout{1}.stats.nfevals = vnfevals;
543 varargout{1}.stats.npds = vnpds;
544 varargout{1}.stats.ndecomps = vndecomps;
545 varargout{1}.stats.nlinsols = vnlinsols;
547 elseif (nargout == 2)
548 varargout{1} = vretvaltime; %# Time stamps are first output argument
549 varargout{2} = vretvalresult; %# Results are second output argument
550 elseif (nargout == 5)
551 varargout{1} = vretvaltime; %# Same as (nargout == 2)
552 varargout{2} = vretvalresult; %# Same as (nargout == 2)
553 varargout{3} = []; %# LabMat doesn't accept lines like
554 varargout{4} = []; %# varargout{3} = varargout{4} = [];
556 if (vhaveeventfunction)
557 varargout{3} = vevent{3}; %# Time info when an event occured
558 varargout{4} = vevent{4}; %# Results when an event occured
559 varargout{5} = vevent{2}; %# Index info which event occured
564 %! # We are using the "Van der Pol" implementation for all tests that
565 %! # are done for this function. We also define a Jacobian, Events,
566 %! # pseudo-Mass implementation. For further tests we also define a
567 %! # reference solution (computed at high accuracy) and an OutputFcn
568 %!function [ydot] = fpol (vt, vy, varargin) %# The Van der Pol
569 %! ydot = [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)];
570 %!function [vjac] = fjac (vt, vy, varargin) %# its Jacobian
571 %! vjac = [0, 1; -1 - 2 * vy(1) * vy(2), 1 - vy(1)^2];
572 %!function [vjac] = fjcc (vt, vy, varargin) %# sparse type
573 %! vjac = sparse ([0, 1; -1 - 2 * vy(1) * vy(2), 1 - vy(1)^2]);
574 %!function [vval, vtrm, vdir] = feve (vt, vy, varargin)
575 %! vval = fpol (vt, vy, varargin); %# We use the derivatives
576 %! vtrm = zeros (2,1); %# that's why component 2
577 %! vdir = ones (2,1); %# seems to not be exact
578 %!function [vval, vtrm, vdir] = fevn (vt, vy, varargin)
579 %! vval = fpol (vt, vy, varargin); %# We use the derivatives
580 %! vtrm = ones (2,1); %# that's why component 2
581 %! vdir = ones (2,1); %# seems to not be exact
582 %!function [vmas] = fmas (vt, vy)
583 %! vmas = [1, 0; 0, 1]; %# Dummy mass matrix for tests
584 %!function [vmas] = fmsa (vt, vy)
585 %! vmas = sparse ([1, 0; 0, 1]); %# A sparse dummy matrix
586 %!function [vref] = fref () %# The computed reference sol
587 %! vref = [0.32331666704577, -1.83297456798624];
588 %!function [vout] = fout (vt, vy, vflag, varargin)
589 %! if (regexp (char (vflag), 'init') == 1)
590 %! if (any (size (vt) ~= [2, 1])) error ('"fout" step "init"'); end
591 %! elseif (isempty (vflag))
592 %! if (any (size (vt) ~= [1, 1])) error ('"fout" step "calc"'); end
594 %! elseif (regexp (char (vflag), 'done') == 1)
595 %! if (any (size (vt) ~= [1, 1])) error ('"fout" step "done"'); end
596 %! else error ('"fout" invalid vflag');
599 %! %# Turn off output of warning messages for all tests, turn them on
600 %! %# again if the last test is called
601 %!error %# input argument number one
602 %! warning ('off', 'OdePkg:InvalidArgument');
603 %! B = ode45 (1, [0 25], [3 15 1]);
604 %!error %# input argument number two
605 %! B = ode45 (@fpol, 1, [3 15 1]);
606 %!error %# input argument number three
607 %! B = ode45 (@flor, [0 25], 1);
608 %!test %# one output argument
609 %! vsol = ode45 (@fpol, [0 2], [2 0]);
610 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
611 %! assert (isfield (vsol, 'solver'));
612 %! assert (vsol.solver, 'ode45');
613 %!test %# two output arguments
614 %! [vt, vy] = ode45 (@fpol, [0 2], [2 0]);
615 %! assert ([vt(end), vy(end,:)], [2, fref], 1e-3);
616 %!test %# five output arguments and no Events
617 %! [vt, vy, vxe, vye, vie] = ode45 (@fpol, [0 2], [2 0]);
618 %! assert ([vt(end), vy(end,:)], [2, fref], 1e-3);
619 %! assert ([vie, vxe, vye], []);
620 %!test %# anonymous function instead of real function
621 %! fvdb = @(vt,vy) [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)];
622 %! vsol = ode45 (fvdb, [0 2], [2 0]);
623 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
624 %!test %# extra input arguments passed through
625 %! vsol = ode45 (@fpol, [0 2], [2 0], 12, 13, 'KL');
626 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
627 %!test %# empty OdePkg structure *but* extra input arguments
629 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt, 12, 13, 'KL');
630 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
631 %!error %# strange OdePkg structure
632 %! vopt = struct ('foo', 1);
633 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
634 %!test %# Solve vdp in fixed step sizes
635 %! vsol = ode45 (@fpol, [0:0.1:2], [2 0]);
636 %! assert (vsol.x(:), [0:0.1:2]');
637 %! assert (vsol.y(end,:), fref, 1e-3);
638 %!test %# Solve in backward direction starting at t=0
639 %! vref = [-1.205364552835178, 0.951542399860817];
640 %! vsol = ode45 (@fpol, [0 -2], [2 0]);
641 %! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-3);
642 %!test %# Solve in backward direction starting at t=2
643 %! vref = [-1.205364552835178, 0.951542399860817];
644 %! vsol = ode45 (@fpol, [2 -2], fref);
645 %! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-3);
646 %!test %# Solve another anonymous function in backward direction
647 %! vref = [-1, 0.367879437558975];
648 %! vsol = ode45 (@(t,y) y, [0 -1], 1);
649 %! assert ([vsol.x(end), vsol.y(end,:)], vref, 1e-3);
650 %!test %# Solve another anonymous function below zero
651 %! vref = [0, 14.77810590694212];
652 %! vsol = ode45 (@(t,y) y, [-2 0], 2);
653 %! assert ([vsol.x(end), vsol.y(end,:)], vref, 1e-3);
654 %!test %# AbsTol option
655 %! vopt = odeset ('AbsTol', 1e-5);
656 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
657 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
658 %!test %# AbsTol and RelTol option
659 %! vopt = odeset ('AbsTol', 1e-8, 'RelTol', 1e-8);
660 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
661 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
662 %!test %# RelTol and NormControl option -- higher accuracy
663 %! vopt = odeset ('RelTol', 1e-8, 'NormControl', 'on');
664 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
665 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-5);
666 %!test %# Keeps initial values while integrating
667 %! vopt = odeset ('NonNegative', 2);
668 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
669 %! assert ([vsol.x(end), vsol.y(end,:)], [2, 2, 0], 0.5);
670 %!test %# Details of OutputSel and Refine can't be tested
671 %! vopt = odeset ('OutputFcn', @fout, 'OutputSel', 1, 'Refine', 5);
672 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
673 %!test %# Details of OutputSave can't be tested
674 %! vopt = odeset ('OutputSave', 1, 'OutputSel', 1);
675 %! vsla = ode45 (@fpol, [0 2], [2 0], vopt);
676 %! vopt = odeset ('OutputSave', 2);
677 %! vslb = ode45 (@fpol, [0 2], [2 0], vopt);
678 %! assert (length (vsla.x) > length (vslb.x))
679 %!test %# Stats must add further elements in vsol
680 %! vopt = odeset ('Stats', 'on');
681 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
682 %! assert (isfield (vsol, 'stats'));
683 %! assert (isfield (vsol.stats, 'nsteps'));
684 %!test %# InitialStep option
685 %! vopt = odeset ('InitialStep', 1e-8);
686 %! vsol = ode45 (@fpol, [0 0.2], [2 0], vopt);
687 %! assert ([vsol.x(2)-vsol.x(1)], [1e-8], 1e-9);
688 %!test %# MaxStep option
689 %! vopt = odeset ('MaxStep', 1e-2);
690 %! vsol = ode45 (@fpol, [0 0.2], [2 0], vopt);
691 %! assert ([vsol.x(5)-vsol.x(4)], [1e-2], 1e-3);
692 %!test %# Events option add further elements in vsol
693 %! vopt = odeset ('Events', @feve);
694 %! vsol = ode45 (@fpol, [0 10], [2 0], vopt);
695 %! assert (isfield (vsol, 'ie'));
696 %! assert (vsol.ie(1), 2);
697 %! assert (isfield (vsol, 'xe'));
698 %! assert (isfield (vsol, 'ye'));
699 %!test %# Events option, now stop integration
700 %! vopt = odeset ('Events', @fevn, 'NormControl', 'on');
701 %! vsol = ode45 (@fpol, [0 10], [2 0], vopt);
702 %! assert ([vsol.ie, vsol.xe, vsol.ye], ...
703 %! [2.0, 2.496110, -0.830550, -2.677589], .5e-1);
704 %!test %# Events option, five output arguments
705 %! vopt = odeset ('Events', @fevn, 'NormControl', 'on');
706 %! [vt, vy, vxe, vye, vie] = ode45 (@fpol, [0 10], [2 0], vopt);
707 %! assert ([vie, vxe, vye], ...
708 %! [2.0, 2.496110, -0.830550, -2.677589], 1e-1);
709 %!test %# Jacobian option
710 %! vopt = odeset ('Jacobian', @fjac);
711 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
712 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
713 %!test %# Jacobian option and sparse return value
714 %! vopt = odeset ('Jacobian', @fjcc);
715 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
716 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
718 %! %# test for JPattern option is missing
719 %! %# test for Vectorized option is missing
720 %! %# test for NewtonTol option is missing
721 %! %# test for MaxNewtonIterations option is missing
723 %!test %# Mass option as function
724 %! vopt = odeset ('Mass', @fmas);
725 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
726 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
727 %!test %# Mass option as matrix
728 %! vopt = odeset ('Mass', eye (2,2));
729 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
730 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
731 %!test %# Mass option as sparse matrix
732 %! vopt = odeset ('Mass', sparse (eye (2,2)));
733 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
734 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
735 %!test %# Mass option as function and sparse matrix
736 %! vopt = odeset ('Mass', @fmsa);
737 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
738 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
739 %!test %# Mass option as function and MStateDependence
740 %! vopt = odeset ('Mass', @fmas, 'MStateDependence', 'strong');
741 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
742 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
743 %!test %# Set BDF option to something else than default
744 %! vopt = odeset ('BDF', 'on');
745 %! [vt, vy] = ode45 (@fpol, [0 2], [2 0], vopt);
746 %! assert ([vt(end), vy(end,:)], [2, fref], 1e-3);
748 %! %# test for MvPattern option is missing
749 %! %# test for InitialSlope option is missing
750 %! %# test for MaxOrder option is missing
752 %! warning ('on', 'OdePkg:InvalidArgument');
754 %# Local Variables: ***