1 %# Copyright (C) 2009-2012, Sebastian Schoeps <schoeps AT math DOT uni-wuppertal DOT de>
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{}] =} odebwe (@var{@@fun}, @var{slot}, @var{init}, [@var{opt}], [@var{par1}, @var{par2}, @dots{}])
19 %# @deftypefnx {Command} {[@var{sol}] =} odebwe (@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}]] =} odebwe (@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 stiff ordinary differential equations (stiff ODEs) or stiff differential algebraic equations (stiff DAEs) with the Backward Euler method.
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)];
34 %# vjac = @@(vt,vy) [0, 1; -1 - 2 * vy(1) * vy(2), 1 - vy(1)^2];
35 %# vopt = odeset ("RelTol", 1e-3, "AbsTol", 1e-3, \
36 %# "NormControl", "on", "OutputFcn", @@odeplot, \
38 %# odebwe (fvdb, [0 20], [2 0], vopt);
44 function [varargout] = odebwe (vfun, vslot, vinit, varargin)
46 if (nargin == 0) %# Check number and types of all input arguments
48 error ('OdePkg:InvalidArgument', ...
49 'Number of input arguments must be greater than zero');
54 elseif ~(isa (vfun, 'function_handle') || isa (vfun, 'inline'))
55 error ('OdePkg:InvalidArgument', ...
56 'First input argument must be a valid function handle');
58 elseif (~isvector (vslot) || length (vslot) < 2)
59 error ('OdePkg:InvalidArgument', ...
60 'Second input argument must be a valid vector');
62 elseif (~isvector (vinit) || ~isnumeric (vinit))
63 error ('OdePkg:InvalidArgument', ...
64 'Third input argument must be a valid numerical value');
68 if (~isstruct (varargin{1}))
69 %# varargin{1:len} are parameters for vfun
71 vfunarguments = varargin;
73 elseif (length (varargin) > 1)
74 %# varargin{1} is an OdePkg options structure vopt
75 vodeoptions = odepkg_structure_check (varargin{1}, 'odebwe');
76 vfunarguments = {varargin{2:length(varargin)}};
78 else %# if (isstruct (varargin{1}))
79 vodeoptions = odepkg_structure_check (varargin{1}, 'odebwe');
84 else %# if (nargin == 3)
89 %# Start preprocessing, have a look which options are set in
90 %# vodeoptions, check if an invalid or unused option is set
91 vslot = vslot(:).'; %# Create a row vector
92 vinit = vinit(:).'; %# Create a row vector
93 if (length (vslot) > 2) %# Step size checking
94 vstepsizefixed = true;
96 vstepsizefixed = false;
99 %# The adaptive method require a second estimate for
100 %# the comparsion, while the fixed step size algorithm
108 %# Get the default options that can be set with 'odeset' temporarily
111 %# Implementation of the option RelTol has been finished. This option
112 %# can be set by the user to another value than default value.
113 if (isempty (vodeoptions.RelTol) && ~vstepsizefixed)
114 vodeoptions.RelTol = 1e-6;
115 warning ('OdePkg:InvalidArgument', ...
116 'Option "RelTol" not set, new value %f is used', vodeoptions.RelTol);
117 elseif (~isempty (vodeoptions.RelTol) && vstepsizefixed)
118 warning ('OdePkg:InvalidArgument', ...
119 'Option "RelTol" will be ignored if fixed time stamps are given');
122 %# Implementation of the option AbsTol has been finished. This option
123 %# can be set by the user to another value than default value.
124 if (isempty (vodeoptions.AbsTol) && ~vstepsizefixed)
125 vodeoptions.AbsTol = 1e-6;
126 warning ('OdePkg:InvalidArgument', ...
127 'Option "AbsTol" not set, new value %f is used', vodeoptions.AbsTol);
128 elseif (~isempty (vodeoptions.AbsTol) && vstepsizefixed)
129 warning ('OdePkg:InvalidArgument', ...
130 'Option "AbsTol" will be ignored if fixed time stamps are given');
132 vodeoptions.AbsTol = vodeoptions.AbsTol(:); %# Create column vector
135 %# Implementation of the option NormControl has been finished. This
136 %# option can be set by the user to another value than default value.
137 if (strcmp (vodeoptions.NormControl, 'on')) vnormcontrol = true;
138 else vnormcontrol = false; end
140 %# Implementation of the option OutputFcn has been finished. This
141 %# option can be set by the user to another value than default value.
142 if (isempty (vodeoptions.OutputFcn) && nargout == 0)
143 vodeoptions.OutputFcn = @odeplot;
144 vhaveoutputfunction = true;
145 elseif (isempty (vodeoptions.OutputFcn)), vhaveoutputfunction = false;
146 else vhaveoutputfunction = true;
149 %# Implementation of the option OutputSel has been finished. This
150 %# option can be set by the user to another value than default value.
151 if (~isempty (vodeoptions.OutputSel)), vhaveoutputselection = true;
152 else vhaveoutputselection = false; end
154 %# Implementation of the option OutputSave has been finished. This
155 %# option can be set by the user to another value than default value.
156 if (isempty (vodeoptions.OutputSave)), vodeoptions.OutputSave = 1;
159 %# Implementation of the option Stats has been finished. This option
160 %# can be set by the user to another value than default value.
162 %# Implementation of the option InitialStep has been finished. This
163 %# option can be set by the user to another value than default value.
164 if (isempty (vodeoptions.InitialStep) && ~vstepsizefixed)
165 vodeoptions.InitialStep = (vslot(1,2) - vslot(1,1)) / 10;
166 warning ('OdePkg:InvalidArgument', ...
167 'Option "InitialStep" not set, new value %f is used', vodeoptions.InitialStep);
170 %# Implementation of the option MaxNewtonIterations has been finished. This option
171 %# can be set by the user to another value than default value.
172 if isempty (vodeoptions.MaxNewtonIterations)
173 vodeoptions.MaxNewtonIterations = 10;
174 warning ('OdePkg:InvalidArgument', ...
175 'Option "MaxNewtonIterations" not set, new value %f is used', vodeoptions.MaxNewtonIterations);
178 %# Implementation of the option NewtonTol has been finished. This option
179 %# can be set by the user to another value than default value.
180 if isempty (vodeoptions.NewtonTol)
181 vodeoptions.NewtonTol = 1e-7;
182 warning ('OdePkg:InvalidArgument', ...
183 'Option "NewtonTol" not set, new value %f is used', vodeoptions.NewtonTol);
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 = (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 %# Implementation of the option Jacobian has been finished. This option
200 %# can be set by the user to another value than default value.
201 if (~isempty (vodeoptions.Jacobian) && isnumeric (vodeoptions.Jacobian))
202 vhavejachandle = false; vjac = vodeoptions.Jacobian; %# constant jac
203 elseif (isa (vodeoptions.Jacobian, 'function_handle'))
204 vhavejachandle = true; %# jac defined by a function handle
205 else %# no Jacobian - we will use numerical differentiation
206 vhavejachandle = false;
209 %# Implementation of the option Mass has been finished. This option
210 %# can be set by the user to another value than default value.
211 if (~isempty (vodeoptions.Mass) && isnumeric (vodeoptions.Mass))
212 vhavemasshandle = false; vmass = vodeoptions.Mass; %# constant mass
213 elseif (isa (vodeoptions.Mass, 'function_handle'))
214 vhavemasshandle = true; %# mass defined by a function handle
215 else %# no mass matrix - creating a diag-matrix of ones for mass
216 vhavemasshandle = false; vmass = sparse (eye (length (vinit)) );
219 %# Implementation of the option MStateDependence has been finished.
220 %# This option can be set by the user to another value than default
222 if (strcmp (vodeoptions.MStateDependence, 'none'))
223 vmassdependence = false;
224 else vmassdependence = true;
227 %# Other options that are not used by this solver. Print a warning
228 %# message to tell the user that the option(s) is/are ignored.
229 if (~isequal (vodeoptions.NonNegative, vodetemp.NonNegative))
230 warning ('OdePkg:InvalidArgument', ...
231 'Option "NonNegative" will be ignored by this solver');
233 if (~isequal (vodeoptions.Refine, vodetemp.Refine))
234 warning ('OdePkg:InvalidArgument', ...
235 'Option "Refine" will be ignored by this solver');
237 if (~isequal (vodeoptions.JPattern, vodetemp.JPattern))
238 warning ('OdePkg:InvalidArgument', ...
239 'Option "JPattern" will be ignored by this solver');
241 if (~isequal (vodeoptions.Vectorized, vodetemp.Vectorized))
242 warning ('OdePkg:InvalidArgument', ...
243 'Option "Vectorized" will be ignored by this solver');
245 if (~isequal (vodeoptions.MvPattern, vodetemp.MvPattern))
246 warning ('OdePkg:InvalidArgument', ...
247 'Option "MvPattern" will be ignored by this solver');
249 if (~isequal (vodeoptions.MassSingular, vodetemp.MassSingular))
250 warning ('OdePkg:InvalidArgument', ...
251 'Option "MassSingular" will be ignored by this solver');
253 if (~isequal (vodeoptions.InitialSlope, vodetemp.InitialSlope))
254 warning ('OdePkg:InvalidArgument', ...
255 'Option "InitialSlope" will be ignored by this solver');
257 if (~isequal (vodeoptions.MaxOrder, vodetemp.MaxOrder))
258 warning ('OdePkg:InvalidArgument', ...
259 'Option "MaxOrder" will be ignored by this solver');
261 if (~isequal (vodeoptions.BDF, vodetemp.BDF))
262 warning ('OdePkg:InvalidArgument', ...
263 'Option "BDF" will be ignored by this solver');
266 %# Starting the initialisation of the core solver odebwe
267 vtimestamp = vslot(1,1); %# timestamp = start time
268 vtimelength = length (vslot); %# length needed if fixed steps
269 vtimestop = vslot(1,vtimelength); %# stop time = last value
270 vdirection = sign (vtimestop); %# Flag for direction to solve
273 vstepsize = vodeoptions.InitialStep;
274 vminstepsize = (vtimestop - vtimestamp) / (1/eps);
275 else %# If step size is given then use the fixed time steps
276 vstepsize = vslot(1,2) - vslot(1,1);
277 vminstepsize = sign (vstepsize) * eps;
280 vretvaltime = vtimestamp; %# first timestamp output
281 vretvalresult = vinit; %# first solution output
283 %# Initialize the OutputFcn
284 if (vhaveoutputfunction)
285 if (vhaveoutputselection)
286 vretout = vretvalresult(vodeoptions.OutputSel);
288 vretout = vretvalresult;
290 feval (vodeoptions.OutputFcn, vslot.', ...
291 vretout.', 'init', vfunarguments{:});
294 %# Initialize the EventFcn
295 if (vhaveeventfunction)
296 odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
297 vretvalresult.', 'init', vfunarguments{:});
300 %# Initialize parameters and counters
301 vcntloop = 2; vcntcycles = 1; vu = vinit; vcntsave = 2;
302 vunhandledtermination = true; vpow = 1/2; vnpds = 0;
303 vcntiter = 0; vcntnewt = 0; vndecomps = 0; vnlinsols = 0;
305 %# the following option enables the simplified Newton method
306 %# which evaluates the Jacobian only once instead of the
307 %# standard method that updates the Jacobian in each iteration
308 vsimplified = false; % or true
310 %# The solver main loop - stop if the endpoint has been reached
311 while ((vdirection * (vtimestamp) < vdirection * (vtimestop)) && ...
312 (vdirection * (vstepsize) >= vdirection * (vminstepsize)))
314 %# Hit the endpoint of the time slot exactely
315 if ((vtimestamp + vstepsize) > vdirection * vtimestop)
316 vstepsize = vtimestop - vdirection * vtimestamp;
319 %# Run the time integration for each estimator
320 %# from vtimestamp -> vtimestamp+vstepsize
321 for j = 1:vestimators
322 %# Initial value (result of the previous timestep)
324 %# Initial guess for Newton-Raphson
326 %# We do not use a higher order approximation for the
327 %# comparsion, but two steps by the Backward Euler
330 % Initialize the time stepping parameters
331 vthestep = vstepsize / j;
332 vthetime = vtimestamp + i*vthestep;
334 vresnrm = inf (1, vodeoptions.MaxNewtonIterations);
336 %# Start the Newton iteration
337 while (vnewtit < vodeoptions.MaxNewtonIterations) && ...
338 (vresnrm (vnewtit) > vodeoptions.NewtonTol)
340 %# Compute the Jacobian of the non-linear equation,
341 %# that is the matrix pencil of the mass matrix and
342 %# the right-hand-side's Jacobian. Perform a (sparse)
343 %# LU-Decomposition afterwards.
344 if ( (vnewtit==1) || (~vsimplified) )
345 %# Get the mass matrix from the left-hand-side
346 if (vhavemasshandle) %# Handle only the dynamic mass matrix,
347 if (vmassdependence) %# constant mass matrices have already
348 vmass = feval ... %# been set before (if any)
349 (vodeoptions.Mass, vthetime, y(j,:)', vfunarguments{:});
350 else %# if (vmassdependence == false)
351 vmass = feval ... %# then we only have the time argument
352 (vodeoptions.Mass, y(j,:)', vfunarguments{:});
355 %# Get the Jacobian of the right-hand-side's function
356 if (vhavejachandle) %# Handle only the dynamic jacobian
357 vjac = feval(vodeoptions.Jacobian, vthetime,...
358 y(j,:)', vfunarguments{:});
359 elseif isempty(vodeoptions.Jacobian) %# If no Jacobian is given
360 vjac = feval(@jacobian, vfun, vthetime,y(j,:)',...
361 vfunarguments); %# then we differentiate
364 vfulljac = vmass/vthestep - vjac;
365 %# one could do a matrix decomposition of vfulljac here,
366 %# but the choice of decomposition depends on the problem
367 %# and therefore we use the backslash-operator in row 374
370 %# Compute the residual
371 vres = vmass/vthestep*(y(j,:)-y0)' - feval(vfun,vthetime,y(j,:)',vfunarguments{:});
372 vresnrm(vnewtit+1) = norm(vres,inf);
373 %# Solve the linear system
374 y(j,:) = vfulljac\(-vres+vfulljac*y(j,:)');
375 %# the backslash operator decomposes the matrix
376 %# and solves the system in a single step.
377 vndecomps = vndecomps + 1;
378 vnlinsols = vnlinsols + 1;
379 %# Prepare next iteration
380 vnewtit = vnewtit + 1;
383 %# Leave inner loop if Newton diverged
384 if vresnrm(vnewtit)>vodeoptions.NewtonTol
387 %# Save intermediate solution as initial value
388 %# for the next intermediate step
390 %# Count all Newton iterations
391 vcntnewt = vcntnewt + (vnewtit-1);
394 %# Leave outer loop if Newton diverged
395 if vresnrm(vnewtit)>vodeoptions.NewtonTol
398 end %# for estimators
400 % if all Newton iterations converged
401 if vresnrm(vnewtit)<vodeoptions.NewtonTol
402 %# First order approximation using step size h
404 %# If adaptive: first order approximation using step
405 %# size h/2, if fixed: y1=y2=y3
406 y2 = y(vestimators,:);
407 %# Second order approximation by ("Richardson")
408 %# extrapolation using h and h/2
412 %# If Newton did not converge, repeat step with reduced
413 %# step size, otherwise calculate the absolute local
414 %# truncation error and the acceptable error
415 if vresnrm(vnewtit)>vodeoptions.NewtonTol
416 vdelta = 2; vtau = 1;
417 elseif (~vstepsizefixed)
419 vdelta = abs (y3 - y1)';
420 vtau = max (vodeoptions.RelTol * abs (vu.'), vodeoptions.AbsTol);
422 vdelta = norm ((y3 - y1)', Inf);
423 vtau = max (vodeoptions.RelTol * max (norm (vu.', Inf), 1.0), ...
426 else %# if (vstepsizefixed == true)
427 vdelta = 1; vtau = 2;
430 %# If the error is acceptable then update the vretval variables
431 if (all (vdelta <= vtau))
432 vtimestamp = vtimestamp + vstepsize;
433 vu = y2; % or y3 if we want the extrapolation....
435 %# Save the solution every vodeoptions.OutputSave steps
436 if (mod (vcntloop-1,vodeoptions.OutputSave) == 0)
437 vretvaltime(vcntsave,:) = vtimestamp;
438 vretvalresult(vcntsave,:) = vu;
439 vcntsave = vcntsave + 1;
441 vcntloop = vcntloop + 1; vcntiter = 0;
443 %# Call plot only if a valid result has been found, therefore this
444 %# code fragment has moved here. Stop integration if plot function
446 if (vhaveoutputfunction)
447 if (vhaveoutputselection)
448 vpltout = vu(vodeoptions.OutputSel);
452 vpltret = feval (vodeoptions.OutputFcn, vtimestamp, ...
453 vpltout.', [], vfunarguments{:});
454 if vpltret %# Leave loop
455 vunhandledtermination = false; break;
459 %# Call event only if a valid result has been found, therefore this
460 %# code fragment has moved here. Stop integration if veventbreak is
462 if (vhaveeventfunction)
464 odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
465 vu(:), [], vfunarguments{:});
466 if (~isempty (vevent{1}) && vevent{1} == 1)
467 vretvaltime(vcntloop-1,:) = vevent{3}(end,:);
468 vretvalresult(vcntloop-1,:) = vevent{4}(end,:);
469 vunhandledtermination = false; break;
472 end %# If the error is acceptable ...
474 %# Update the step size for the next integration step
476 %# 20080425, reported by Marco Caliari
477 %# vdelta cannot be negative (because of the absolute value that
478 %# has been introduced) but it could be 0, then replace the zeros
479 %# with the maximum value of vdelta
480 vdelta(find (vdelta == 0)) = max (vdelta);
481 %# It could happen that max (vdelta) == 0 (ie. that the original
482 %# vdelta was 0), in that case we double the previous vstepsize
483 vdelta(find (vdelta == 0)) = max (vtau) .* (0.4 ^ (1 / vpow));
486 vstepsize = min (vodeoptions.MaxStep, ...
487 min (0.8 * vstepsize * (vtau ./ vdelta) .^ vpow));
489 vstepsize = max (vodeoptions.MaxStep, ...
490 max (0.8 * vstepsize * (vtau ./ vdelta) .^ vpow));
493 else %# if (vstepsizefixed)
494 if (vresnrm(vnewtit)>vodeoptions.NewtonTol)
495 vunhandledtermination = true;
497 elseif (vcntloop <= vtimelength)
498 vstepsize = vslot(vcntloop) - vslot(vcntloop-1);
499 else %# Get out of the main integration loop
504 %# Update counters that count the number of iteration cycles
505 vcntcycles = vcntcycles + 1; %# Needed for cost statistics
506 vcntiter = vcntiter + 1; %# Needed to find iteration problems
508 %# Stop solving because the last 1000 steps no successful valid
509 %# value has been found
510 if (vcntiter >= 5000)
511 error (['Solving has not been successful. The iterative', ...
512 ' integration loop exited at time t = %f before endpoint at', ...
513 ' tend = %f was reached. This happened because the iterative', ...
514 ' integration loop does not find a valid solution at this time', ...
515 ' stamp. Try to reduce the value of "InitialStep" and/or', ...
516 ' "MaxStep" with the command "odeset".\n'], vtimestamp, vtimestop);
521 %# Check if integration of the ode has been successful
522 if (vdirection * vtimestamp < vdirection * vtimestop)
523 if (vunhandledtermination == true)
524 error ('OdePkg:InvalidArgument', ...
525 ['Solving has not been successful. The iterative', ...
526 ' integration loop exited at time t = %f', ...
527 ' before endpoint at tend = %f was reached. This may', ...
528 ' happen if the stepsize grows smaller than defined in', ...
529 ' vminstepsize. Try to reduce the value of "InitialStep" and/or', ...
530 ' "MaxStep" with the command "odeset".\n'], vtimestamp, vtimestop);
532 warning ('OdePkg:InvalidArgument', ...
533 ['Solver has been stopped by a call of "break" in', ...
534 ' the main iteration loop at time t = %f before endpoint at', ...
535 ' tend = %f was reached. This may happen because the @odeplot', ...
536 ' function returned "true" or the @event function returned "true".'], ...
537 vtimestamp, vtimestop);
541 %# Postprocessing, do whatever when terminating integration algorithm
542 if (vhaveoutputfunction) %# Cleanup plotter
543 feval (vodeoptions.OutputFcn, vtimestamp, ...
544 vu.', 'done', vfunarguments{:});
546 if (vhaveeventfunction) %# Cleanup event function handling
547 odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
548 vu.', 'done', vfunarguments{:});
550 %# Save the last step, if not already saved
551 if (mod (vcntloop-2,vodeoptions.OutputSave) ~= 0)
552 vretvaltime(vcntsave,:) = vtimestamp;
553 vretvalresult(vcntsave,:) = vu;
556 %# Print additional information if option Stats is set
557 if (strcmp (vodeoptions.Stats, 'on'))
559 vnsteps = vcntloop-2; %# vcntloop from 2..end
560 vnfailed = (vcntcycles-1)-(vcntloop-2)+1; %# vcntcycl from 1..end
561 vnfevals = vcntnewt; %# number of rhs evaluations
562 if isempty(vodeoptions.Jacobian) %# additional evaluations due
563 vnfevals = vnfevals + vcntnewt*(1+length(vinit)); %# to differentiation
565 %# Print cost statistics if no output argument is given
567 vmsg = fprintf (1, 'Number of successful steps: %d\n', vnsteps);
568 vmsg = fprintf (1, 'Number of failed attempts: %d\n', vnfailed);
569 vmsg = fprintf (1, 'Number of function calls: %d\n', vnfevals);
575 if (nargout == 1) %# Sort output variables, depends on nargout
576 varargout{1}.x = vretvaltime; %# Time stamps are saved in field x
577 varargout{1}.y = vretvalresult; %# Results are saved in field y
578 varargout{1}.solver = 'odebwe'; %# Solver name is saved in field solver
579 if (vhaveeventfunction)
580 varargout{1}.ie = vevent{2}; %# Index info which event occured
581 varargout{1}.xe = vevent{3}; %# Time info when an event occured
582 varargout{1}.ye = vevent{4}; %# Results when an event occured
585 varargout{1}.stats = struct;
586 varargout{1}.stats.nsteps = vnsteps;
587 varargout{1}.stats.nfailed = vnfailed;
588 varargout{1}.stats.nfevals = vnfevals;
589 varargout{1}.stats.npds = vnpds;
590 varargout{1}.stats.ndecomps = vndecomps;
591 varargout{1}.stats.nlinsols = vnlinsols;
593 elseif (nargout == 2)
594 varargout{1} = vretvaltime; %# Time stamps are first output argument
595 varargout{2} = vretvalresult; %# Results are second output argument
596 elseif (nargout == 5)
597 varargout{1} = vretvaltime; %# Same as (nargout == 2)
598 varargout{2} = vretvalresult; %# Same as (nargout == 2)
599 varargout{3} = []; %# LabMat doesn't accept lines like
600 varargout{4} = []; %# varargout{3} = varargout{4} = [];
602 if (vhaveeventfunction)
603 varargout{3} = vevent{3}; %# Time info when an event occured
604 varargout{4} = vevent{4}; %# Results when an event occured
605 varargout{5} = vevent{2}; %# Index info which event occured
610 function df = jacobian(vfun,vthetime,vtheinput,vfunarguments);
611 %# Internal function for the numerical approximation of a jacobian
612 vlen = length(vtheinput);
614 vfun0 = feval(vfun,vthetime,vtheinput,vfunarguments{:});
615 df = zeros(vlen,vlen);
617 vbuffer = vtheinput(j);
620 elseif (abs(vbuffer)>1)
623 vh = sign(vbuffer)*vsigma;
625 vtheinput(j) = vbuffer + vh;
626 df(:,j) = (feval(vfun,vthetime,vtheinput,...
627 vfunarguments{:}) - vfun0) / vh;
628 vtheinput(j) = vbuffer;
631 %! # We are using the "Van der Pol" implementation for all tests that
632 %! # are done for this function. We also define a Jacobian, Events,
633 %! # pseudo-Mass implementation. For further tests we also define a
634 %! # reference solution (computed at high accuracy) and an OutputFcn
635 %!function [ydot] = fpol (vt, vy, varargin) %# The Van der Pol
636 %! ydot = [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)];
637 %!function [vjac] = fjac (vt, vy, varargin) %# its Jacobian
638 %! vjac = [0, 1; -1 - 2 * vy(1) * vy(2), 1 - vy(1)^2];
639 %!function [vjac] = fjcc (vt, vy, varargin) %# sparse type
640 %! vjac = sparse ([0, 1; -1 - 2 * vy(1) * vy(2), 1 - vy(1)^2]);
641 %!function [vval, vtrm, vdir] = feve (vt, vy, varargin)
642 %! vval = fpol (vt, vy, varargin); %# We use the derivatives
643 %! vtrm = zeros (2,1); %# that's why component 2
644 %! vdir = ones (2,1); %# seems to not be exact
645 %!function [vval, vtrm, vdir] = fevn (vt, vy, varargin)
646 %! vval = fpol (vt, vy, varargin); %# We use the derivatives
647 %! vtrm = ones (2,1); %# that's why component 2
648 %! vdir = ones (2,1); %# seems to not be exact
649 %!function [vmas] = fmas (vt, vy)
650 %! vmas = [1, 0; 0, 1]; %# Dummy mass matrix for tests
651 %!function [vmas] = fmsa (vt, vy)
652 %! vmas = sparse ([1, 0; 0, 1]); %# A sparse dummy matrix
653 %!function [vref] = fref () %# The computed reference sol
654 %! vref = [0.32331666704577, -1.83297456798624];
655 %!function [vout] = fout (vt, vy, vflag, varargin)
656 %! if (regexp (char (vflag), 'init') == 1)
657 %! if (any (size (vt) ~= [2, 1])) error ('"fout" step "init"'); end
658 %! elseif (isempty (vflag))
659 %! if (any (size (vt) ~= [1, 1])) error ('"fout" step "calc"'); end
661 %! elseif (regexp (char (vflag), 'done') == 1)
662 %! if (any (size (vt) ~= [1, 1])) error ('"fout" step "done"'); end
663 %! else error ('"fout" invalid vflag');
666 %! %# Turn off output of warning messages for all tests, turn them on
667 %! %# again if the last test is called
668 %!error %# input argument number one
669 %! warning ('off', 'OdePkg:InvalidArgument');
670 %! B = odebwe (1, [0 25], [3 15 1]);
671 %!error %# input argument number two
672 %! B = odebwe (@fpol, 1, [3 15 1]);
673 %!error %# input argument number three
674 %! B = odebwe (@flor, [0 25], 1);
675 %!test %# one output argument
676 %! vsol = odebwe (@fpol, [0 2], [2 0]);
677 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
678 %! assert (isfield (vsol, 'solver'));
679 %! assert (vsol.solver, 'odebwe');
680 %!test %# two output arguments
681 %! [vt, vy] = odebwe (@fpol, [0 2], [2 0]);
682 %! assert ([vt(end), vy(end,:)], [2, fref], 1e-3);
683 %!test %# five output arguments and no Events
684 %! [vt, vy, vxe, vye, vie] = odebwe (@fpol, [0 2], [2 0]);
685 %! assert ([vt(end), vy(end,:)], [2, fref], 1e-3);
686 %! assert ([vie, vxe, vye], []);
687 %!test %# anonymous function instead of real function
688 %! fvdb = @(vt,vy) [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)];
689 %! vsol = odebwe (fvdb, [0 2], [2 0]);
690 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
691 %!test %# extra input arguments passed trhough
692 %! vsol = odebwe (@fpol, [0 2], [2 0], 12, 13, 'KL');
693 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
694 %!test %# empty OdePkg structure *but* extra input arguments
696 %! vsol = odebwe (@fpol, [0 2], [2 0], vopt, 12, 13, 'KL');
697 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
698 %!error %# strange OdePkg structure
699 %! vopt = struct ('foo', 1);
700 %! vsol = odebwe (@fpol, [0 2], [2 0], vopt);
701 %!test %# Solve vdp in fixed step sizes
702 %! vsol = odebwe (@fpol, [0:0.001:2], [2 0]);
703 %! assert (vsol.x(:), [0:0.001:2]');
704 %! assert (vsol.y(end,:), fref, 1e-2);
705 %!test %# Solve in backward direction starting at t=0
706 %! %# vref = [-1.2054034414, 0.9514292694];
707 %! vsol = odebwe (@fpol, [0 -2], [2 0]);
708 %! %# assert ([vsol.x(end), vsol.y(end,:)], [-2, fref], 1e-3);
709 %!test %# Solve in backward direction starting at t=2
710 %! %# vref = [-1.2154183302, 0.9433018000];
711 %! vsol = odebwe (@fpol, [2 -2], [0.3233166627 -1.8329746843]);
712 %! %# assert ([vsol.x(end), vsol.y(end,:)], [-2, fref], 1e-3);
713 %!test %# AbsTol option
714 %! vopt = odeset ('AbsTol', 1e-5);
715 %! vsol = odebwe (@fpol, [0 2], [2 0], vopt);
716 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-2);
717 %!test %# AbsTol and RelTol option
718 %! vopt = odeset ('AbsTol', 1e-6, 'RelTol', 1e-6);
719 %! vsol = odebwe (@fpol, [0 2], [2 0], vopt);
720 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-2);
721 %!test %# RelTol and NormControl option -- higher accuracy
722 %! vopt = odeset ('RelTol', 1e-6, 'NormControl', 'on');
723 %! vsol = odebwe (@fpol, [0 2], [2 0], vopt);
724 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
726 %! %# test for NonNegative option is missing
727 %! %# test for OutputSel and Refine option is missing
729 %!test %# Details of OutputSave can't be tested
730 %! vopt = odeset ('OutputSave', 1, 'OutputSel', 1);
731 %! vsla = odebwe (@fpol, [0 2], [2 0], vopt);
732 %! vopt = odeset ('OutputSave', 2);
733 %! vslb = odebwe (@fpol, [0 2], [2 0], vopt);
734 %! assert (length (vsla.x) > length (vslb.x))
735 %!test %# Stats must add further elements in vsol
736 %! vopt = odeset ('Stats', 'on');
737 %! vsol = odebwe (@fpol, [0 2], [2 0], vopt);
738 %! assert (isfield (vsol, 'stats'));
739 %! assert (isfield (vsol.stats, 'nsteps'));
740 %!test %# InitialStep option
741 %! vopt = odeset ('InitialStep', 1e-8);
742 %! vsol = odebwe (@fpol, [0 0.2], [2 0], vopt);
743 %! assert ([vsol.x(2)-vsol.x(1)], [1e-8], 1e-9);
744 %!test %# MaxStep option
745 %! vopt = odeset ('MaxStep', 1e-2);
746 %! vsol = odebwe (@fpol, [0 0.2], [2 0], vopt);
747 %! assert ([vsol.x(5)-vsol.x(4)], [1e-2], 1e-2);
748 %!test %# Events option add further elements in vsol
749 %! vopt = odeset ('Events', @feve);
750 %! vsol = odebwe (@fpol, [0 10], [2 0], vopt);
751 %! assert (isfield (vsol, 'ie'));
752 %! assert (vsol.ie, [2; 1; 2; 1]);
753 %! assert (isfield (vsol, 'xe'));
754 %! assert (isfield (vsol, 'ye'));
755 %!test %# Events option, now stop integration
756 %! vopt = odeset ('Events', @fevn, 'NormControl', 'on');
757 %! vsol = odebwe (@fpol, [0 10], [2 0], vopt);
758 %! assert ([vsol.ie, vsol.xe, vsol.ye], ...
759 %! [2.0, 2.496110, -0.830550, -2.677589], 1e-3);
760 %!test %# Events option, five output arguments
761 %! vopt = odeset ('Events', @fevn, 'NormControl', 'on');
762 %! [vt, vy, vxe, vye, vie] = odebwe (@fpol, [0 10], [2 0], vopt);
763 %! assert ([vie, vxe, vye], ...
764 %! [2.0, 2.496110, -0.830550, -2.677589], 1e-3);
765 %!test %# Jacobian option
766 %! vopt = odeset ('Jacobian', @fjac);
767 %! vsol = odebwe (@fpol, [0 2], [2 0], vopt);
768 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
769 %!test %# Jacobian option and sparse return value
770 %! vopt = odeset ('Jacobian', @fjcc);
771 %! vsol = odebwe (@fpol, [0 2], [2 0], vopt);
772 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
774 %! %# test for JPattern option is missing
775 %! %# test for Vectorized option is missing
777 %!test %# Mass option as function
778 %! vopt = odeset ('Mass', @fmas);
779 %! vsol = odebwe (@fpol, [0 2], [2 0], vopt);
780 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
781 %!test %# Mass option as matrix
782 %! vopt = odeset ('Mass', eye (2,2));
783 %! vsol = odebwe (@fpol, [0 2], [2 0], vopt);
784 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
785 %!test %# Mass option as sparse matrix
786 %! vopt = odeset ('Mass', sparse (eye (2,2)));
787 %! vsol = odebwe (@fpol, [0 2], [2 0], vopt);
788 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
789 %!test %# Mass option as function and sparse matrix
790 %! vopt = odeset ('Mass', @fmsa);
791 %! vsol = odebwe (@fpol, [0 2], [2 0], vopt);
792 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
794 %! %# test for MStateDependence option is missing
795 %! %# test for MvPattern option is missing
796 %! %# test for InitialSlope option is missing
797 %! %# test for MaxOrder option is missing
798 %! %# test for BDF option is missing
800 %! warning ('on', 'OdePkg:InvalidArgument');
802 %# Local Variables: ***