Mercurial > octave

annotate scripts/ode/private/integrate_adaptive.m @ 22626:869c02fde46c stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Further clean-up of ode functions. * scripts/ode/private/AbsRel_Norm.m: Renamed to AbsRel_norm.m * scripts/ode/module.mk: Add AbsRel_norm.m to build system. * ode23.m, ode45.m: Remove extra comma from Copyright statement. Add ode45 to @seealso links in docstring. Add FIXME notes for questionable code. Use numel instead of length for clarity. Use space after function name and opening parenthesis. Wrap long lines to less than 80 characters. Change odemergeopts function call to match new prototype. Use single quotes to simplify strings that contain double quotes. Use 2-space indent in %!function blocks. * odeget.m: Remove extra comma from Copyright statement. Use double quote in preference to single quote. Delete whitespace at end of lines. * odeplot.m: Re-write docstring. Include @seealso links in docstring. Declare all persistent variables in a single declaration. Use in-place += operator for efficiency. Add FIXME notes for questionable code. * odeset.m: Remove extra comma from Copyright statement. Add additional calling form with 1 output and 0 inputs to docstring. Add FIXME notes for questionable code. Delete whitespace at end of lines. * odeset.m(print_options): Use single quotes to simplify strings with double quotes. Put default value of option first in list. * integrate_adaptive.m: Wrap long lines < 80 characters. Delete whitespace at end of lines. Correct indentation of declared values after '='. * kahan.m: Reise docstring. * ode_event_handler.m: Use retval in docstring to match functin prototype. Delete unnecessary comments. * odedefaults.m: Remove extra comma from Copyright statement. Match variable names in docstring to function prototype. Use space after function name and before opening parenthesis. Delete whitespace at end of lines. * odemergeopts.m: Remove extra comma from Copyright statement. Delete extra space in function prototype. Change function prototype to have "caller" as first argument to match rest of Octave. * runge_kutta_23.m: Clean up declaration of persistent variables. Correct misspellings in comments. * runge_kutta_45_dorpri.m: Put description of input arguments before output arguments in docstring. Clean up declaration of persistent variables. * runge_kutta_interpolate.m: Use double quotes in preference to single quotes. Eliminate line continuations for code that would fit on a single line. Remove obsolete code that calls non-existent functions. Capitalize Hermite in comments. Cleanup declaration of persistent variables. * starting_stepsize.m: Replace calls to AbsRel_Norm with AbsRel_norm.
author Rik <rik@octave.org>
date Sat, 15 Oct 2016 22:39:29 -0700
parents 177e0c71bcc0
children 7efa2d0e22c9
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
22323
bac0d6f07a3e maint: Update copyright notices for 2016.
John W. Eaton <jwe@octave.org>
parents: 22299
diff changeset
1 ## Copyright (C) 2016 Carlo de Falco
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
2 ## Copyright (C) 2013 Roberto Porcu' <roberto.porcu@polimi.it>
20536
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
3 ##
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
4 ## 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
5 ##
20536
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
6 ## Octave is free software; you can redistribute it and/or modify it
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
7 ## under the terms of the GNU General Public License as published by
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
8 ## the Free Software Foundation; either version 3 of the License, or (at
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
9 ## 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
10 ##
20536
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
11 ## 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
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
14 ## 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
15 ##
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
16 ## 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
17 ## along with Octave; see the file COPYING. If not, see
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
18 ## <http://www.gnu.org/licenses/>.
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 ## -*- texinfo -*-
20852
516bb87ea72e 2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents: 20847
diff changeset
21 ## @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
22 ##
20548
25623ef2ff4f doc: Rewrite docstrings for ode* family of functions.
Rik <rik@octave.org>
parents: 20543
diff changeset
23 ## 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
24 ## 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
25 ## timestep.
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
26 ##
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
27 ## 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
28 ## 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
29 ## @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
30 ## 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
31 ## 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
32 ## 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
33 ##
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
34 ## 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
35 ## 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
36 ## 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
37 ##
20548
25623ef2ff4f doc: Rewrite docstrings for ode* family of functions.
Rik <rik@octave.org>
parents: 20543
diff changeset
38 ## 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
39 ## to compute the adaptive timesteps.
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
40 ##
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
41 ## 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
42 ## defines the ODE:
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
43 ##
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
44 ## @ifhtml
22299
9fc91bb2aec3 doc: grammarcheck documentation for 4.2 release.
Rik <rik@octave.org>
parents: 22123
diff changeset
45 ##
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
46 ## @example
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
47 ## @math{y' = f(t,y)}
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
48 ## @end example
22299
9fc91bb2aec3 doc: grammarcheck documentation for 4.2 release.
Rik <rik@octave.org>
parents: 22123
diff changeset
49 ##
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
50 ## @end ifhtml
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
51 ## @ifnothtml
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
52 ## @math{y' = f(t,y)}.
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
53 ## @end ifnothtml
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
54 ##
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
55 ## 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
56 ## 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
57 ## 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
58 ##
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
59 ## 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
60 ## 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
61 ## stepper.
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
62 ##
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
63 ## @end deftypefn
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
64
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
65 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
66 options)
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
67
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
68 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
69
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
70 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
71 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
72
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
73 ## 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
74 dt = options.InitialStep;
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
75 if (isempty (dt))
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
76 dt = starting_stepsize (order, func, t, x,
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
77 options.AbsTol, options.RelTol,
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
78 strcmp (options.NormControl, "on"),
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
79 options.funarguments);
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
80 endif
20602
756b052037fb avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20601
diff changeset
81
22605
177e0c71bcc0 make sure the additional function arguments are always passed to odefun.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 22596
diff changeset
82 dir = options.direction;
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
83 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
84
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
85 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
86
20552
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20548
diff changeset
87 ## 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
88 facmin = 0.8;
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
89 facmax = 1.5;
20552
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20548
diff changeset
90 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
91
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
92 ## Initialize the OutputFcn
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
93 if (options.haveoutputfunction)
22593
dba5074bdc79 simplify options management in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents: 22538
diff changeset
94 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
95 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
96 else
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
97 solution.retout = x;
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
98 endif
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
99 feval (options.OutputFcn, tspan, solution.retout,
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
100 "init", options.funarguments{:});
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
101 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
102
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
103 ## Initialize the EventFcn
22593
dba5074bdc79 simplify options management in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents: 22538
diff changeset
104 if (! isempty (options.Events))
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
105 ode_event_handler (options.Events, tspan(end), x,
22538
21c89e691804 Stop adaptive ODE integration if timestep gets too small.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 22323
diff changeset
106 "init", 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
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
109 if (options.havenonnegative)
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
110 nn = options.NonNegative;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
111 endif
20602
756b052037fb avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20601
diff changeset
112
21955
0f3e875d9078 Fix statistics on solution solving for ode45, ode23 (bug #48243).
jcorno <jacopo.corno@gmail.com>
parents: 21546
diff changeset
113 solution.cntloop = 0;
0f3e875d9078 Fix statistics on solution solving for ode45, ode23 (bug #48243).
jcorno <jacopo.corno@gmail.com>
parents: 21546
diff changeset
114 solution.cntcycles = 0;
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
115 solution.cntsave = 2;
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
116 solution.unhandledtermination = true;
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
117 ireject = 0;
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
118
20602
756b052037fb avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20601
diff changeset
119 k_vals = [];
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
120 iout = istep = 1;
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
121
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
122 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
123
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
124 ## 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
125 [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
126 [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
127 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
128
20636
6e81f4b37e13 Performance improvements for ODE functions.
Rik <rik@octave.org>
parents: 20634
diff changeset
129 solution.cntcycles += 1;
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
130
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
131 if (options.havenonnegative)
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
132 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
133 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
134 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
135
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
136 ## FIXME: Take strcmp out of while loop and calculate just once
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
137 err = AbsRel_norm (x_new, x_old, options.AbsTol, options.RelTol,
22593
dba5074bdc79 simplify options management in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents: 22538
diff changeset
138 strcmp (options.NormControl, "on"), x_est);
20602
756b052037fb avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20601
diff changeset
139
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
140 ## 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
141 if (err <= 1)
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
142
20636
6e81f4b37e13 Performance improvements for ODE functions.
Rik <rik@octave.org>
parents: 20634
diff changeset
143 solution.cntloop += 1;
6e81f4b37e13 Performance improvements for ODE functions.
Rik <rik@octave.org>
parents: 20634
diff changeset
144 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
145
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
146 ## if output time steps are fixed
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
147 if (fixed_times)
20543
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20536
diff changeset
148
20604
a260a6acb70f fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20602
diff changeset
149 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
150 & (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
151 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
152
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
153 if (! isempty (t_caught))
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
154 t(t_caught) = tspan(t_caught);
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
155 iout = max (t_caught);
20601
43822bda4f65 fix indexing bug introduced with a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20600
diff changeset
156 x(:, t_caught) = ...
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
157 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
158 t(t_caught), new_k_vals, dt, func,
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
159 options.funarguments);
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
160
20636
6e81f4b37e13 Performance improvements for ODE functions.
Rik <rik@octave.org>
parents: 20634
diff changeset
161 istep += 1;
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
162
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
163 ## 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
164 ## Stop integration if eventbreak is true.
22593
dba5074bdc79 simplify options management in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents: 22538
diff changeset
165 if (! isempty (options.Events))
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
166 break_loop = false;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
167 for idenseout = 1:numel (t_caught)
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
168 id = t_caught(idenseout);
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
169 td = t(id);
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
170 solution.event = ...
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
171 ode_event_handler (options.Events, t(id), x(:, id), [],
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
172 options.funarguments{:});
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
173 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
174 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
175 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
176 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
177 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
178 solution.unhandledtermination = false;
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
179 break_loop = true;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
180 break;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
181 endif
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
182 endfor
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
183 if (break_loop)
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
184 break;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
185 endif
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
186 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
187
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
188 ## Call OutputFcn only if a valid result has been found.
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
189 ## Stop integration if function returns false.
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
190 if (options.haveoutputfunction)
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
191 cnt = options.Refine + 1;
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
192 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
193 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
194 [x_old, x(:, t_caught), x_new] .',
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
195 approxtime, "linear") .';
22593
dba5074bdc79 simplify options management in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents: 22538
diff changeset
196 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
197 approxvals = approxvals(options.OutputSel, :);
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
198 endif
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
199 for ii = 1:numel (approxtime)
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
200 pltret = feval (options.OutputFcn, approxtime(ii),
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
201 approxvals(:, ii), [],
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
202 options.funarguments{:});
20604
a260a6acb70f fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20602
diff changeset
203 endfor
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
204 if (pltret) # Leave main loop
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
205 solution.unhandledtermination = false;
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
206 break;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
207 endif
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
208 endif
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
209
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
210 endif
20602
756b052037fb avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20601
diff changeset
211
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
212 else # not fixed times
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
213
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
214 t(++istep) = t_new;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
215 x(:, istep) = x_new;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
216 iout = istep;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
217
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
218 ## 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
219 ## Stop integration if eventbreak is true.
22593
dba5074bdc79 simplify options management in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents: 22538
diff changeset
220 if (! isempty (options.Events))
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
221 solution.event = ...
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
222 ode_event_handler (options.Events, t(istep), x(:, istep), [],
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
223 options.funarguments{:});
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
224 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
225 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
226 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
227 solution.unhandledtermination = false;
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
228 break;
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
229 endif
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
230 endif
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
231
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
232 ## Call OutputFcn only if a valid result has been found.
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
233 ## Stop integration if function returns false.
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
234 if (options.haveoutputfunction)
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
235 cnt = options.Refine + 1;
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
236 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
237 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
238 [x_old, x_new] .',
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
239 approxtime, "linear") .';
22593
dba5074bdc79 simplify options management in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents: 22538
diff changeset
240 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
241 approxvals = approxvals(options.OutputSel, :);
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
242 endif
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
243 for ii = 1:numel (approxtime)
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
244 pltret = feval (options.OutputFcn, approxtime(ii),
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
245 approxvals(:, ii), [], options.funarguments{:});
20604
a260a6acb70f fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20602
diff changeset
246 endfor
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
247 if (pltret) # Leave main loop
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
248 solution.unhandledtermination = false;
20552
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20548
diff changeset
249 break;
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
250 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
251 endif
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
252
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
253 endif
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
254
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
255 ## move to next time-step
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
256 t_old = t_new;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
257 x_old = x_new;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
258 k_vals = new_k_vals;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
259
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
260 else # error condition
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
261
20636
6e81f4b37e13 Performance improvements for ODE functions.
Rik <rik@octave.org>
parents: 20634
diff changeset
262 ireject += 1;
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
263
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
264 ## 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
265 ## value has been found.
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
266 if (ireject >= 5_000)
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20604
diff changeset
267 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
268 " 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
269 " 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
270 " 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
271 " 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
272 " 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
273 " 'MaxStep' with the command 'odeset'.\n"],
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
274 t_old, tspan(end));
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
275 endif
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
276
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
277 endif
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 ## 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
280 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
281 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
282 dt = dir * min (abs (dt), options.MaxStep);
22538
21c89e691804 Stop adaptive ODE integration if timestep gets too small.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 22323
diff changeset
283 if (! (abs (dt) > eps (t (end))))
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
284 break;
22538
21c89e691804 Stop adaptive ODE integration if timestep gets too small.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 22323
diff changeset
285 endif
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
286
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22605
diff changeset
287 ## 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
288 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
289
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
290 endwhile
20602
756b052037fb avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20601
diff changeset
291
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
292 ## 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
293 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
294 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
295 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
296 [" 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
297 " 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
298 " 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
299 " 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
300 " 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
301 " 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
302 t(end), tspan(end));
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
303 else
20601
43822bda4f65 fix indexing bug introduced with a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20600
diff changeset
304 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
305 ["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
306 " 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
307 " 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
308 " 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
309 " 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
310 t(end), tspan(end));
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
311 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
312 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
313
20601
43822bda4f65 fix indexing bug introduced with a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20600
diff changeset
314 ## Set up return structure
43822bda4f65 fix indexing bug introduced with a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20600
diff changeset
315 solution.t = t(:);
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
316 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
317
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
318 endfunction
20552
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20548
diff changeset
319