Mercurial > octave

annotate scripts/ode/private/integrate_adaptive.m @ 29358:0a5b15007766 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
update Octave Project Developers copyright for the new year In files that have the "Octave Project Developers" copyright notice, update for 2021.
author John W. Eaton <jwe@octave.org>
date Wed, 10 Feb 2021 09:52:15 -0500
parents bd51beb6205e
children 796f54d4ddbf
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 ##
29358
0a5b15007766 update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents: 27923
diff changeset
3 ## Copyright (C) 2013-2021 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 -*-
20852
516bb87ea72e 2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents: 20847
diff changeset
27 ## @deftypefn {} {@var{solution} =} integrate_adaptive (@var{@@stepper}, @var{order}, @var{@@func}, @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
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
71 function solution = integrate_adaptive (stepper, order, func, tspan, x0,
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
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
76 t_new = t_old = t = tspan(1);
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
77 x_new = x_old = 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))
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
82 dt = starting_stepsize (order, func, t, x,
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
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
98 ## Initialize the OutputFcn
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
99 if (options.haveoutputfunction)
22593
dba5074bdc79 simplify options management in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents: 22538
diff changeset
100 if (! isempty (options.OutputSel))
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
101 solution.retout = 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
102 else
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
103 solution.retout = x;
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
104 endif
22639
7efa2d0e22c9 More Matlab-compatible implementation of OutputFcn. Clean up odeplot.m.
Rik <rik@octave.org>
parents: 22626
diff changeset
105 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
106 options.funarguments{:});
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
107 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
108
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
109 ## 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
110 have_EventFcn = false;
22593
dba5074bdc79 simplify options management in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents: 22538
diff changeset
111 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
112 have_EventFcn = true;
22894
544dac72d414 Initialize ODE Event function with start time (bug #49846).
Rik <rik@octave.org>
parents: 22669
diff changeset
113 ode_event_handler (options.Events, tspan(1), x,
22538
21c89e691804 Stop adaptive ODE integration if timestep gets too small.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 22323
diff changeset
114 "init", options.funarguments{:});
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
115 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
116
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
117 if (options.havenonnegative)
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
118 nn = options.NonNegative;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
119 endif
20602
756b052037fb avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20601
diff changeset
120
21955
0f3e875d9078 Fix statistics on solution solving for ode45, ode23 (bug #48243).
jcorno <jacopo.corno@gmail.com>
parents: 21546
diff changeset
121 solution.cntloop = 0;
0f3e875d9078 Fix statistics on solution solving for ode45, ode23 (bug #48243).
jcorno <jacopo.corno@gmail.com>
parents: 21546
diff changeset
122 solution.cntcycles = 0;
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
123 solution.cntsave = 2;
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
124 solution.unhandledtermination = true;
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
125 ireject = 0;
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
126
22645
34bf558de127 integrate_adaptive.m: Take strcmp call out of for loop for performance.
Rik <rik@octave.org>
parents: 22639
diff changeset
127 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
128 k_vals = [];
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
129 iout = istep = 1;
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
130
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
131 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
132
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
133 ## 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
134 [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
135 [t_new, x_new, x_est, new_k_vals] = ...
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
136 stepper (func, 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
137
20636
6e81f4b37e13 Performance improvements for ODE functions.
Rik <rik@octave.org>
parents: 20634
diff changeset
138 solution.cntcycles += 1;
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
139
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
140 if (options.havenonnegative)
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
141 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
142 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
143 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
144
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
145 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
146 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
147
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
148 ## 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
149 if (err <= 1)
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
150
20636
6e81f4b37e13 Performance improvements for ODE functions.
Rik <rik@octave.org>
parents: 20634
diff changeset
151 solution.cntloop += 1;
6e81f4b37e13 Performance improvements for ODE functions.
Rik <rik@octave.org>
parents: 20634
diff changeset
152 ireject = 0; # Clear reject counter
20602
756b052037fb avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20601
diff changeset
153
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
154 ## if output time steps are fixed
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
155 if (fixed_times)
20543
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20536
diff changeset
156
20604
a260a6acb70f fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20602
diff changeset
157 t_caught = find ((dir * tspan(iout:end) > dir * t_old)
a260a6acb70f fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20602
diff changeset
158 & (dir * tspan(iout:end) <= dir * t_new));
20601
43822bda4f65 fix indexing bug introduced with a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20600
diff changeset
159 t_caught = t_caught + iout - 1;
20602
756b052037fb avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20601
diff changeset
160
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
161 if (! isempty (t_caught))
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
162 t(t_caught) = tspan(t_caught);
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
163 iout = max (t_caught);
20601
43822bda4f65 fix indexing bug introduced with a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20600
diff changeset
164 x(:, t_caught) = ...
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
165 runge_kutta_interpolate (order, [t_old t_new], [x_old x_new],
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
166 t(t_caught), new_k_vals, dt, func,
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
167 options.funarguments);
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
168
20636
6e81f4b37e13 Performance improvements for ODE functions.
Rik <rik@octave.org>
parents: 20634
diff changeset
169 istep += 1;
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
170
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
171 ## Call Events function only if a valid result has been found.
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
172 ## Stop integration if eventbreak is true.
23079
17dc6c7ef427 Ignore ODE event function at t==0 for Matlab compatibility (bug #49919).
Rik <rik@octave.org>
parents: 22916
diff changeset
173 if (have_EventFcn)
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
174 break_loop = false;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
175 for idenseout = 1:numel (t_caught)
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
176 id = t_caught(idenseout);
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
177 td = t(id);
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
178 solution.event = ...
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
179 ode_event_handler (options.Events, t(id), x(:, id), [],
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
180 options.funarguments{:});
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
181 if (! isempty (solution.event{1}) && solution.event{1} == 1)
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
182 t(id) = solution.event{3}(end);
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
183 t = t(1:id);
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
184 x(:, id) = solution.event{4}(end, :).';
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
185 x = x(:,1:id);
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
186 solution.unhandledtermination = false;
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
187 break_loop = true;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
188 break;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
189 endif
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
190 endfor
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
191 if (break_loop)
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
192 break;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
193 endif
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
194 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
195
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
196 ## Call OutputFcn only if a valid result has been found.
22639
7efa2d0e22c9 More Matlab-compatible implementation of OutputFcn. Clean up odeplot.m.
Rik <rik@octave.org>
parents: 22626
diff changeset
197 ## Stop integration if function returns true.
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
198 if (options.haveoutputfunction)
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
199 cnt = options.Refine + 1;
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
200 approxtime = linspace (t_old, t_new, cnt);
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
201 approxvals = interp1 ([t_old, t(t_caught), t_new],
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
202 [x_old, x(:, t_caught), x_new] .',
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
203 approxtime, "linear") .';
22669
bcfcbb5653e6 integrate_adaptive.m: Further refinement to orientation of approxvals (cset 7458bacc631a).
Rik <rik@octave.org>
parents: 22650
diff changeset
204 if (isvector (approxvals) && ! isscalar (approxtime))
22650
7458bacc631a integrate_adaptive.m: Fix orientation of approxvals when y is a vector.
Rik <rik@octave.org>
parents: 22645
diff changeset
205 approxvals = approxvals.';
7458bacc631a integrate_adaptive.m: Fix orientation of approxvals when y is a vector.
Rik <rik@octave.org>
parents: 22645
diff changeset
206 endif
22593
dba5074bdc79 simplify options management in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents: 22538
diff changeset
207 if (! isempty (options.OutputSel))
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
208 approxvals = approxvals(options.OutputSel, :);
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
209 endif
22639
7efa2d0e22c9 More Matlab-compatible implementation of OutputFcn. Clean up odeplot.m.
Rik <rik@octave.org>
parents: 22626
diff changeset
210 stop_solve = false;
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
211 for ii = 1:numel (approxtime)
22639
7efa2d0e22c9 More Matlab-compatible implementation of OutputFcn. Clean up odeplot.m.
Rik <rik@octave.org>
parents: 22626
diff changeset
212 stop_solve = feval (options.OutputFcn,
7efa2d0e22c9 More Matlab-compatible implementation of OutputFcn. Clean up odeplot.m.
Rik <rik@octave.org>
parents: 22626
diff changeset
213 approxtime(ii), approxvals(:, ii), [],
7efa2d0e22c9 More Matlab-compatible implementation of OutputFcn. Clean up odeplot.m.
Rik <rik@octave.org>
parents: 22626
diff changeset
214 options.funarguments{:});
7efa2d0e22c9 More Matlab-compatible implementation of OutputFcn. Clean up odeplot.m.
Rik <rik@octave.org>
parents: 22626
diff changeset
215 if (stop_solve)
7efa2d0e22c9 More Matlab-compatible implementation of OutputFcn. Clean up odeplot.m.
Rik <rik@octave.org>
parents: 22626
diff changeset
216 break; # break from inner loop
7efa2d0e22c9 More Matlab-compatible implementation of OutputFcn. Clean up odeplot.m.
Rik <rik@octave.org>
parents: 22626
diff changeset
217 endif
20604
a260a6acb70f fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20602
diff changeset
218 endfor
22639
7efa2d0e22c9 More Matlab-compatible implementation of OutputFcn. Clean up odeplot.m.
Rik <rik@octave.org>
parents: 22626
diff changeset
219 if (stop_solve) # Leave main loop
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
220 solution.unhandledtermination = false;
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
221 break;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
222 endif
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
223 endif
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
224
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
225 endif
20602
756b052037fb avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20601
diff changeset
226
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
227 else # not fixed times
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
228
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
229 t(++istep) = t_new;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
230 x(:, istep) = x_new;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
231 iout = istep;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
232
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
233 ## Call Events function only if a valid result has been found.
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
234 ## Stop integration if eventbreak is true.
23079
17dc6c7ef427 Ignore ODE event function at t==0 for Matlab compatibility (bug #49919).
Rik <rik@octave.org>
parents: 22916
diff changeset
235 if (have_EventFcn)
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
236 solution.event = ...
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
237 ode_event_handler (options.Events, t(istep), x(:, istep), [],
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
238 options.funarguments{:});
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
239 if (! isempty (solution.event{1}) && solution.event{1} == 1)
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
240 t(istep) = solution.event{3}(end);
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
241 x(:, istep) = solution.event{4}(end, :).';
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
242 solution.unhandledtermination = false;
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
243 break;
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
244 endif
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
245 endif
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
246
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
247 ## Call OutputFcn only if a valid result has been found.
22639
7efa2d0e22c9 More Matlab-compatible implementation of OutputFcn. Clean up odeplot.m.
Rik <rik@octave.org>
parents: 22626
diff changeset
248 ## Stop integration if function returns true.
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
249 if (options.haveoutputfunction)
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
250 cnt = options.Refine + 1;
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
251 approxtime = linspace (t_old, t_new, cnt);
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
252 approxvals = interp1 ([t_old, t_new],
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
253 [x_old, x_new] .',
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
254 approxtime, "linear") .';
22669
bcfcbb5653e6 integrate_adaptive.m: Further refinement to orientation of approxvals (cset 7458bacc631a).
Rik <rik@octave.org>
parents: 22650
diff changeset
255 if (isvector (approxvals) && ! isscalar (approxtime))
22650
7458bacc631a integrate_adaptive.m: Fix orientation of approxvals when y is a vector.
Rik <rik@octave.org>
parents: 22645
diff changeset
256 approxvals = approxvals.';
7458bacc631a integrate_adaptive.m: Fix orientation of approxvals when y is a vector.
Rik <rik@octave.org>
parents: 22645
diff changeset
257 endif
22593
dba5074bdc79 simplify options management in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents: 22538
diff changeset
258 if (! isempty (options.OutputSel))
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
259 approxvals = approxvals(options.OutputSel, :);
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
260 endif
22639
7efa2d0e22c9 More Matlab-compatible implementation of OutputFcn. Clean up odeplot.m.
Rik <rik@octave.org>
parents: 22626
diff changeset
261 stop_solve = false;
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
262 for ii = 1:numel (approxtime)
22639
7efa2d0e22c9 More Matlab-compatible implementation of OutputFcn. Clean up odeplot.m.
Rik <rik@octave.org>
parents: 22626
diff changeset
263 stop_solve = feval (options.OutputFcn,
7efa2d0e22c9 More Matlab-compatible implementation of OutputFcn. Clean up odeplot.m.
Rik <rik@octave.org>
parents: 22626
diff changeset
264 approxtime(ii), approxvals(:, ii), [],
7efa2d0e22c9 More Matlab-compatible implementation of OutputFcn. Clean up odeplot.m.
Rik <rik@octave.org>
parents: 22626
diff changeset
265 options.funarguments{:});
7efa2d0e22c9 More Matlab-compatible implementation of OutputFcn. Clean up odeplot.m.
Rik <rik@octave.org>
parents: 22626
diff changeset
266 if (stop_solve)
7efa2d0e22c9 More Matlab-compatible implementation of OutputFcn. Clean up odeplot.m.
Rik <rik@octave.org>
parents: 22626
diff changeset
267 break; # break from inner loop
7efa2d0e22c9 More Matlab-compatible implementation of OutputFcn. Clean up odeplot.m.
Rik <rik@octave.org>
parents: 22626
diff changeset
268 endif
20604
a260a6acb70f fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20602
diff changeset
269 endfor
22639
7efa2d0e22c9 More Matlab-compatible implementation of OutputFcn. Clean up odeplot.m.
Rik <rik@octave.org>
parents: 22626
diff changeset
270 if (stop_solve) # Leave main loop
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
271 solution.unhandledtermination = false;
20552
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20548
diff changeset
272 break;
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
273 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
274 endif
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
275
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
276 endif
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
277
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
278 ## move to next time-step
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
279 t_old = t_new;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
280 x_old = x_new;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
281 k_vals = new_k_vals;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
282
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
283 else # error condition
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
284
20636
6e81f4b37e13 Performance improvements for ODE functions.
Rik <rik@octave.org>
parents: 20634
diff changeset
285 ireject += 1;
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
286
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
287 ## 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
288 ## value has been found.
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
289 if (ireject >= 5_000)
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
290 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
291 " 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
292 " 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
293 " 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
294 " 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
295 " 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
296 " 'MaxStep' with the command 'odeset'.\n"],
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
297 t_old, tspan(end));
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
298 endif
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
299
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
300 endif
20602
756b052037fb avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20601
diff changeset
301
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
302 ## 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
303 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
304 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
305 dt = dir * min (abs (dt), options.MaxStep);
22645
34bf558de127 integrate_adaptive.m: Take strcmp call out of for loop for performance.
Rik <rik@octave.org>
parents: 22639
diff changeset
306 if (! (abs (dt) > eps (t(end))))
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
307 break;
22538
21c89e691804 Stop adaptive ODE integration if timestep gets too small.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 22323
diff changeset
308 endif
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
309
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
310 ## 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
311 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
312
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
313 endwhile
20602
756b052037fb avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20601
diff changeset
314
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
315 ## Check if integration of the ode has been successful
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
316 if (dir * 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
317 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
318 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
319 [" 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
320 " 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
321 " 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
322 " 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
323 " 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
324 " and/or 'MaxStep' with the command 'odeset'."],
21c89e691804 Stop adaptive ODE integration if timestep gets too small.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 22323
diff changeset
325 t(end), tspan(end));
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
326 else
20601
43822bda4f65 fix indexing bug introduced with a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20600
diff changeset
327 warning ("integrate_adaptive:unexpected_termination",
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
328 ["Solver was stopped by a call of 'break'", ...
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
329 " in the main iteration loop at time", ...
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
330 " t = %f before the endpoint at tend = %f was reached. ", ...
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
331 " This may happen because the @odeplot function", ...
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
332 " returned 'true' or the @event function returned 'true'."],
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
333 t(end), tspan(end));
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
334 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
335 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
336
20601
43822bda4f65 fix indexing bug introduced with a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20600
diff changeset
337 ## Set up return structure
43822bda4f65 fix indexing bug introduced with a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20600
diff changeset
338 solution.t = t(:);
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
339 solution.x = x.';
20602
756b052037fb avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20601
diff changeset
340
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
341 endfunction