Mercurial > octave-nkf

comparison scripts/ode/ode45.m @ 20635:a22d8a2eb0e5

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
fix adaptive strategy in ode solvers. * script/ode/ode45.m: remove unused option OutputSave * script/ode/private/integrate_adaptive.m: rewrite algorithm to be more compatible. * script/ode/private/runge_kutta_45_dorpri.m: use kahan summation for time increment.
author Carlo de Falco <carlo.defalco@polimi.it>
date Sun, 11 Oct 2015 18:44:58 +0200
parents 45151de7423f
children 01586012300e
comparison
equal deleted inserted replaced
20634:1c5a86b7f838 20635:a22d8a2eb0e5
236 vodeoptions.vhaveoutputselection = true; 236 vodeoptions.vhaveoutputselection = true;
237 else 237 else
238 vodeoptions.vhaveoutputselection = false; 238 vodeoptions.vhaveoutputselection = false;
239 endif 239 endif
240 240
241 ## Implementation of the option OutputSave has been finished. 241 ## "OutputSave" option will be ignored.
242 ## This option can be set by the user to another value than default value. 242 if (! isempty (vodeoptions.OutputSave))
243 if (isempty (vodeoptions.OutputSave)) 243 warning ("OdePkg:InvalidArgument",
244 vodeoptions.OutputSave = 1; 244 "Option 'OutputSave' will be ignored.");
245 endif 245 endif
246 246
247 ## Implementation of the option Refine has been finished. 247 ## Implementation of the option Refine has been finished.
248 ## This option can be set by the user to another value than default value. 248 ## This option can be set by the user to another value than default value.
249 if (vodeoptions.Refine > 0) 249 if (vodeoptions.Refine > 0)
404 endif 404 endif
405 405
406 ## Print additional information if option Stats is set 406 ## Print additional information if option Stats is set
407 if (strcmp (vodeoptions.Stats, "on")) 407 if (strcmp (vodeoptions.Stats, "on"))
408 vhavestats = true; 408 vhavestats = true;
409 vnsteps = solution.vcntloop-2; # vcntloop from 2..end 409 vnsteps = solution.vcntloop-2; # vcntloop from 2..end
410 vnfailed = (solution.vcntcycles-1)-(vnsteps)+1; # vcntcycl from 1..end 410 vnfailed = (solution.vcntcycles-1)-(vnsteps)+1; # vcntcycl from 1..end
411 vnfevals = 7*(solution.vcntcycles-1); # number of ode evaluations 411 vnfevals = 7*(solution.vcntcycles-1); # number of ode evaluations
412 vndecomps = 0; # number of LU decompositions 412 vndecomps = 0; # number of LU decompositions
413 vnpds = 0; # number of partial derivatives 413 vnpds = 0; # number of partial derivatives
414 vnlinsols = 0; # no. of linear systems solutions 414 vnlinsols = 0; # no. of linear systems solutions
415 ## Print cost statistics if no output argument is given 415 ## Print cost statistics if no output argument is given
416 if (nargout == 0) 416 if (nargout == 0)
417 printf ("Number of successful steps: %d\n", vnsteps); 417 printf ("Number of successful steps: %d\n", vnsteps);
418 printf ("Number of failed attempts: %d\n", vnfailed); 418 printf ("Number of failed attempts: %d\n", vnfailed);
419 printf ("Number of function calls: %d\n", vnfevals); 419 printf ("Number of function calls: %d\n", vnfevals);
596 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt); 596 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
597 %! assert ([vsol.x(end), vsol.y(end,:)], [2, 2, 0], 0.5); 597 %! assert ([vsol.x(end), vsol.y(end,:)], [2, 2, 0], 0.5);
598 %!test # Details of OutputSel and Refine can't be tested 598 %!test # Details of OutputSel and Refine can't be tested
599 %! vopt = odeset ("OutputFcn", @fout, "OutputSel", 1, "Refine", 5); 599 %! vopt = odeset ("OutputFcn", @fout, "OutputSel", 1, "Refine", 5);
600 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt); 600 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
601 %!test # Details of OutputSave can't be tested
602 %! vopt = odeset ("OutputSave", 1, "OutputSel", 1);
603 %! vsla = ode45 (@fpol, [0 2], [2 0], vopt);
604 %! vopt = odeset ("OutputSave", 2);
605 %! vslb = ode45 (@fpol, [0 2], [2 0], vopt);
606 %! assert (length (vsla.x) + 1 >= 2 * length (vslb.x))
607 %!test # Stats must add further elements in vsol 601 %!test # Stats must add further elements in vsol
608 %! vopt = odeset ("Stats", "on"); 602 %! vopt = odeset ("Stats", "on");
609 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt); 603 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
610 %! assert (isfield (vsol, "stats")); 604 %! assert (isfield (vsol, "stats"));
611 %! assert (isfield (vsol.stats, "nsteps")); 605 %! assert (isfield (vsol.stats, "nsteps"));