Mercurial > octave

annotate scripts/ode/private/integrate_adaptive.m @ 31548:c8ad083a5802 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Clean up m-files before Octave 8.1 release. * external.txi, oop.txi, Table.h, documentation.cc, gui-preferences-ed.h, lo-specfun.cc, range.tst : Eliminate triple newlines. * Map.m, MemoizedFunction.m, delaunayn.m, inputParser.m, __publish_latex_output__.m, publish.m, unpack.m, fminbnd.m, __add_default_menu__.m, gammainc.m, gallery.m, hadamard.m, weboptions.m: Add newline after keyword "function" or before keyword "endfunction" for readability. * getaudiodata.m, pkg.m : Add semicolon to end of line for error() statement. * movegui.m: Combine mutliple calls to set() into one for performance. * __unimplemented__.m (missing_functions): Remove missing functions that have been implemented. * __vectorize__.m, check_default_input.m, betaincinv.m, gammaincinv.m: Remove semicolon at end of line with "function" declaration. * weboptions.m: Remove semicolon at end of line with "if" keyword. * integrate_adaptive.m, factor.m: Use keyword "endif" rather than bare "end".
author Rik <rik@octave.org>
date Fri, 25 Nov 2022 21:23:54 -0800
parents 7286327ec4b6
children 597f3ee61a48
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
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;
31417
7286327ec4b6 Improve "Events" location in ode45, ode23, ode23s (bug #63162)
Ken Marek <marek_ka@mercer.edu>
parents: 31416
diff changeset
124 options.Events = @(t,y) options.Events (t, y, options.funarguments{:});
7286327ec4b6 Improve "Events" location in ode45, ode23, ode23s (bug #63162)
Ken Marek <marek_ka@mercer.edu>
parents: 31416
diff changeset
125 ode_event_handler (options.Events, tspan(1), ode_x, ...
7286327ec4b6 Improve "Events" location in ode45, ode23, ode23s (bug #63162)
Ken Marek <marek_ka@mercer.edu>
parents: 31416
diff changeset
126 [], order, "init");
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
127 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
128
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
129 if (options.havenonnegative)
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
130 nn = options.NonNegative;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
131 endif
20602
756b052037fb avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20601
diff changeset
132
21955
0f3e875d9078 Fix statistics on solution solving for ode45, ode23 (bug #48243).
jcorno <jacopo.corno@gmail.com>
parents: 21546
diff changeset
133 solution.cntloop = 0;
0f3e875d9078 Fix statistics on solution solving for ode45, ode23 (bug #48243).
jcorno <jacopo.corno@gmail.com>
parents: 21546
diff changeset
134 solution.cntcycles = 0;
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
135 solution.cntsave = 2;
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
136 solution.unhandledtermination = true;
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
137 ireject = 0;
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
138
22645
34bf558de127 integrate_adaptive.m: Take strcmp call out of for loop for performance.
Rik <rik@octave.org>
parents: 22639
diff changeset
139 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
140 k_vals = [];
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
141 iout = istep = 1;
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
142
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
143 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
144
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
145 ## 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
146 [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
147 [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
148 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
149
20636
6e81f4b37e13 Performance improvements for ODE functions.
Rik <rik@octave.org>
parents: 20634
diff changeset
150 solution.cntcycles += 1;
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
151
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
152 if (options.havenonnegative)
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
153 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
154 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
155 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
156
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
157 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
158 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
159
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
160 ## 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
161 if (err <= 1)
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
162
20636
6e81f4b37e13 Performance improvements for ODE functions.
Rik <rik@octave.org>
parents: 20634
diff changeset
163 solution.cntloop += 1;
6e81f4b37e13 Performance improvements for ODE functions.
Rik <rik@octave.org>
parents: 20634
diff changeset
164 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
165 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
166 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
167
31263
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
168 istep++;
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_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
170 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
171 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
172
31263
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
173 ## 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
174 if (have_EventFcn)
31417
7286327ec4b6 Improve "Events" location in ode45, ode23, ode23s (bug #63162)
Ken Marek <marek_ka@mercer.edu>
parents: 31416
diff changeset
175 solution.event = ode_event_handler ([], t_new, x_new, ...
7286327ec4b6 Improve "Events" location in ode45, ode23, ode23s (bug #63162)
Ken Marek <marek_ka@mercer.edu>
parents: 31416
diff changeset
176 new_k_vals, [], []);
31263
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
177 ## 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
178 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
179 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
180 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
181 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
182 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
183 endif
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
184 endif
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
185
31263
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
186 ## 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
187 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
188 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
189 & (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
190 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
191 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
192 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
193 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
194 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
195 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
196 runge_kutta_interpolate (order, [t_old t_new], [x_old x_new], ...
31416
33d5c1c41bbc Clean up unused portions of runge_kutta_interpolate (bug #63162)
Ken Marek <marek_ka@mercer.edu>
parents: 31414
diff changeset
197 output_t(t_caught), new_k_vals);
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], ...
31416
33d5c1c41bbc Clean up unused portions of runge_kutta_interpolate (bug #63162)
Ken Marek <marek_ka@mercer.edu>
parents: 31414
diff changeset
213 tadd, new_k_vals);
31263
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
214 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
215 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
216 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
217 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
218 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
219 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
220 output_t(iout) = ode_t(istep);
31548
c8ad083a5802 maint: Clean up m-files before Octave 8.1 release.
Rik <rik@octave.org>
parents: 31417
diff changeset
221 endif
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
222
31263
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
223 ## 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
224 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
225 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
226 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
227 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
228 xadd = xadd(options.OutputSel, :);
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
229 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
230 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
231 [], 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
232 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
233 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
234 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
235 endif
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 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
238 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
239 endif
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
240
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 ## move to next time-step
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
243 t_old = t_new;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
244 x_old = x_new;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
245 k_vals = new_k_vals;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
246
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
247 else # error condition
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
248
20636
6e81f4b37e13 Performance improvements for ODE functions.
Rik <rik@octave.org>
parents: 20634
diff changeset
249 ireject += 1;
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
250
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
251 ## 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
252 ## value has been found.
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
253 if (ireject >= 5_000)
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
254 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
255 " 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
256 " 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
257 " 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
258 " 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
259 " 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
260 " 'MaxStep' with the command 'odeset'.\n"],
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
261 t_old, tspan(end));
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
262 endif
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
263
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
264 endif
20602
756b052037fb avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20601
diff changeset
265
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
266 ## 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
267 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
268 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
269 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
270 if (! (abs (dt) > eps (ode_t(end))))
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
271 break;
22538
21c89e691804 Stop adaptive ODE integration if timestep gets too small.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 22323
diff changeset
272 endif
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
273
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
274 ## 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
275 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
276
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
277 endwhile
20602
756b052037fb avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20601
diff changeset
278
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
279 ## 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
280 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
281 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
282 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
283 [" 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
284 " 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
285 " 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
286 " 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
287 " 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
288 " 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
289 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
290 endif
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
20601
43822bda4f65 fix indexing bug introduced with a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20600
diff changeset
293 ## 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
294 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
295 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
296 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
297 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
298
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
299 endfunction