Mercurial > octave

comparison scripts/ode/ode45.m @ 31263:449ed6f427cb

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) * scripts/ode/ode23.m: Remove disabling of Refine option with struct output. Modify solution struct to output two sets of solution variables: output_t, output_x and ode_t and ode_x, and transpose struct output variables for improved Matlab compatibility. Update BISTs and perform minor code formatting. * scripts/ode/ode23s.m: Make same changes as ode23.m. * scripts/ode/ode45.m: Make same changes as ode23.m. Remove comment indicating that Refine is not implemented. * scripts/ode/private/integrate_adaptive.m: Update internal handling of variables t and x, separating them into ode_t & ode_x for internal integration and output_t & output_x for function output or calls to OutputFcn. Replace prior attempt at Refine option with new implementation. Specify time output or Refine != 0 are both interpolated from internal variables (ode_t, ode_x) for output of non-struct variables and/or for use with OutputFcn. Improve event handling when multiple Events (including at least one terminal Event) are detected in a single simulation step so that all Events up to and including the first terminal one are reported, and final data point is set to that of terminal Event. Send multiple data points in a single call to OutputFcn if they are all interpolated from a single integration step. Remove warning for termination when term signal is received from Events or OutputFcn. Return both internal variables (ode_t, ode_x) and interpolated variables (output_t, output_x) to allow calling function to correctly return either struct or separate variables. * scripts/ode/private/ode_event_handler.m: Sort multiple Events in ascending time order when they are encountered in one integration step. Remove any events after the time of a terminal Event. * scripts/ode/odeset.m: Update docstring to remove indication that Refine is not implemented * scripts/ode/odeplot.m: Update docstring to indicate that input t can be a scalar or vector. Add file test. * etc/NEWS.8.md: Add descriptions of changes under General improvements and Matlab compatibility.
author Ken Marek <marek_ka@mercer.edu>
date Wed, 05 Oct 2022 16:53:01 -0400
parents e1788b1a315f
children c332a2f2959f
comparison
equal deleted inserted replaced
31262:b8f4ec18e728 31263:449ed6f427cb
154 ## Start preprocessing, have a look which options are set in odeopts, 154 ## Start preprocessing, have a look which options are set in odeopts,
155 ## check if an invalid or unused option is set 155 ## check if an invalid or unused option is set
156 [defaults, classes, attributes] = odedefaults (numel (init), 156 [defaults, classes, attributes] = odedefaults (numel (init),
157 trange(1), trange(end)); 157 trange(1), trange(end));
158 158
159 ## FIXME: Refine is not correctly implemented yet
160 defaults = odeset (defaults, "Refine", 4); 159 defaults = odeset (defaults, "Refine", 4);
161 160
162 persistent ode45_ignore_options = ... 161 persistent ode45_ignore_options = ...
163 {"BDF", "InitialSlope", "Jacobian", "JPattern", 162 {"BDF", "InitialSlope", "Jacobian", "JPattern",
164 "MassSingular", "MaxOrder", "MvPattern", "Vectorized"}; 163 "MassSingular", "MaxOrder", "MvPattern", "Vectorized"};
194 193
195 if (isempty (odeopts.InitialStep)) 194 if (isempty (odeopts.InitialStep))
196 odeopts.InitialStep = odeopts.direction * ... 195 odeopts.InitialStep = odeopts.direction * ...
197 starting_stepsize (order, fcn, trange(1), init, 196 starting_stepsize (order, fcn, trange(1), init,
198 odeopts.AbsTol, odeopts.RelTol, 197 odeopts.AbsTol, odeopts.RelTol,
199 strcmpi (odeopts.NormControl, "on"), 198 strcmpi (odeopts.NormControl,
200 odeopts.funarguments); 199 "on"), odeopts.funarguments);
201 endif 200 endif
202 201
203 if (! isempty (odeopts.Mass)) 202 if (! isempty (odeopts.Mass))
204 if (isnumeric (odeopts.Mass)) 203 if (isnumeric (odeopts.Mass))
205 havemasshandle = false; 204 havemasshandle = false;
227 fcn = @(t,x) mass (t, odeopts.funarguments{:}) ... 226 fcn = @(t,x) mass (t, odeopts.funarguments{:}) ...
228 \ fcn (t, x, odeopts.funarguments{:}); 227 \ fcn (t, x, odeopts.funarguments{:});
229 endif 228 endif
230 endif 229 endif
231 230
232 if (nargout == 1) 231 if (numel (trange) > 2)
233 ## Single output requires auto-selected intermediate times,
234 ## which is obtained by NOT specifying specific solution times.
235 trange = [trange(1); trange(end)];
236 odeopts.Refine = []; # disable Refine when single output requested
237 elseif (numel (trange) > 2)
238 odeopts.Refine = []; # disable Refine when specific times requested 232 odeopts.Refine = []; # disable Refine when specific times requested
239 endif 233 endif
240 234
241 solution = integrate_adaptive (@runge_kutta_45_dorpri, 235 solution = integrate_adaptive (@runge_kutta_45_dorpri,
242 order, fcn, trange, init, odeopts); 236 order, fcn, trange, init, odeopts);
244 ## Postprocessing, do whatever when terminating integration algorithm 238 ## Postprocessing, do whatever when terminating integration algorithm
245 if (odeopts.haveoutputfunction) # Cleanup plotter 239 if (odeopts.haveoutputfunction) # Cleanup plotter
246 feval (odeopts.OutputFcn, [], [], "done", odeopts.funarguments{:}); 240 feval (odeopts.OutputFcn, [], [], "done", odeopts.funarguments{:});
247 endif 241 endif
248 if (! isempty (odeopts.Events)) # Cleanup event function handling 242 if (! isempty (odeopts.Events)) # Cleanup event function handling
249 ode_event_handler (odeopts.Events, solution.t(end), 243 ode_event_handler (odeopts.Events, solution.ode_t(end), ...
250 solution.x(end,:).', "done", odeopts.funarguments{:}); 244 solution.ode_x(end,:).', "done", ...
245 odeopts.funarguments{:});
251 endif 246 endif
252 247
253 ## Print additional information if option Stats is set 248 ## Print additional information if option Stats is set
254 if (strcmpi (odeopts.Stats, "on")) 249 if (strcmpi (odeopts.Stats, "on"))
255 nsteps = solution.cntloop; # cntloop from 2..end 250 nsteps = solution.cntloop; # cntloop from 2..end
263 printf ("Number of failed attempts: %d\n", nfailed); 258 printf ("Number of failed attempts: %d\n", nfailed);
264 printf ("Number of function calls: %d\n", nfevals); 259 printf ("Number of function calls: %d\n", nfevals);
265 endif 260 endif
266 261
267 if (nargout == 2) 262 if (nargout == 2)
268 varargout{1} = solution.t; # Time stamps are first output argument 263 varargout{1} = solution.output_t; # Time stamps are first output argument
269 varargout{2} = solution.x; # Results are second output argument 264 varargout{2} = solution.output_x; # Results are second output argument
270 elseif (nargout == 1) 265 elseif (nargout == 1)
271 varargout{1}.x = solution.t.'; # Time stamps saved in field x (row vector) 266 varargout{1}.x = solution.ode_t.'; #Time stamps saved in field x (row vect.)
272 varargout{1}.y = solution.x.'; # Results are saved in field y (row vector) 267 varargout{1}.y = solution.ode_x.'; #Results are saved in field y (row vect.)
273 varargout{1}.solver = solver; # Solver name is saved in field solver 268 varargout{1}.solver = solver; # Solver name is saved in field solver
274 if (! isempty (odeopts.Events)) 269 if (! isempty (odeopts.Events))
275 varargout{1}.xe = solution.event{3}; # Time info when an event occurred 270 varargout{1}.xe = solution.event{3}.'; # Time info when an event occurred
276 varargout{1}.ye = solution.event{4}; # Results when an event occurred 271 varargout{1}.ye = solution.event{4}.'; # Results when an event occurred
277 varargout{1}.ie = solution.event{2}; # Index info which event occurred 272 varargout{1}.ie = solution.event{2}.'; # Index info which event occurred
278 endif 273 endif
279 if (strcmpi (odeopts.Stats, "on")) 274 if (strcmpi (odeopts.Stats, "on"))
280 varargout{1}.stats = struct (); 275 varargout{1}.stats = struct ();
281 varargout{1}.stats.nsteps = nsteps; 276 varargout{1}.stats.nsteps = nsteps;
282 varargout{1}.stats.nfailed = nfailed; 277 varargout{1}.stats.nfailed = nfailed;
285 varargout{1}.stats.ndecomps = ndecomps; 280 varargout{1}.stats.ndecomps = ndecomps;
286 varargout{1}.stats.nlinsols = nlinsols; 281 varargout{1}.stats.nlinsols = nlinsols;
287 endif 282 endif
288 elseif (nargout > 2) 283 elseif (nargout > 2)
289 varargout = cell (1,5); 284 varargout = cell (1,5);
290 varargout{1} = solution.t; 285 varargout{1} = solution.output_t;
291 varargout{2} = solution.x; 286 varargout{2} = solution.output_x;
292 if (! isempty (odeopts.Events)) 287 if (! isempty (odeopts.Events))
293 varargout{3} = solution.event{3}; # Time info when an event occurred 288 varargout{3} = solution.event{3}; # Time info when an event occurred
294 varargout{4} = solution.event{4}; # Results when an event occurred 289 varargout{4} = solution.event{4}; # Results when an event occurred
295 varargout{5} = solution.event{2}; # Index info which event occurred 290 varargout{5} = solution.event{2}; # Index info which event occurred
296 endif 291 endif
353 %! if (strcmp (flag, "init")) 348 %! if (strcmp (flag, "init"))
354 %! if (! isequal (size (t), [2, 1])) 349 %! if (! isequal (size (t), [2, 1]))
355 %! error ('fout: step "init"'); 350 %! error ('fout: step "init"');
356 %! endif 351 %! endif
357 %! elseif (isempty (flag)) 352 %! elseif (isempty (flag))
358 %! if (! isequal (size (t), [1, 1])) 353 %! # Multiple steps can be sent in one function call
354 %! if (! isequal ( size (t), size (y)))
359 %! error ('fout: step "calc"'); 355 %! error ('fout: step "calc"');
360 %! endif 356 %! endif
361 %! elseif (strcmp (flag, "done")) 357 %! elseif (strcmp (flag, "done"))
362 %! if (! isempty (t)) 358 %! if (! isempty (t))
363 %! warning ('fout: step "done"'); 359 %! warning ('fout: step "done"');
369 %! 365 %!
370 %!test # two output arguments 366 %!test # two output arguments
371 %! [t, y] = ode45 (@fpol, [0 2], [2 0]); 367 %! [t, y] = ode45 (@fpol, [0 2], [2 0]);
372 %! assert ([t(end), y(end,:)], [2, fref], 1e-2); 368 %! assert ([t(end), y(end,:)], [2, fref], 1e-2);
373 %!test # not too many steps 369 %!test # not too many steps
374 %! [t, y] = ode45 (@fpol, [0 2], [2 0]); 370 %! [t, y] = ode45 (@fpol, [0 2], [2 0], odeset("Refine",1));
375 %! assert (size (t) < 20); 371 %! assert (size (t) < 20);
372 %!test # correct number of steps with Refine
373 %! [t1, y1] = ode45 (@fpol, [0 2], [2 0], odeset ("Refine", 1));
374 %! [t2, y2] = ode45 (@fpol, [0 2], [2 0], odeset ("Refine", 4));
375 %! [t3, y3] = ode45 (@fpol, [0 2], [2 0]); #default Refine=4
376 %! s = ode45 (@fpol, [0 2], [2 0], odeset ("Refine", 4));
377 %! assert (length (t2) == length (t3));
378 %! assert (length (t2) == 4*length (t1) - 3);
379 %! assert (length (s.x) == length (t1));
376 %!test # anonymous function instead of real function 380 %!test # anonymous function instead of real function
377 %! fvdp = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)]; 381 %! fvdp = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)];
378 %! [t, y] = ode45 (fvdp, [0 2], [2 0]); 382 %! [t, y] = ode45 (fvdp, [0 2], [2 0]);
379 %! assert ([t(end), y(end,:)], [2, fref], 1e-2); 383 %! assert ([t(end), y(end,:)], [2, fref], 1e-2);
380 %!test # string instead of function 384 %!test # string instead of function
390 %!test # Solve another anonymous function below zero 394 %!test # Solve another anonymous function below zero
391 %! vref = [0, 14.77810590694212]; 395 %! vref = [0, 14.77810590694212];
392 %! [t, y] = ode45 (@(t,y) y, [-2 0], 2); 396 %! [t, y] = ode45 (@(t,y) y, [-2 0], 2);
393 %! assert ([t(end), y(end,:)], vref, 1e-1); 397 %! assert ([t(end), y(end,:)], vref, 1e-1);
394 %!test # InitialStep option 398 %!test # InitialStep option
395 %! opt = odeset ("InitialStep", 1e-8); 399 %! opt = odeset ("InitialStep", 1e-8, "Refine", 1);
396 %! [t, y] = ode45 (@fpol, [0 0.2], [2 0], opt); 400 %! [t, y] = ode45 (@fpol, [0 0.2], [2 0], opt);
397 %! assert ([t(2)-t(1)], [1e-8], 1e-9); 401 %! assert ([t(2)-t(1)], [1e-8], 1e-9);
398 %!test # MaxStep option 402 %!test # MaxStep option
399 %! opt = odeset ("MaxStep", 1e-3); 403 %! opt = odeset ("MaxStep", 1e-3);
400 %! sol = ode45 (@fpol, [0 0.2], [2 0], opt); 404 %! sol = ode45 (@fpol, [0 0.2], [2 0], opt);
445 %! assert ([sol.x(end); sol.y(:,end)], [2; fref'], 1e-5); 449 %! assert ([sol.x(end); sol.y(:,end)], [2; fref'], 1e-5);
446 %!test # Keeps initial values while integrating 450 %!test # Keeps initial values while integrating
447 %! opt = odeset ("NonNegative", 2); 451 %! opt = odeset ("NonNegative", 2);
448 %! sol = ode45 (@fpol, [0 2], [2 0], opt); 452 %! sol = ode45 (@fpol, [0 2], [2 0], opt);
449 %! assert ([sol.x(end); sol.y(:,end)], [2; 2; 0], 0.5); 453 %! assert ([sol.x(end); sol.y(:,end)], [2; 2; 0], 0.5);
450 %!test # Details of OutputSel and Refine can't be tested 454 %!test # Details of OutputSel can't be tested
451 %! opt = odeset ("OutputFcn", @fout, "OutputSel", 1, "Refine", 5); 455 %! opt = odeset ("OutputFcn", @fout, "OutputSel", 1);
452 %! sol = ode45 (@fpol, [0 2], [2 0], opt); 456 %! sol = ode45 (@fpol, [0 2], [2 0], opt);
453 %!test # Stats must add further elements in sol 457 %!test # Stats must add further elements in sol
454 %! opt = odeset ("Stats", "on"); 458 %! opt = odeset ("Stats", "on");
455 %! stat_str = evalc ("sol = ode45 (@fpol, [0 2], [2 0], opt);"); 459 %! stat_str = evalc ("sol = ode45 (@fpol, [0 2], [2 0], opt);");
456 %! assert (strncmp (stat_str, "Number of successful steps:", 27)); 460 %! assert (strncmp (stat_str, "Number of successful steps:", 27));
462 %! assert (isfield (sol, "ie")); 466 %! assert (isfield (sol, "ie"));
463 %! assert (sol.ie(1), 2); 467 %! assert (sol.ie(1), 2);
464 %! assert (isfield (sol, "xe")); 468 %! assert (isfield (sol, "xe"));
465 %! assert (isfield (sol, "ye")); 469 %! assert (isfield (sol, "ye"));
466 %!test # Events option, now stop integration 470 %!test # Events option, now stop integration
467 %! warning ("off", "integrate_adaptive:unexpected_termination", "local");
468 %! opt = odeset ("Events", @fevn, "NormControl", "on"); 471 %! opt = odeset ("Events", @fevn, "NormControl", "on");
469 %! sol = ode45 (@fpol, [0 10], [2 0], opt); 472 %! sol = ode45 (@fpol, [0 10], [2 0], opt);
470 %! assert ([sol.ie, sol.xe, sol.ye], 473 %! assert ([sol.ie, sol.xe, sol.ye.'],
471 %! [2.0, 2.496110, -0.830550, -2.677589], 6e-1); 474 %! [2.0, 2.496110, -0.830550, -2.677589], 6e-1);
472 %!test # Events option, five output arguments 475 %!test # Events option, five output arguments
473 %! warning ("off", "integrate_adaptive:unexpected_termination", "local");
474 %! opt = odeset ("Events", @fevn, "NormControl", "on"); 476 %! opt = odeset ("Events", @fevn, "NormControl", "on");
475 %! [t, y, vxe, ye, vie] = ode45 (@fpol, [0 10], [2 0], opt); 477 %! [t, y, vxe, ye, vie] = ode45 (@fpol, [0 10], [2 0], opt);
476 %! assert ([vie, vxe, ye], 478 %! assert ([vie, vxe, ye],
477 %! [2.0, 2.496110, -0.830550, -2.677589], 6e-1); 479 %! [2.0, 2.496110, -0.830550, -2.677589], 6e-1);
478 %!test # Mass option as function 480 %!test # Mass option as function
509 ## "Vectorized" 511 ## "Vectorized"
510 512
511 %!test # Check that imaginary part of solution does not get inverted 513 %!test # Check that imaginary part of solution does not get inverted
512 %! sol = ode45 (@(x,y) 1, [0 1], 1i); 514 %! sol = ode45 (@(x,y) 1, [0 1], 1i);
513 %! assert (imag (sol.y), ones (size (sol.y))); 515 %! assert (imag (sol.y), ones (size (sol.y)));
514 %! [x, y] = ode45 (@(x,y) 1, [0 1], 1i); 516 %! [x, y] = ode45 (@(x,y) 1, [0 1], 1i, odeset ("Refine", 1));
515 %! assert (imag (y), ones (size (y))); 517 %! assert (imag (y), ones (size (y)));
518
519 ## FIXME: convert to demo or a visible=off test with failable assert/error
520 ## statemments
521 ##%!test # Make sure odeplot works (default OutputFcn when no return value)
522 ##%! ode45 (@fpol, [0 2], [2 0]);
523 ##%! close all
516 524
517 %!error <Invalid call> ode45 () 525 %!error <Invalid call> ode45 ()
518 %!error <Invalid call> ode45 (1) 526 %!error <Invalid call> ode45 (1)
519 %!error <Invalid call> ode45 (1,2) 527 %!error <Invalid call> ode45 (1,2)
520 %!error <TRANGE must be a numeric> ode45 (@fpol, {[0 25]}, [3 15 1]) 528 %!error <TRANGE must be a numeric> ode45 (@fpol, {[0 25]}, [3 15 1])