Mercurial > octave-nkf

comparison scripts/ode/ode45.m @ 20568:fcb792acab9b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Moving ode45, odeset, odeget, and levenshtein from odepkg to core. * libinterp/corefcn/levenshtein.cc: move function from odepkg into core * libinterp/corefcn/module.mk: include levenshtein.cc * scripts/ode: move ode45, odeset, odeget, and all dependencies from odepkg into core * scripts/module.mk: include them * doc/interpreter/diffeq.txi: add documentation for ode45, odeset, odeget * NEWS: announce functions included with this changeset * scripts/help/__unimplemented__.m: removed new functions
author jcorno <jacopo.corno@gmail.com>
date Thu, 24 Sep 2015 12:58:46 +0200
parents
children 3339c9bdfe6a
comparison
equal deleted inserted replaced
20567:2480bbcd1333 20568:fcb792acab9b
1 ## Copyright (C) 2006-2012, Thomas Treichl <treichl@users.sourceforge.net>
2 ## Copyright (C) 2013, Roberto Porcu' <roberto.porcu@polimi.it>
3 ## Copyright (C) 2014, Jacopo Corno <jacopo.corno@gmail.com>
4 ##
5 ## This file is part of Octave.
6 ##
7 ## Octave is free software; you can redistribute it and/or modify it
8 ## under the terms of the GNU General Public License as published by
9 ## the Free Software Foundation; either version 3 of the License, or (at
10 ## your option) any later version.
11 ##
12 ## Octave is distributed in the hope that it will be useful, but
13 ## WITHOUT ANY WARRANTY; without even the implied warranty of
14 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 ## General Public License for more details.
16 ##
17 ## You should have received a copy of the GNU General Public License
18 ## along with Octave; see the file COPYING. If not, see
19 ## <http://www.gnu.org/licenses/>.
20
21 ## -*- texinfo -*-
22 ## @deftypefn {Function File} {[@var{sol}] =} ode45 (@var{fun}, @var{slot}, @var{init}, [@var{opt}], [@var{par1}, @var{par2}, @dots{}])
23 ## @deftypefnx {Function File} {[@var{t}, @var{y}, [@var{xe}, @var{ye}, @var{ie}]] =} ode45 (@var{fun}, @var{slot}, @var{init}, [@var{opt}], [@var{par1}, @var{par2}, @dots{}])
24 ##
25 ## This function file can be used to solve a set of non--stiff ordinary
26 ## differential equations (non--stiff ODEs) with the well known explicit
27 ## Dormand-Prince method of order 4.
28 ##
29 ## This function can be called with two output arguments: @var{t} and @var{y}.
30 ## Variable @var{t} is a column vector and contains the time stamps, instead
31 ## @var{y} is a matrix in which each column refers to a different unknown of
32 ## the problem and the rows number is the same of @var{t} rows number so that each
33 ## row of @var{y} contains the values of all unknowns at the time value contained
34 ## in the corresponding row in @var{t}.
35 ##
36 ## The first input argument must be a function_handle or an inline function that
37 ## defines the set of ODE: @code{y' = f(t,y)}. As described above, this function
38 ## must take two input arguments, where the first is the time and the second
39 ## the unknowns, and must have just one output argument.
40 ##
41 ## The second input argument must contain time informations. Usually it should
42 ## be a vector with at least two elements which define the initial and the final
43 ## time instants; if the elements are more than two, then the solution will be
44 ## evaluated also at these intermediate time instants unless the integrate function
45 ## called is the @command{integrate_n_steps}. If there is only one time value,
46 ## then it will give an error unless the options structure has no empty fields
47 ## named @var{"TimeStepNumber"} and @var{"TimeStepSize"}. If the option
48 ## @var{"TimeStepSize"} is not empty, then the stepper called will be
49 ## @command{integrate_const}, if also @var{"TimeStepNumber"} is not empty it will
50 ## be called the integrate function @command{integrate_n_steps}, otherwise it will
51 ## be called @command{integrate_adaptive}. For this last possibility the user can
52 ## set the tolerance for the timestep computation by setting a value to the option
53 ## @var{"Tau"}, that as default value has @math{1.e-6}.
54 ##
55 ## The third input argument must contain the initial value for the unknown.
56 ## If this is a vector then the solution @var{y} will be a matrix in which each
57 ## column is the solution for the corresponding initial value in @var{init}.
58 ##
59 ## The fourth input argument is not mandatory and it should contain a structure
60 ## with valid ODE fields.
61 ##
62 ## For example, solve an anonymous implementation of the Van der Pol equation
63 ## @example
64 ## fvdp = @@(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)];
65 ## [T,Y] = ode45 (fvdp, [0 20], [2 0]);
66 ## @end example
67 ## @end deftypefn
68 ##
69
70 function [varargout] = ode45 (vfun, vslot, vinit, varargin)
71
72 vorder = 5; % runge_kutta_45_dorpri uses local extrapolation
73 vsolver = "ode45";
74
75 if (nargin == 0) ## Check number and types of all input arguments
76 help (vsolver);
77 error ("OdePkg:InvalidArgument", ...
78 "number of input arguments must be greater than zero");
79 endif
80
81 if (nargin < 3)
82 print_usage;
83 endif
84
85 if (nargin >= 4)
86 if (~isstruct (varargin{1}))
87 ## varargin{1:len} are parameters for vfun
88 vodeoptions = odeset;
89 vodeoptions.vfunarguments = varargin;
90 elseif (length (varargin) > 1)
91 ## varargin{1} is an OdePkg options structure vopt
92 vodeoptions = odepkg_structure_check (varargin{1}, "ode45");
93 vodeoptions.vfunarguments = {varargin{2:length(varargin)}};
94 else ## if (isstruct (varargin{1}))
95 vodeoptions = odepkg_structure_check (varargin{1}, "ode45");
96 vodeoptions.vfunarguments = {};
97 endif
98 else ## if (nargin == 3)
99 vodeoptions = odeset;
100 vodeoptions.vfunarguments = {};
101 endif
102
103 if (~isvector (vslot) || ~isnumeric (vslot))
104 error ("OdePkg:InvalidArgument", ...
105 "second input argument must be a valid vector");
106 endif
107
108 if (length (vslot) < 2 && ...
109 (isempty (vodeoptions.TimeStepSize) ...
110 || isempty (vodeoptions.TimeStepNumber)))
111 error ("OdePkg:InvalidArgument", ...
112 "second input argument must be a valid vector");
113 elseif (vslot(2) == vslot(1))
114 error ("OdePkg:InvalidArgument", ...
115 "second input argument must be a valid vector");
116 else
117 vodeoptions.vdirection = sign (vslot(2) - vslot(1));
118 endif
119 vslot = vslot(:);
120
121 if (~isvector (vinit) || ~isnumeric (vinit))
122 error ("OdePkg:InvalidArgument", ...
123 "third input argument must be a valid numerical value");
124 endif
125 vinit = vinit(:);
126
127 if ~(isa (vfun, "function_handle") || isa (vfun, "inline"))
128 error ("OdePkg:InvalidArgument", ...
129 "first input argument must be a valid function handle");
130 endif
131
132 ## Start preprocessing, have a look which options are set in
133 ## vodeoptions, check if an invalid or unused option is set
134 if (isempty (vodeoptions.TimeStepNumber) ...
135 && isempty (vodeoptions.TimeStepSize))
136 integrate_func = "adaptive";
137 vodeoptions.vstepsizefixed = false;
138 elseif (~isempty (vodeoptions.TimeStepNumber) ...
139 && ~isempty (vodeoptions.TimeStepSize))
140 integrate_func = "n_steps";
141 vodeoptions.vstepsizefixed = true;
142 if (sign (vodeoptions.TimeStepSize) != vodeoptions.vdirection)
143 warning ("OdePkg:InvalidArgument", ...
144 "option ''TimeStepSize'' has a wrong sign", ...
145 "it will be corrected automatically");
146 vodeoptions.TimeStepSize = (-1)*vodeoptions.TimeStepSize;
147 endif
148 elseif (isempty (vodeoptions.TimeStepNumber) && ~isempty (vodeoptions.TimeStepSize))
149 integrate_func = "const";
150 vodeoptions.vstepsizefixed = true;
151 if (sign (vodeoptions.TimeStepSize) != vodeoptions.vdirection)
152 warning ("OdePkg:InvalidArgument", ...
153 "option ''TimeStepSize'' has a wrong sign", ...
154 "it will be corrected automatically");
155 vodeoptions.TimeStepSize = (-1)*vodeoptions.TimeStepSize;
156 endif
157 else
158 warning ("OdePkg:InvalidArgument", ...
159 "assuming an adaptive integrate function");
160 integrate_func = "adaptive";
161 endif
162
163 ## Get the default options that can be set with "odeset" temporarily
164 vodetemp = odeset;
165
166 ## Implementation of the option RelTol has been finished. This option
167 ## can be set by the user to another value than default value.
168 if (isempty (vodeoptions.RelTol) && ~vodeoptions.vstepsizefixed)
169 vodeoptions.RelTol = 1e-3;
170 warning ("OdePkg:InvalidArgument", ...
171 "Option ''RelTol'' not set, new value %f is used", ...
172 vodeoptions.RelTol);
173 elseif (~isempty (vodeoptions.RelTol) && vodeoptions.vstepsizefixed)
174 warning ("OdePkg:InvalidArgument", ...
175 "Option ''RelTol'' will be ignored if fixed time stamps are given");
176 endif
177
178 ## Implementation of the option AbsTol has been finished. This option
179 ## can be set by the user to another value than default value.
180 if (isempty (vodeoptions.AbsTol) && ~vodeoptions.vstepsizefixed)
181 vodeoptions.AbsTol = 1e-6;
182 warning ("OdePkg:InvalidArgument", ...
183 "Option ''AbsTol'' not set, new value %f is used", ...
184 vodeoptions.AbsTol);
185 elseif (~isempty (vodeoptions.AbsTol) && vodeoptions.vstepsizefixed)
186 warning ("OdePkg:InvalidArgument", ...
187 "Option ''AbsTol'' will be ignored if fixed time stamps are given");
188 else
189 vodeoptions.AbsTol = vodeoptions.AbsTol(:); ## Create column vector
190 endif
191
192 ## Implementation of the option NormControl has been finished. This
193 ## option can be set by the user to another value than default value.
194 if (strcmp (vodeoptions.NormControl, "on"))
195 vodeoptions.vnormcontrol = true;
196 else
197 vodeoptions.vnormcontrol = false;
198 endif
199
200 ## Implementation of the option NonNegative has been finished. This
201 ## option can be set by the user to another value than default value.
202 if (~isempty (vodeoptions.NonNegative))
203 if (isempty (vodeoptions.Mass))
204 vodeoptions.vhavenonnegative = true;
205 else
206 vodeoptions.vhavenonnegative = false;
207 warning ("OdePkg:InvalidArgument", ...
208 "Option 'NonNegative' will be ignored if mass matrix is set");
209 endif
210 else
211 vodeoptions.vhavenonnegative = false;
212 endif
213
214 ## Implementation of the option OutputFcn has been finished. This
215 ## option can be set by the user to another value than default value.
216 if (isempty (vodeoptions.OutputFcn) && nargout == 0)
217 vodeoptions.OutputFcn = @odeplot;
218 vodeoptions.vhaveoutputfunction = true;
219 elseif (isempty (vodeoptions.OutputFcn))
220 vodeoptions.vhaveoutputfunction = false;
221 else
222 vodeoptions.vhaveoutputfunction = true;
223 endif
224
225 ## Implementation of the option OutputSel has been finished. This
226 ## option can be set by the user to another value than default value.
227 if (~isempty (vodeoptions.OutputSel))
228 vodeoptions.vhaveoutputselection = true;
229 else
230 vodeoptions.vhaveoutputselection = false;
231 endif
232
233 ## Implementation of the option OutputSave has been finished. This
234 ## option can be set by the user to another value than default value.
235 if (isempty (vodeoptions.OutputSave))
236 vodeoptions.OutputSave = 1;
237 endif
238
239 ## Implementation of the option Refine has been finished. This option
240 ## can be set by the user to another value than default value.
241 if (vodeoptions.Refine > 0)
242 vodeoptions.vhaverefine = true;
243 else
244 vodeoptions.vhaverefine = false;
245 endif
246
247 ## Implementation of the option Stats has been finished. This option
248 ## can be set by the user to another value than default value.
249
250 ## Implementation of the option InitialStep has been finished. This
251 ## option can be set by the user to another value than default value.
252 if (isempty (vodeoptions.InitialStep) && strcmp (integrate_func, "adaptive"))
253 vodeoptions.InitialStep = vodeoptions.vdirection* ...
254 starting_stepsize (vorder, vfun, vslot(1), vinit, vodeoptions.AbsTol, ...
255 vodeoptions.RelTol, vodeoptions.vnormcontrol);
256 warning ("OdePkg:InvalidArgument", ...
257 "option ''InitialStep'' not set, estimated value %f is used", ...
258 vodeoptions.InitialStep);
259 elseif(isempty (vodeoptions.InitialStep))
260 vodeoptions.InitialStep = odeget (vodeoptions, "TimeStepSize");
261 endif
262
263 ## Implementation of the option MaxStep has been finished. This option
264 ## can be set by the user to another value than default value.
265 if (isempty (vodeoptions.MaxStep) && ~vodeoptions.vstepsizefixed)
266 vodeoptions.MaxStep = abs (vslot(end) - vslot(1)) / 10;
267 warning ("OdePkg:InvalidArgument", ...
268 "Option ''MaxStep'' not set, new value %f is used", ...
269 vodeoptions.MaxStep);
270 endif
271
272 ## Implementation of the option Events has been finished. This option
273 ## can be set by the user to another value than default value.
274 if (~isempty (vodeoptions.Events))
275 vodeoptions.vhaveeventfunction = true;
276 else
277 vodeoptions.vhaveeventfunction = false;
278 endif
279
280 ## The options 'Jacobian', 'JPattern' and 'Vectorized' will be ignored
281 ## by this solver because this solver uses an explicit Runge-Kutta
282 ## method and therefore no Jacobian calculation is necessary
283 if (~isequal (vodeoptions.Jacobian, vodetemp.Jacobian))
284 warning ("OdePkg:InvalidArgument", ...
285 "option ''Jacobian'' will be ignored by this solver");
286 endif
287
288 if (~isequal (vodeoptions.JPattern, vodetemp.JPattern))
289 warning ("OdePkg:InvalidArgument", ...
290 "option ''JPattern'' will be ignored by this solver");
291 endif
292
293 if (~isequal (vodeoptions.Vectorized, vodetemp.Vectorized))
294 warning ("OdePkg:InvalidArgument", ...
295 "option ''Vectorized'' will be ignored by this solver");
296 endif
297
298 if (~isequal (vodeoptions.NewtonTol, vodetemp.NewtonTol))
299 warning ("OdePkg:InvalidArgument", ...
300 "option ''NewtonTol'' will be ignored by this solver");
301 endif
302
303 if (~isequal (vodeoptions.MaxNewtonIterations,vodetemp.MaxNewtonIterations))
304 warning ("OdePkg:InvalidArgument", ...
305 "option ''MaxNewtonIterations'' will be ignored by this solver");
306 endif
307
308 ## Implementation of the option Mass has been finished. This option
309 ## can be set by the user to another value than default value.
310 if (~isempty (vodeoptions.Mass) && isnumeric (vodeoptions.Mass))
311 vhavemasshandle = false;
312 vmass = vodeoptions.Mass; ## constant mass
313 elseif (isa (vodeoptions.Mass, "function_handle"))
314 vhavemasshandle = true; ## mass defined by a function handle
315 else ## no mass matrix - creating a diag-matrix of ones for mass
316 vhavemasshandle = false; ## vmass = diag (ones (length (vinit), 1), 0);
317 endif
318
319 ## Implementation of the option MStateDependence has been finished.
320 ## This option can be set by the user to another value than default
321 ## value.
322 if (strcmp (vodeoptions.MStateDependence, "none"))
323 vmassdependence = false;
324 else
325 vmassdependence = true;
326 endif
327
328 ## Other options that are not used by this solver. Print a warning
329 ## message to tell the user that the option(s) is/are ignored.
330 if (~isequal (vodeoptions.MvPattern, vodetemp.MvPattern))
331 warning ("OdePkg:InvalidArgument", ...
332 "option ''MvPattern'' will be ignored by this solver");
333 endif
334
335 if (~isequal (vodeoptions.MassSingular, vodetemp.MassSingular))
336 warning ("OdePkg:InvalidArgument", ...
337 "option ''MassSingular'' will be ignored by this solver");
338 endif
339
340 if (~isequal (vodeoptions.InitialSlope, vodetemp.InitialSlope))
341 warning ("OdePkg:InvalidArgument", ...
342 "option ''InitialSlope'' will be ignored by this solver");
343 endif
344
345 if (~isequal (vodeoptions.MaxOrder, vodetemp.MaxOrder))
346 warning ("OdePkg:InvalidArgument", ...
347 "option ''MaxOrder'' will be ignored by this solver");
348 endif
349
350 if (~isequal (vodeoptions.BDF, vodetemp.BDF))
351 warning ("OdePkg:InvalidArgument", ...
352 "option ''BDF'' will be ignored by this solver");
353 endif
354
355 ## Starting the initialisation of the core solver ode45
356 SubOpts = vodeoptions;
357
358 if (vhavemasshandle) ## Handle only the dynamic mass matrix,
359 if (vmassdependence) ## constant mass matrices have already
360 vmass = @(t,x) vodeoptions.Mass (t, x, vodeoptions.vfunarguments{:});
361 vfun = @(t,x) vmass (t, x, vodeoptions.vfunarguments{:}) ...
362 \ vfun (t, x, vodeoptions.vfunarguments{:});
363 else ## if (vmassdependence == false)
364 vmass = @(t) vodeoptions.Mass (t, vodeoptions.vfunarguments{:});
365 vfun = @(t,x) vmass (t, vodeoptions.vfunarguments{:}) ...
366 \ vfun (t, x, vodeoptions.vfunarguments{:});
367 endif
368 endif
369
370 switch integrate_func
371 case "adaptive"
372 solution = integrate_adaptive (@runge_kutta_45_dorpri, ...
373 vorder, vfun, vslot, vinit, SubOpts);
374 case "n_steps"
375 solution = integrate_n_steps (@runge_kutta_45_dorpri, ...
376 vfun, vslot(1), vinit, ...
377 vodeoptions.TimeStepSize, ...
378 vodeoptions.TimeStepNumber, SubOpts);
379 case "const"
380 solution = integrate_const (@runge_kutta_45_dorpri, ...
381 vfun, vslot, vinit, ...
382 vodeoptions.TimeStepSize, SubOpts);
383 endswitch
384
385 ## Postprocessing, do whatever when terminating integration algorithm
386 if (vodeoptions.vhaveoutputfunction) ## Cleanup plotter
387 feval (vodeoptions.OutputFcn, solution.t(end), ...
388 solution.x(end,:)', "done", vodeoptions.vfunarguments{:});
389 endif
390 if (vodeoptions.vhaveeventfunction) ## Cleanup event function handling
391 odepkg_event_handle (vodeoptions.Events, solution.t(end), ...
392 solution.x(end,:)', "done", vodeoptions.vfunarguments{:});
393 endif
394
395 ## Print additional information if option Stats is set
396 if (strcmp (vodeoptions.Stats, "on"))
397 vhavestats = true;
398 vnsteps = solution.vcntloop-2; ## vcntloop from 2..end
399 vnfailed = (solution.vcntcycles-1)-(solution.vcntloop-2)+1; ## vcntcycl from 1..end
400 vnfevals = 7*(solution.vcntcycles-1); ## number of ode evaluations
401 vndecomps = 0; ## number of LU decompositions
402 vnpds = 0; ## number of partial derivatives
403 vnlinsols = 0; ## no. of solutions of linear systems
404 ## Print cost statistics if no output argument is given
405 if (nargout == 0)
406 vmsg = fprintf (1, "Number of successful steps: %d\n", vnsteps);
407 vmsg = fprintf (1, "Number of failed attempts: %d\n", vnfailed);
408 vmsg = fprintf (1, "Number of function calls: %d\n", vnfevals);
409 endif
410 else
411 vhavestats = false;
412 endif
413
414 if (nargout == 1) ## Sort output variables, depends on nargout
415 varargout{1}.x = solution.t; ## Time stamps are saved in field x
416 varargout{1}.y = solution.x; ## Results are saved in field y
417 varargout{1}.solver = vsolver; ## Solver name is saved in field solver
418 if (vodeoptions.vhaveeventfunction)
419 varargout{1}.ie = solution.vevent{2}; ## Index info which event occured
420 varargout{1}.xe = solution.vevent{3}; ## Time info when an event occured
421 varargout{1}.ye = solution.vevent{4}; ## Results when an event occured
422 endif
423 if (vhavestats)
424 varargout{1}.stats = struct;
425 varargout{1}.stats.nsteps = vnsteps;
426 varargout{1}.stats.nfailed = vnfailed;
427 varargout{1}.stats.nfevals = vnfevals;
428 varargout{1}.stats.npds = vnpds;
429 varargout{1}.stats.ndecomps = vndecomps;
430 varargout{1}.stats.nlinsols = vnlinsols;
431 endif
432 elseif (nargout == 2)
433 varargout{1} = solution.t; ## Time stamps are first output argument
434 varargout{2} = solution.x; ## Results are second output argument
435 elseif (nargout == 5)
436 varargout{1} = solution.t; ## Same as (nargout == 2)
437 varargout{2} = solution.x; ## Same as (nargout == 2)
438 varargout{3} = []; ## LabMat doesn't accept lines like
439 varargout{4} = []; ## varargout{3} = varargout{4} = [];
440 varargout{5} = [];
441 if (vodeoptions.vhaveeventfunction)
442 varargout{3} = solution.vevent{3}; ## Time info when an event occured
443 varargout{4} = solution.vevent{4}; ## Results when an event occured
444 varargout{5} = solution.vevent{2}; ## Index info which event occured
445 endif
446 endif
447
448 endfunction
449
450 %! # We are using the "Van der Pol" implementation for all tests that
451 %! # are done for this function.
452 %! # For further tests we also define a reference solution (computed at high accuracy)
453 %!function [ydot] = fpol (vt, vy) ## The Van der Pol
454 %! ydot = [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)];
455 %!endfunction
456 %!function [vref] = fref () ## The computed reference sol
457 %! vref = [0.32331666704577, -1.83297456798624];
458 %!endfunction
459 %!function [vjac] = fjac (vt, vy, varargin) ## its Jacobian
460 %! vjac = [0, 1; -1 - 2 * vy(1) * vy(2), 1 - vy(1)^2];
461 %!function [vjac] = fjcc (vt, vy, varargin) ## sparse type
462 %! vjac = sparse ([0, 1; -1 - 2 * vy(1) * vy(2), 1 - vy(1)^2])
463 %!function [vval, vtrm, vdir] = feve (vt, vy, varargin)
464 %! vval = fpol (vt, vy, varargin); ## We use the derivatives
465 %! vtrm = zeros (2,1); ## that's why component 2
466 %! vdir = ones (2,1); ## seems to not be exact
467 %!function [vval, vtrm, vdir] = fevn (vt, vy, varargin)
468 %! vval = fpol (vt, vy, varargin); ## We use the derivatives
469 %! vtrm = ones (2,1); ## that's why component 2
470 %! vdir = ones (2,1); ## seems to not be exact
471 %!function [vmas] = fmas (vt, vy, varargin)
472 %! vmas = [1, 0; 0, 1]; ## Dummy mass matrix for tests
473 %!function [vmas] = fmsa (vt, vy, varargin)
474 %! vmas = sparse ([1, 0; 0, 1]); ## A sparse dummy matrix
475 %!function [vout] = fout (vt, vy, vflag, varargin)
476 %! if (regexp (char (vflag), 'init') == 1)
477 %! if (any (size (vt) ~= [2, 1])) error ('"fout" step "init"'); end
478 %! elseif (isempty (vflag))
479 %! if (any (size (vt) ~= [1, 1])) error ('"fout" step "calc"'); end
480 %! vout = false;
481 %! elseif (regexp (char (vflag), 'done') == 1)
482 %! if (any (size (vt) ~= [1, 1])) error ('"fout" step "done"'); end
483 %! else error ('"fout" invalid vflag');
484 %! end
485 %!
486 %! ## Turn off output of warning messages for all tests, turn them on
487 %! ## again if the last test is called
488 %!error ## ouput argument
489 %! warning ('off', 'OdePkg:InvalidArgument');
490 %! B = ode45 (1, [0 25], [3 15 1]);
491 %!error ## input argument number one
492 %! [vt, vy] = ode45 (1, [0 25], [3 15 1]);
493 %!error ## input argument number two
494 %! [vt, vy] = ode45 (@fpol, 1, [3 15 1]);
495 %!test ## two output arguments
496 %! [vt, vy] = ode45 (@fpol, [0 2], [2 0]);
497 %! assert ([vt(end), vy(end,:)], [2, fref], 1e-2);
498 %!test ## not too many steps
499 %! [vt, vy] = ode45 (@fpol, [0 2], [2 0]);
500 %! assert (size (vt) < 20);
501 %!test ## anonymous function instead of real function
502 %! fvdb = @(vt,vy) [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)];
503 %! [vt, vy] = ode45 (fvdb, [0 2], [2 0]);
504 %! assert ([vt(end), vy(end,:)], [2, fref], 1e-2);
505 %!test ## extra input arguments passed through
506 %! [vt, vy] = ode45 (@fpol, [0 2], [2 0], 12, 13, 'KL');
507 %! assert ([vt(end), vy(end,:)], [2, fref], 1e-2);
508 %!test ## empty OdePkg structure *but* extra input arguments
509 %! vopt = odeset;
510 %! [vt, vy] = ode45 (@fpol, [0 2], [2 0], vopt, 12, 13, 'KL');
511 %! assert ([vt(end), vy(end,:)], [2, fref], 1e-2);
512 %!error ## strange OdePkg structure
513 %! vopt = struct ('foo', 1);
514 %! [vt, vy] = ode45 (@fpol, [0 2], [2 0], vopt);
515 %!test ## Solve vdp in fixed step sizes
516 %! vopt = odeset('TimeStepSize', 0.1);
517 %! [vt, vy] = ode45 (@fpol, [0,2], [2 0], vopt);
518 %! assert (vt(:), [0:0.1:2]', 1e-2);
519 %!test ## Solve another anonymous function below zero
520 %! vref = [0, 14.77810590694212];
521 %! [vt, vy] = ode45 (@(t,y) y, [-2 0], 2);
522 %! assert ([vt(end), vy(end,:)], vref, 1e-1);
523 %!test ## InitialStep option
524 %! vopt = odeset ('InitialStep', 1e-8);
525 %! [vt, vy] = ode45 (@fpol, [0 0.2], [2 0], vopt);
526 %! assert ([vt(2)-vt(1)], [1e-8], 1e-9);
527 %!test ## MaxStep option
528 %! vopt = odeset ('MaxStep', 1e-3);
529 %! vsol = ode45 (@fpol, [0 0.2], [2 0], vopt);
530 %! assert ([vsol.x(5)-vsol.x(4)], [1e-3], 1e-3);
531 %!test ## Solve with intermidiate step
532 %! vsol = ode45 (@fpol, [0 1 2], [2 0]);
533 %! assert (any((vsol.x-1) == 0));
534 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
535 %!test ## Solve in backward direction starting at t=0
536 %! vref = [-1.205364552835178, 0.951542399860817];
537 %! vsol = ode45 (@fpol, [0 -2], [2 0]);
538 %! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-2);
539 %!test ## Solve in backward direction starting at t=2
540 %! vref = [-1.205364552835178, 0.951542399860817];
541 %! vsol = ode45 (@fpol, [2 -2], fref);
542 %! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-2);
543 %!test ## Solve in backward direction starting at t=2, with intermidiate step
544 %! vref = [-1.205364552835178, 0.951542399860817];
545 %! vsol = ode45 (@fpol, [2 0 -2], fref);
546 %! idx = find(vsol.x < 0, 1, 'first') - 1;
547 %! assert ([vsol.x(idx), vsol.y(idx,:)], [0 2 0], 1e-2);
548 %! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-2);
549 %!test ## Solve another anonymous function in backward direction
550 %! vref = [-1, 0.367879437558975];
551 %! vsol = ode45 (@(t,y) y, [0 -1], 1);
552 %! assert ([vsol.x(end), vsol.y(end,:)], vref, 1e-3);
553 %!test ## Solve another anonymous function below zero
554 %! vref = [0, 14.77810590694212];
555 %! vsol = ode45 (@(t,y) y, [-2 0], 2);
556 %! assert ([vsol.x(end), vsol.y(end,:)], vref, 1e-3);
557 %!test ## Solve in backward direction starting at t=0 with MaxStep option
558 %! vref = [-1.205364552835178, 0.951542399860817];
559 %! vopt = odeset ('MaxStep', 1e-3);
560 %! vsol = ode45 (@fpol, [0 -2], [2 0], vopt);
561 %! assert ([abs(vsol.x(8)-vsol.x(7))], [1e-3], 1e-3);
562 %! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-3);
563 %!test ## AbsTol option
564 %! vopt = odeset ('AbsTol', 1e-5);
565 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
566 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
567 %!test ## AbsTol and RelTol option
568 %! vopt = odeset ('AbsTol', 1e-8, 'RelTol', 1e-8);
569 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
570 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
571 %!test ## RelTol and NormControl option -- higher accuracy
572 %! vopt = odeset ('RelTol', 1e-8, 'NormControl', 'on');
573 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
574 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-5);
575 %!test ## Keeps initial values while integrating
576 %! vopt = odeset ('NonNegative', 2);
577 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
578 %! assert ([vsol.x(end), vsol.y(end,:)], [2, 2, 0], 0.5);
579 %!test ## Details of OutputSel and Refine can't be tested
580 %! vopt = odeset ('OutputFcn', @fout, 'OutputSel', 1, 'Refine', 5);
581 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
582 %!test ## Details of OutputSave can't be tested
583 %! vopt = odeset ('OutputSave', 1, 'OutputSel', 1);
584 %! vsla = ode45 (@fpol, [0 2], [2 0], vopt);
585 %! vopt = odeset ('OutputSave', 2);
586 %! vslb = ode45 (@fpol, [0 2], [2 0], vopt);
587 %! assert (length (vsla.x) + 1 >= 2 * length (vslb.x))
588 %!test ## Stats must add further elements in vsol
589 %! vopt = odeset ('Stats', 'on');
590 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
591 %! assert (isfield (vsol, 'stats'));
592 %! assert (isfield (vsol.stats, 'nsteps'));
593 %!test ## Events option add further elements in vsol
594 %! vopt = odeset ('Events', @feve);
595 %! vsol = ode45 (@fpol, [0 10], [2 0], vopt);
596 %! assert (isfield (vsol, 'ie'));
597 %! assert (vsol.ie(1), 2);
598 %! assert (isfield (vsol, 'xe'));
599 %! assert (isfield (vsol, 'ye'));
600 %!test ## Events option, now stop integration
601 %! vopt = odeset ('Events', @fevn, 'NormControl', 'on');
602 %! vsol = ode45 (@fpol, [0 10], [2 0], vopt);
603 %! assert ([vsol.ie, vsol.xe, vsol.ye], ...
604 %! [2.0, 2.496110, -0.830550, -2.677589], 6e-1);
605 %!test ## Events option, five output arguments
606 %! vopt = odeset ('Events', @fevn, 'NormControl', 'on');
607 %! [vt, vy, vxe, vye, vie] = ode45 (@fpol, [0 10], [2 0], vopt);
608 %! assert ([vie, vxe, vye], ...
609 %! [2.0, 2.496110, -0.830550, -2.677589], 6e-1);
610 %!test ## Jacobian option
611 %! vopt = odeset ('Jacobian', @fjac);
612 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
613 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
614 %!test ## Jacobian option and sparse return value
615 %! vopt = odeset ('Jacobian', @fjcc);
616 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
617 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
618 %!
619 %! ## test for JPattern option is missing
620 %! ## test for Vectorized option is missing
621 %! ## test for NewtonTol option is missing
622 %! ## test for MaxNewtonIterations option is missing
623 %!
624 %!test ## Mass option as function
625 %! vopt = odeset ('Mass', @fmas);
626 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
627 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
628 %!test ## Mass option as matrix
629 %! vopt = odeset ('Mass', eye (2,2));
630 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
631 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
632 %!test ## Mass option as sparse matrix
633 %! vopt = odeset ('Mass', sparse (eye (2,2)));
634 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
635 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
636 %!test ## Mass option as function and sparse matrix
637 %! vopt = odeset ('Mass', @fmsa);
638 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
639 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
640 %!test ## Mass option as function and MStateDependence
641 %! vopt = odeset ('Mass', @fmas, 'MStateDependence', 'strong');
642 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
643 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
644 %!test ## Set BDF option to something else than default
645 %! vopt = odeset ('BDF', 'on');
646 %! [vt, vy] = ode45 (@fpol, [0 2], [2 0], vopt);
647 %! assert ([vt(end), vy(end,:)], [2, fref], 1e-3);
648 %!
649 %! ## test for MvPattern option is missing
650 %! ## test for InitialSlope option is missing
651 %! ## test for MaxOrder option is missing
652 %!
653 %! warning ('on', 'OdePkg:InvalidArgument');
654
655 ## Local Variables: ***
656 ## mode: octave ***
657 ## End: ***