octave: scripts/ode/private/integrate

annotate scripts/ode/private/integrate_adaptive.m @ 31414:5341ff8d21d0

Backed out changeset 88fff8521d76 These changes were committed under an incorrect ID.

author	John W. Eaton <jwe@octave.org>
date	Fri, 11 Nov 2022 08:31:11 -0500
parents	88fff8521d76
children	33d5c1c41bbc

rev	line source
27923 bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	1 ########################################################################
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	2 ##
30564 796f54d4ddbf update Octave Project Developers copyright for the new year John W. Eaton <jwe@octave.org> parents: 29358 diff changeset	3 ## Copyright (C) 2013-2022 The Octave Project Developers
27918 b442ec6dda5c use centralized file for copyright info for individual contributors John W. Eaton <jwe@octave.org> parents: 26376 diff changeset	4 ##
27923 bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	5 ## See the file COPYRIGHT.md in the top-level directory of this
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	6 ## distribution or <https://octave.org/copyright/>.
20536 6256f6e366ac Fix copyright text in private ode functions Carlo de Falco <carlo.defalco@polimi.it> parents: 20533 diff changeset	7 ##
6256f6e366ac Fix copyright text in private ode functions Carlo de Falco <carlo.defalco@polimi.it> parents: 20533 diff changeset	8 ## This file is part of Octave.
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	9 ##
24534 194eb4bd202b maint: Update punctuation for GPL v3 license text. Rik <rik@octave.org> parents: 23220 diff changeset	10 ## Octave is free software: you can redistribute it and/or modify it
20536 6256f6e366ac Fix copyright text in private ode functions Carlo de Falco <carlo.defalco@polimi.it> parents: 20533 diff changeset	11 ## under the terms of the GNU General Public License as published by
24534 194eb4bd202b maint: Update punctuation for GPL v3 license text. Rik <rik@octave.org> parents: 23220 diff changeset	12 ## the Free Software Foundation, either version 3 of the License, or
22755 3a2b891d0b33 maint: Standardize Copyright formatting. Rik <rik@octave.org> parents: 22669 diff changeset	13 ## (at your option) any later version.
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	14 ##
20536 6256f6e366ac Fix copyright text in private ode functions Carlo de Falco <carlo.defalco@polimi.it> parents: 20533 diff changeset	15 ## Octave is distributed in the hope that it will be useful, but
6256f6e366ac Fix copyright text in private ode functions Carlo de Falco <carlo.defalco@polimi.it> parents: 20533 diff changeset	16 ## WITHOUT ANY WARRANTY; without even the implied warranty of
22755 3a2b891d0b33 maint: Standardize Copyright formatting. Rik <rik@octave.org> parents: 22669 diff changeset	17 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3a2b891d0b33 maint: Standardize Copyright formatting. Rik <rik@octave.org> parents: 22669 diff changeset	18 ## GNU General Public License for more details.
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	19 ##
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	20 ## You should have received a copy of the GNU General Public License
20536 6256f6e366ac Fix copyright text in private ode functions Carlo de Falco <carlo.defalco@polimi.it> parents: 20533 diff changeset	21 ## along with Octave; see the file COPYING. If not, see
24534 194eb4bd202b maint: Update punctuation for GPL v3 license text. Rik <rik@octave.org> parents: 23220 diff changeset	22 ## <https://www.gnu.org/licenses/>.
27923 bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	23 ##
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	24 ########################################################################
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	25
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	26 ## -- texinfo --
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	27 ## @deftypefn {} {@var{solution} =} integrate_adaptive (@var{@@stepper}, @var{order}, @var{@@fcn}, @var{tspan}, @var{x0}, @var{options})
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	28 ##
20548 25623ef2ff4f doc: Rewrite docstrings for ode* family of functions. Rik <rik@octave.org> parents: 20543 diff changeset	29 ## This function file can be called by an ODE solver function in order to
22626 869c02fde46c Further clean-up of ode functions. Rik <rik@octave.org> parents: 22605 diff changeset	30 ## integrate the set of ODEs on the interval @var{[t0, t1]} with an adaptive
869c02fde46c Further clean-up of ode functions. Rik <rik@octave.org> parents: 22605 diff changeset	31 ## timestep.
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	32 ##
22626 869c02fde46c Further clean-up of ode functions. Rik <rik@octave.org> parents: 22605 diff changeset	33 ## The function returns a structure @var{solution} with two fields: @var{t}
20621 b92f8e148936 maint: Continued clean-up of functions in ode/private dir. Rik <rik@octave.org> parents: 20604 diff changeset	34 ## and @var{y}. @var{t} is a column vector and contains the time stamps.
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	35 ## @var{y} is a matrix in which each column refers to a different unknown
20621 b92f8e148936 maint: Continued clean-up of functions in ode/private dir. Rik <rik@octave.org> parents: 20604 diff changeset	36 ## of the problem and the row number is the same as the @var{t} row number.
b92f8e148936 maint: Continued clean-up of functions in ode/private dir. Rik <rik@octave.org> parents: 20604 diff changeset	37 ## Thus, each row of the matrix @var{y} contains the values of all unknowns at
b92f8e148936 maint: Continued clean-up of functions in ode/private dir. Rik <rik@octave.org> parents: 20604 diff changeset	38 ## the time value contained in the corresponding row in @var{t}.
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	39 ##
20621 b92f8e148936 maint: Continued clean-up of functions in ode/private dir. Rik <rik@octave.org> parents: 20604 diff changeset	40 ## The first input argument must be a function handle or inline function
b92f8e148936 maint: Continued clean-up of functions in ode/private dir. Rik <rik@octave.org> parents: 20604 diff changeset	41 ## representing the stepper, i.e., the function responsible for step-by-step
20548 25623ef2ff4f doc: Rewrite docstrings for ode* family of functions. Rik <rik@octave.org> parents: 20543 diff changeset	42 ## integration. This function discriminates one method from the others.
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	43 ##
20548 25623ef2ff4f doc: Rewrite docstrings for ode* family of functions. Rik <rik@octave.org> parents: 20543 diff changeset	44 ## The second input argument is the order of the stepper. It is needed
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	45 ## to compute the adaptive timesteps.
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	46 ##
20621 b92f8e148936 maint: Continued clean-up of functions in ode/private dir. Rik <rik@octave.org> parents: 20604 diff changeset	47 ## The third input argument is a function handle or inline function that
b92f8e148936 maint: Continued clean-up of functions in ode/private dir. Rik <rik@octave.org> parents: 20604 diff changeset	48 ## defines the ODE:
b92f8e148936 maint: Continued clean-up of functions in ode/private dir. Rik <rik@octave.org> parents: 20604 diff changeset	49 ##
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	50 ## @ifhtml
22299 9fc91bb2aec3 doc: grammarcheck documentation for 4.2 release. Rik <rik@octave.org> parents: 22123 diff changeset	51 ##
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	52 ## @example
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	53 ## @math{y' = f(t,y)}
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	54 ## @end example
22299 9fc91bb2aec3 doc: grammarcheck documentation for 4.2 release. Rik <rik@octave.org> parents: 22123 diff changeset	55 ##
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	56 ## @end ifhtml
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	57 ## @ifnothtml
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	58 ## @math{y' = f(t,y)}.
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	59 ## @end ifnothtml
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	60 ##
20621 b92f8e148936 maint: Continued clean-up of functions in ode/private dir. Rik <rik@octave.org> parents: 20604 diff changeset	61 ## The fourth input argument is the time vector which defines the integration
20716 1ecee53513d7 doc: Peridodic grammar check of documentation. Rik <rik@octave.org> parents: 20636 diff changeset	62 ## interval, i.e., @var{[tspan(1), tspan(end)]} and all intermediate elements
20621 b92f8e148936 maint: Continued clean-up of functions in ode/private dir. Rik <rik@octave.org> parents: 20604 diff changeset	63 ## are taken as times at which the solution is required.
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	64 ##
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	65 ## The fifth argument represents the initial conditions for the ODEs and the
21546 f7f97d7e9294 doc: Wrap m-file docstrings to 79 characters + newline (80 total). Rik <rik@octave.org> parents: 20852 diff changeset	66 ## last input argument contains some options that may be needed for the
f7f97d7e9294 doc: Wrap m-file docstrings to 79 characters + newline (80 total). Rik <rik@octave.org> parents: 20852 diff changeset	67 ## stepper.
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	68 ##
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	69 ## @end deftypefn
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	70
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	71 function solution = integrate_adaptive (stepper, order, fcn, tspan, x0,
20621 b92f8e148936 maint: Continued clean-up of functions in ode/private dir. Rik <rik@octave.org> parents: 20604 diff changeset	72 options)
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	73
20600 a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	74 fixed_times = numel (tspan) > 2;
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	75
31263 449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	76 t_new = t_old = ode_t = output_t = tspan(1);
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	77 x_new = x_old = ode_x = output_x = x0(:);
20602 756b052037fb avoid stepping beyond end of thspan in ode solvers Carlo de Falco <carlo.defalco@polimi.it> parents: 20601 diff changeset	78
20621 b92f8e148936 maint: Continued clean-up of functions in ode/private dir. Rik <rik@octave.org> parents: 20604 diff changeset	79 ## Get first initial timestep
22605 177e0c71bcc0 make sure the additional function arguments are always passed to odefun. Carlo de Falco <carlo.defalco@polimi.it> parents: 22596 diff changeset	80 dt = options.InitialStep;
20621 b92f8e148936 maint: Continued clean-up of functions in ode/private dir. Rik <rik@octave.org> parents: 20604 diff changeset	81 if (isempty (dt))
31263 449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	82 dt = starting_stepsize (order, fcn, ode_t, ode_x,
22626 869c02fde46c Further clean-up of ode functions. Rik <rik@octave.org> parents: 22605 diff changeset	83 options.AbsTol, options.RelTol,
869c02fde46c Further clean-up of ode functions. Rik <rik@octave.org> parents: 22605 diff changeset	84 strcmp (options.NormControl, "on"),
869c02fde46c Further clean-up of ode functions. Rik <rik@octave.org> parents: 22605 diff changeset	85 options.funarguments);
20621 b92f8e148936 maint: Continued clean-up of functions in ode/private dir. Rik <rik@octave.org> parents: 20604 diff changeset	86 endif
20602 756b052037fb avoid stepping beyond end of thspan in ode solvers Carlo de Falco <carlo.defalco@polimi.it> parents: 20601 diff changeset	87
22605 177e0c71bcc0 make sure the additional function arguments are always passed to odefun. Carlo de Falco <carlo.defalco@polimi.it> parents: 22596 diff changeset	88 dir = options.direction;
20600 a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	89 dt = dir * min (abs (dt), options.MaxStep);
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	90
20600 a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	91 options.comp = 0.0;
20602 756b052037fb avoid stepping beyond end of thspan in ode solvers Carlo de Falco <carlo.defalco@polimi.it> parents: 20601 diff changeset	92
20552 eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir. Rik <rik@octave.org> parents: 20548 diff changeset	93 ## Factor multiplying the stepsize guess
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	94 facmin = 0.8;
20600 a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	95 facmax = 1.5;
20552 eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir. Rik <rik@octave.org> parents: 20548 diff changeset	96 fac = 0.38^(1/(order+1)); # formula taken from Hairer
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	97
31263 449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	98 ## Initialize Refine value
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	99 refine = options.Refine;
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	100 if isempty (refine)
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	101 refine = 1;
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	102 elseif ((refine != round (refine)) \|\| (refine < 1))
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	103 refine = 1;
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	104 warning ("integrate_adaptive:invalid_refine",
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	105 ["Invalid value of Refine. Refine must be a positive " ...
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	106 "integer. Setting Refine = 1."] );
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	107 endif
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	108
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	109 ## Initialize the OutputFcn
20634 80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files. Rik <rik@octave.org> parents: 20621 diff changeset	110 if (options.haveoutputfunction)
22593 dba5074bdc79 simplify options management in ode solvers Carlo de Falco <carlo.defalco@polimi.it> parents: 22538 diff changeset	111 if (! isempty (options.OutputSel))
31263 449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	112 solution.retout = output_x(options.OutputSel, end);
20602 756b052037fb avoid stepping beyond end of thspan in ode solvers Carlo de Falco <carlo.defalco@polimi.it> parents: 20601 diff changeset	113 else
31263 449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	114 solution.retout = output_x;
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	115 endif
22639 7efa2d0e22c9 More Matlab-compatible implementation of OutputFcn. Clean up odeplot.m. Rik <rik@octave.org> parents: 22626 diff changeset	116 feval (options.OutputFcn, tspan, solution.retout, "init",
7efa2d0e22c9 More Matlab-compatible implementation of OutputFcn. Clean up odeplot.m. Rik <rik@octave.org> parents: 22626 diff changeset	117 options.funarguments{:});
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	118 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	119
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	120 ## Initialize the EventFcn
23079 17dc6c7ef427 Ignore ODE event function at t==0 for Matlab compatibility (bug #49919). Rik <rik@octave.org> parents: 22916 diff changeset	121 have_EventFcn = false;
22593 dba5074bdc79 simplify options management in ode solvers Carlo de Falco <carlo.defalco@polimi.it> parents: 22538 diff changeset	122 if (! isempty (options.Events))
23079 17dc6c7ef427 Ignore ODE event function at t==0 for Matlab compatibility (bug #49919). Rik <rik@octave.org> parents: 22916 diff changeset	123 have_EventFcn = true;
31414 5341ff8d21d0 Backed out changeset 88fff8521d76 John W. Eaton <jwe@octave.org> parents: 31409 diff changeset	124 ode_event_handler (options.Events, tspan(1), ode_x,
5341ff8d21d0 Backed out changeset 88fff8521d76 John W. Eaton <jwe@octave.org> parents: 31409 diff changeset	125 "init", options.funarguments{:});
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	126 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	127
20634 80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files. Rik <rik@octave.org> parents: 20621 diff changeset	128 if (options.havenonnegative)
20600 a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	129 nn = options.NonNegative;
a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	130 endif
20602 756b052037fb avoid stepping beyond end of thspan in ode solvers Carlo de Falco <carlo.defalco@polimi.it> parents: 20601 diff changeset	131
21955 0f3e875d9078 Fix statistics on solution solving for ode45, ode23 (bug #48243). jcorno <jacopo.corno@gmail.com> parents: 21546 diff changeset	132 solution.cntloop = 0;
0f3e875d9078 Fix statistics on solution solving for ode45, ode23 (bug #48243). jcorno <jacopo.corno@gmail.com> parents: 21546 diff changeset	133 solution.cntcycles = 0;
20634 80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files. Rik <rik@octave.org> parents: 20621 diff changeset	134 solution.cntsave = 2;
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files. Rik <rik@octave.org> parents: 20621 diff changeset	135 solution.unhandledtermination = true;
20600 a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	136 ireject = 0;
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	137
22645 34bf558de127 integrate_adaptive.m: Take strcmp call out of for loop for performance. Rik <rik@octave.org> parents: 22639 diff changeset	138 NormControl = strcmp (options.NormControl, "on");
20602 756b052037fb avoid stepping beyond end of thspan in ode solvers Carlo de Falco <carlo.defalco@polimi.it> parents: 20601 diff changeset	139 k_vals = [];
20600 a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	140 iout = istep = 1;
22626 869c02fde46c Further clean-up of ode functions. Rik <rik@octave.org> parents: 22605 diff changeset	141
20600 a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	142 while (dir * t_old < dir * tspan(end))
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	143
20600 a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	144 ## Compute integration step from t_old to t_new = t_old + dt
a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	145 [t_new, options.comp] = kahan (t_old, options.comp, dt);
a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	146 [t_new, x_new, x_est, new_k_vals] = ...
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	147 stepper (fcn, t_old, x_old, dt, options, k_vals, t_new);
20600 a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	148
20636 6e81f4b37e13 Performance improvements for ODE functions. Rik <rik@octave.org> parents: 20634 diff changeset	149 solution.cntcycles += 1;
20600 a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	150
20634 80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files. Rik <rik@octave.org> parents: 20621 diff changeset	151 if (options.havenonnegative)
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files. Rik <rik@octave.org> parents: 20621 diff changeset	152 x_new(nn, end) = abs (x_new(nn, end));
20604 a260a6acb70f fix test failures introduced by a22d8a2eb0e5 Carlo de Falco <carlo.defalco@polimi.it> parents: 20602 diff changeset	153 x_est(nn, end) = abs (x_est(nn, end));
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	154 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	155
22626 869c02fde46c Further clean-up of ode functions. Rik <rik@octave.org> parents: 22605 diff changeset	156 err = AbsRel_norm (x_new, x_old, options.AbsTol, options.RelTol,
22645 34bf558de127 integrate_adaptive.m: Take strcmp call out of for loop for performance. Rik <rik@octave.org> parents: 22639 diff changeset	157 NormControl, x_est);
20602 756b052037fb avoid stepping beyond end of thspan in ode solvers Carlo de Falco <carlo.defalco@polimi.it> parents: 20601 diff changeset	158
20621 b92f8e148936 maint: Continued clean-up of functions in ode/private dir. Rik <rik@octave.org> parents: 20604 diff changeset	159 ## Accept solution only if err <= 1.0
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	160 if (err <= 1)
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	161
20636 6e81f4b37e13 Performance improvements for ODE functions. Rik <rik@octave.org> parents: 20634 diff changeset	162 solution.cntloop += 1;
6e81f4b37e13 Performance improvements for ODE functions. Rik <rik@octave.org> parents: 20634 diff changeset	163 ireject = 0; # Clear reject counter
31263 449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	164 terminal_event = false;
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	165 terminal_output = false;
20602 756b052037fb avoid stepping beyond end of thspan in ode solvers Carlo de Falco <carlo.defalco@polimi.it> parents: 20601 diff changeset	166
31263 449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	167 istep++;
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	168 ode_t(istep) = t_new;
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	169 ode_x(:, istep) = x_new;
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	170 iadd = 0; # Number of output points added this iteration
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	171
31263 449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	172 ## Check for Events
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	173 if (have_EventFcn)
31414 5341ff8d21d0 Backed out changeset 88fff8521d76 John W. Eaton <jwe@octave.org> parents: 31409 diff changeset	174 solution.event = ode_event_handler (options.Events, t_new, x_new, ...
5341ff8d21d0 Backed out changeset 88fff8521d76 John W. Eaton <jwe@octave.org> parents: 31409 diff changeset	175 [], options.funarguments{:});
31263 449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	176 ## Check for terminal Event
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	177 if (! isempty (solution.event{1}) && solution.event{1} == 1)
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	178 ode_t(istep) = solution.event{3}(end);
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	179 ode_x(:, istep) = solution.event{4}(end, :).';
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	180 solution.unhandledtermination = false;
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	181 terminal_event = true;
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	182 endif
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	183 endif
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	184
31263 449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	185 ## Interpolate to specified or Refined points
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	186 if (fixed_times)
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	187 t_caught = find ((dir * tspan(iout:end) > dir * t_old) ...
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	188 & (dir * tspan(iout:end) <= dir * ode_t(istep)));
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	189 t_caught = t_caught + iout - 1;
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	190 iadd = length (t_caught);
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	191 if (! isempty (t_caught))
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	192 output_t(t_caught) = tspan(t_caught);
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	193 iout = max (t_caught);
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	194 output_x(:, t_caught) = ...
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	195 runge_kutta_interpolate (order, [t_old t_new], [x_old x_new], ...
31414 5341ff8d21d0 Backed out changeset 88fff8521d76 John W. Eaton <jwe@octave.org> parents: 31409 diff changeset	196 output_t(t_caught), new_k_vals, dt, ...
5341ff8d21d0 Backed out changeset 88fff8521d76 John W. Eaton <jwe@octave.org> parents: 31409 diff changeset	197 fcn, options.funarguments);
20600 a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	198 endif
31263 449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	199 ## Add a possible additional output value if we found a terminal Event
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	200 if ((terminal_event == true) && ...
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	201 (dir * ode_t(istep) > dir * output_t(iout)))
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	202 iadd += 1;
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	203 iout += 1;
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	204 output_x(:, iout) = ode_x(:, istep);
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	205 output_t(iout) = ode_t(istep);
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	206 endif
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	207 elseif (refine > 1)
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	208 iadd = refine;
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	209 tadd = linspace (t_old, ode_t(istep), refine + 1);
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	210 tadd = tadd(2:end);
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	211 output_x(:, iout + (1:iadd)) = ...
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	212 runge_kutta_interpolate (order, [t_old t_new], [x_old x_new], ...
31414 5341ff8d21d0 Backed out changeset 88fff8521d76 John W. Eaton <jwe@octave.org> parents: 31409 diff changeset	213 tadd, new_k_vals, dt, fcn, ...
5341ff8d21d0 Backed out changeset 88fff8521d76 John W. Eaton <jwe@octave.org> parents: 31409 diff changeset	214 options.funarguments);
31263 449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	215 output_t(iout + (1:iadd)) = tadd;
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	216 iout = length (output_t);
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	217 else # refine = 1
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	218 iadd = 1;
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	219 iout += iadd;
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	220 output_x(:, iout) = ode_x(:, istep);
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	221 output_t(iout) = ode_t(istep);
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	222 end
20600 a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	223
31263 449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	224 ## Call OutputFcn
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	225 if ((options.haveoutputfunction) && (iadd > 0))
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	226 xadd = output_x(:, (iout-iadd+1):end);
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	227 tadd = output_t((iout-iadd+1):end);
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	228 if (! isempty (options.OutputSel))
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	229 xadd = xadd(options.OutputSel, :);
20600 a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	230 endif
31263 449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	231 stop_solve = feval (options.OutputFcn, tadd, xadd, ...
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	232 [], options.funarguments{:});
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	233 if (stop_solve)
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	234 solution.unhandledtermination = false;
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	235 terminal_output = true;
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	236 endif
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	237 endif
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	238 if (terminal_event \|\| terminal_output)
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	239 break; # break from main loop
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	240 endif
20600 a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	241
a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	242
a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	243 ## move to next time-step
a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	244 t_old = t_new;
a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	245 x_old = x_new;
a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	246 k_vals = new_k_vals;
a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	247
22626 869c02fde46c Further clean-up of ode functions. Rik <rik@octave.org> parents: 22605 diff changeset	248 else # error condition
20600 a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	249
20636 6e81f4b37e13 Performance improvements for ODE functions. Rik <rik@octave.org> parents: 20634 diff changeset	250 ireject += 1;
20600 a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	251
22626 869c02fde46c Further clean-up of ode functions. Rik <rik@octave.org> parents: 22605 diff changeset	252 ## Stop solving if, in the last 5,000 steps, no successful valid
869c02fde46c Further clean-up of ode functions. Rik <rik@octave.org> parents: 22605 diff changeset	253 ## value has been found.
20621 b92f8e148936 maint: Continued clean-up of functions in ode/private dir. Rik <rik@octave.org> parents: 20604 diff changeset	254 if (ireject >= 5_000)
b92f8e148936 maint: Continued clean-up of functions in ode/private dir. Rik <rik@octave.org> parents: 20604 diff changeset	255 error (["integrate_adaptive: Solving was not successful. ", ...
b92f8e148936 maint: Continued clean-up of functions in ode/private dir. Rik <rik@octave.org> parents: 20604 diff changeset	256 " The iterative integration loop exited at time", ...
b92f8e148936 maint: Continued clean-up of functions in ode/private dir. Rik <rik@octave.org> parents: 20604 diff changeset	257 " t = %f before the endpoint at tend = %f was reached. ", ...
b92f8e148936 maint: Continued clean-up of functions in ode/private dir. Rik <rik@octave.org> parents: 20604 diff changeset	258 " This happened because the iterative integration loop", ...
b92f8e148936 maint: Continued clean-up of functions in ode/private dir. Rik <rik@octave.org> parents: 20604 diff changeset	259 " did not find a valid solution at this time stamp. ", ...
b92f8e148936 maint: Continued clean-up of functions in ode/private dir. Rik <rik@octave.org> parents: 20604 diff changeset	260 " Try to reduce the value of 'InitialStep' and/or", ...
20600 a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	261 " 'MaxStep' with the command 'odeset'.\n"],
a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	262 t_old, tspan(end));
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	263 endif
20600 a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	264
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	265 endif
20602 756b052037fb avoid stepping beyond end of thspan in ode solvers Carlo de Falco <carlo.defalco@polimi.it> parents: 20601 diff changeset	266
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	267 ## Compute next timestep, formula taken from Hairer
20621 b92f8e148936 maint: Continued clean-up of functions in ode/private dir. Rik <rik@octave.org> parents: 20604 diff changeset	268 err += eps; # avoid divisions by zero
20600 a22d8a2eb0e5 fix adaptive strategy in ode solvers. Carlo de Falco <carlo.defalco@polimi.it> parents: 20565 diff changeset	269 dt = min (facmax, max (facmin, fac (1 / err)^(1 / (order + 1))));
20602 756b052037fb avoid stepping beyond end of thspan in ode solvers Carlo de Falco <carlo.defalco@polimi.it> parents: 20601 diff changeset	270 dt = dir * min (abs (dt), options.MaxStep);
31263 449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	271 if (! (abs (dt) > eps (ode_t(end))))
22626 869c02fde46c Further clean-up of ode functions. Rik <rik@octave.org> parents: 22605 diff changeset	272 break;
22538 21c89e691804 Stop adaptive ODE integration if timestep gets too small. Carlo de Falco <carlo.defalco@polimi.it> parents: 22323 diff changeset	273 endif
22626 869c02fde46c Further clean-up of ode functions. Rik <rik@octave.org> parents: 22605 diff changeset	274
869c02fde46c Further clean-up of ode functions. Rik <rik@octave.org> parents: 22605 diff changeset	275 ## Make sure we don't go past tpan(end)
20602 756b052037fb avoid stepping beyond end of thspan in ode solvers Carlo de Falco <carlo.defalco@polimi.it> parents: 20601 diff changeset	276 dt = dir * min (abs (dt), abs (tspan(end) - t_old));
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	277
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	278 endwhile
20602 756b052037fb avoid stepping beyond end of thspan in ode solvers Carlo de Falco <carlo.defalco@polimi.it> parents: 20601 diff changeset	279
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	280 ## Check if integration of the ode has been successful
31263 449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	281 if (dir * ode_t(end) < dir * tspan(end))
20634 80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files. Rik <rik@octave.org> parents: 20621 diff changeset	282 if (solution.unhandledtermination == true)
22538 21c89e691804 Stop adaptive ODE integration if timestep gets too small. Carlo de Falco <carlo.defalco@polimi.it> parents: 22323 diff changeset	283 warning ("integrate_adaptive:unexpected_termination",
21c89e691804 Stop adaptive ODE integration if timestep gets too small. Carlo de Falco <carlo.defalco@polimi.it> parents: 22323 diff changeset	284 [" Solving was not successful. ", ...
21c89e691804 Stop adaptive ODE integration if timestep gets too small. Carlo de Falco <carlo.defalco@polimi.it> parents: 22323 diff changeset	285 " The iterative integration loop exited at time", ...
21c89e691804 Stop adaptive ODE integration if timestep gets too small. Carlo de Falco <carlo.defalco@polimi.it> parents: 22323 diff changeset	286 " t = %f before the endpoint at tend = %f was reached. ", ...
21c89e691804 Stop adaptive ODE integration if timestep gets too small. Carlo de Falco <carlo.defalco@polimi.it> parents: 22323 diff changeset	287 " This may happen if the stepsize becomes too small. ", ...
21c89e691804 Stop adaptive ODE integration if timestep gets too small. Carlo de Falco <carlo.defalco@polimi.it> parents: 22323 diff changeset	288 " Try to reduce the value of 'InitialStep'", ...
21c89e691804 Stop adaptive ODE integration if timestep gets too small. Carlo de Falco <carlo.defalco@polimi.it> parents: 22323 diff changeset	289 " and/or 'MaxStep' with the command 'odeset'."],
31263 449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	290 ode_t(end), tspan(end));
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	291 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	292 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	293
20601 43822bda4f65 fix indexing bug introduced with a22d8a2eb0e5 Carlo de Falco <carlo.defalco@polimi.it> parents: 20600 diff changeset	294 ## Set up return structure
31263 449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	295 solution.ode_t = ode_t(:);
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	296 solution.ode_x = ode_x.';
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	297 solution.output_t = output_t(:);
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063) Ken Marek <marek_ka@mercer.edu> parents: 30893 diff changeset	298 solution.output_x = output_x.';
20602 756b052037fb avoid stepping beyond end of thspan in ode solvers Carlo de Falco <carlo.defalco@polimi.it> parents: 20601 diff changeset	299
20533 fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core. jcorno <jacopo.corno@gmail.com> parents: diff changeset	300 endfunction

Mercurial > octave

annotate scripts/ode/private/integrate_adaptive.m @ 31414:5341ff8d21d0