Mercurial > octave-nkf

annotate scripts/ode/private/integrate_adaptive.m @ 20635:a22d8a2eb0e5

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
fix adaptive strategy in ode solvers. * script/ode/ode45.m: remove unused option OutputSave * script/ode/private/integrate_adaptive.m: rewrite algorithm to be more compatible. * script/ode/private/runge_kutta_45_dorpri.m: use kahan summation for time increment.
author Carlo de Falco <carlo.defalco@polimi.it>
date Sun, 11 Oct 2015 18:44:58 +0200
parents 87b557ee8e5d
children 43822bda4f65
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
20596
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20586
diff changeset
1 ## Copyright (C) 2015, Carlo de Falco
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
2 ## Copyright (C) 2013, Roberto Porcu' <roberto.porcu@polimi.it>
20570
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20568
diff changeset
3 ##
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20568
diff changeset
4 ## This file is part of Octave.
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
5 ##
20570
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20568
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: 20568
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: 20568
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: 20568
diff changeset
9 ## your option) any later version.
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
10 ##
20570
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20568
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: 20568
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: 20568
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: 20568
diff changeset
14 ## General Public License for more details.
20568
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
20570
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20568
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: 20568
diff changeset
18 ## <http://www.gnu.org/licenses/>.
20568
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 -*-
20580
25623ef2ff4f doc: Rewrite docstrings for ode* family of functions.
Rik <rik@octave.org>
parents: 20575
diff changeset
21 ## @deftypefn {Function File} {[@var{t}, @var{y}] =} integrate_adaptive (@var{@@stepper}, @var{order}, @var{@@fun}, @var{tspan}, @var{x0}, @var{options})
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
22 ##
20580
25623ef2ff4f doc: Rewrite docstrings for ode* family of functions.
Rik <rik@octave.org>
parents: 20575
diff changeset
23 ## This function file can be called by an ODE solver function in order to
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
24 ## integrate the set of ODEs on the interval @var{[t0,t1]} with an
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
25 ## adaptive timestep.
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
26 ##
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
27 ## This function must be called with two output arguments: @var{t} and @var{y}.
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
28 ## Variable @var{t} is a column vector and contains the time stamps, instead
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
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
30 ## of the problem and the rows number is the same of @var{t} rows number so
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
31 ## that each row of @var{y} contains the values of all unknowns at the time
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
32 ## value contained in the corresponding row in @var{t}.
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
33 ##
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
34 ## The first input argument must be a function_handle or an inline function
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
35 ## representing the stepper, that is the function responsible for step-by-step
20580
25623ef2ff4f doc: Rewrite docstrings for ode* family of functions.
Rik <rik@octave.org>
parents: 20575
diff changeset
36 ## integration. This function discriminates one method from the others.
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
37 ##
20580
25623ef2ff4f doc: Rewrite docstrings for ode* family of functions.
Rik <rik@octave.org>
parents: 20575
diff changeset
38 ## The second input argument is the order of the stepper. It is needed
20568
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 ##
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
41 ## The third input argument is a function_handle or an inline function that
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
42 ## defines the set of ODE:
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
43 ## @ifhtml
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
44 ## @example
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
45 ## @math{y' = f(t,y)}
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
46 ## @end example
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
47 ## @end ifhtml
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
48 ## @ifnothtml
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
49 ## @math{y' = f(t,y)}.
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
50 ## @end ifnothtml
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
51 ##
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
52 ## The fourth input argument is the time vector which defines integration
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
53 ## interval, that is @var{[tspan(1),tspan(end)]} and all the intermediate
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
54 ## elements are taken as times at which the solution is required.
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
55 ##
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
56 ## The fifth argument represents the initial conditions for the ODEs and the
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
57 ## last input argument contains some options that may be needed for the stepper.
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 ## @end deftypefn
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
60 ##
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
61 ## @seealso{integrate_const, integrate_n_steps}
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 function solution = integrate_adaptive (stepper, order, func, tspan, x0, options)
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
64
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
65 fixed_times = numel (tspan) > 2;
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
66
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
67 t_new = t_old = t = tspan(1);
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
68 x_new = x_old = x = x0(:);
20575
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20570
diff changeset
69
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
70 ## get first initial timestep
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
71 dt = starting_stepsize (order, func, t, x, options.AbsTol,
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
72 options.RelTol, options.vnormcontrol);
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
73 dt = odeget (options, "InitialStep", dt, "fast_not_empty");
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
74
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
75 dir = odeget (options, "vdirection", [], "fast");
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
76 dt = dir * min (abs (dt), options.MaxStep);
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
77
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
78 options.comp = 0.0;
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
79
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
80 ## Factor multiplying the stepsize guess
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
81 facmin = 0.8;
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
82 facmax = 1.5;
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
83 fac = 0.38^(1/(order+1)); # formula taken from Hairer
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
84
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
85 ## Initialize the OutputFcn
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
86 if (options.vhaveoutputfunction)
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
87 if (options.vhaveoutputselection)
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
88 solution.vretout = x(options.OutputSel,end);
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
89 else
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
90 solution.vretout = x;
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
91 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
92 feval (options.OutputFcn, tspan, solution.vretout,
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
93 "init", options.vfunarguments{:});
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
94 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
95
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
96 ## Initialize the EventFcn
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
97 if (options.vhaveeventfunction)
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
98 odepkg_event_handle (options.Events, tspan(end), x,
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
99 "init", options.vfunarguments{:});
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
100 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
101
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
102 if (options.vhavenonnegative)
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
103 nn = options.NonNegative;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
104 endif
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
105
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
106 solution.vcntloop = 2;
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
107 solution.vcntcycles = 1;
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
108 solution.vcntsave = 2;
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
109 solution.vunhandledtermination = true;
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
110 ireject = 0;
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
111
20575
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20570
diff changeset
112 k_vals = [];
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
113 iout = istep = 1;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
114 while (dir * t_old < dir * tspan(end))
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
115
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
116 ## 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: 20596
diff changeset
117 [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: 20596
diff changeset
118 [t_new, x_new, x_est, new_k_vals] = ...
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
119 stepper (func, t_old, x_old, dt, options, k_vals, t_new);
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
120
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
121 solution.vcntcycles++;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
122
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
123 if (options.vhavenonnegative)
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
124 x(nn,end) = abs (x(nn,end));
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
125 y(nn,end) = abs (y(nn,end));
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
126 y_est(nn,end) = abs (y_est(nn,end));
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
127 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
128
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
129 err = AbsRel_Norm (x_new, x_old, options.AbsTol, options.RelTol,
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
130 options.vnormcontrol, x_est);
20575
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20570
diff changeset
131
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
132 ## Accepted solution only if err <= 1.0
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
133 if (err <= 1)
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
134
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
135 solution.vcntloop++;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
136 ireject = 0;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
137
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
138 ## if output time steps are fixed
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
139 if (fixed_times)
20575
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20570
diff changeset
140
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
141 t_caught = find ((tspan(iout:end) > t_old)
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
142 & (tspan(iout:end) <= t_new));
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
143 if (! isempty (t_caught))
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
144 t(t_caught) = tspan(t_caught);
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
145 iout = max (t_caught);
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
146 x(:, t_caught) = interpolate ([t_old, t_new], [x_old, x_new],
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
147 t(t_caught));
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
148
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
149 istep++;
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
150
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
151 if (options.vhaveeventfunction)
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
152 ## Call event on each dense output timestep.
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
153 ## Stop integration if veventbreak is true
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
154 break_loop = false;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
155 for idenseout = 1:numel (t_caught)
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
156 id = t_caught(idenseout);
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
157 td = t(id);
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
158 solution.vevent = ...
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
159 odepkg_event_handle (options.Events, t(id), x(:, id), [],
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
160 options.vfunarguments{:});
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
161 if (! isempty (solution.vevent{1})
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
162 && solution.vevent{1} == 1)
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
163 t(id) = solution.vevent{3}(end);
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
164 t = t(1:id);
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
165 x(:, id) = solution.vevent{4}(end, :).';
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
166 x = x(:,1:id);
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
167 solution.vunhandledtermination = false;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
168 break_loop = true;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
169 break;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
170 endif
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
171 endfor
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
172 if (break_loop)
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
173 break;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
174 endif
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
175 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
176
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
177 ## Call plot. Stop integration if plot function
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
178 ## returns true.
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
179 if (options.vhaveoutputfunction)
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
180 vcnt = options.Refine + 1;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
181 vapproxtime = linspace (t_old, t_new, vcnt);
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
182 vapproxvals = interp1 ([t_old, t(t_caught), t_new],
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
183 [x_old, x(:, t_caught), x_new],
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
184 vapproxtime, 'linear');
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
185 if (options.vhaveoutputselection)
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
186 vapproxvals = vapproxvals(options.OutputSel);
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
187 endif
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
188 vpltret = feval (options.OutputFcn, vapproxtime,
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
189 vapproxvals, [], options.vfunarguments{:});
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
190 if (vpltret) # Leave main loop
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
191 solution.vunhandledtermination = false;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
192 break;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
193 endif
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
194 endif
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
195
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
196 endif
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
197
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
198 else
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
199
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
200 t(++istep) = t_new;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
201 x(:, istep) = x_new;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
202 iout = istep;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
203
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
204 ## Call event handler on new timestep.
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
205 ## Stop integration if veventbreak is true
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
206 if (options.vhaveeventfunction)
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
207 solution.vevent = ...
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
208 odepkg_event_handle (options.Events, t(istep), x(:, istep), [],
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
209 options.vfunarguments{:});
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
210 if (! isempty (solution.vevent{1})
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
211 && solution.vevent{1} == 1)
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
212 t(istep) = solution.vevent{3}(end);
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
213 x(:, istep) = solution.vevent{4}(end, :).';
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
214 solution.vunhandledtermination = false;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
215 break;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
216 endif
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
217 endif
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
218
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
219 ## Call plot. Stop integration if plot function
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
220 ## returns true.
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
221 if (options.vhaveoutputfunction)
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
222 vcnt = options.Refine + 1;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
223 vapproxtime = linspace (t_old, t_new, vcnt);
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
224 vapproxvals = interp1 ([t_old, t_new],
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
225 [x_old, x_new],
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
226 vapproxtime, 'linear');
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
227 if (options.vhaveoutputselection)
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
228 vapproxvals = vapproxvals(options.OutputSel);
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
229 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
230 vpltret = feval (options.OutputFcn, vapproxtime,
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
231 vapproxvals, [], options.vfunarguments{:});
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
232 if (vpltret) # Leave main loop
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
233 solution.vunhandledtermination = false;
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
234 break;
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
235 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
236 endif
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
237
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
238 endif
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
239
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
240 ## move to next time-step
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
241 t_old = t_new;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
242 x_old = x_new;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
243 k_vals = new_k_vals;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
244
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
245 solution.vcntloop = solution.vcntloop + 1;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
246 vcntiter = 0;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
247
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
248 else
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
249
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
250 ireject++;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
251
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
252 ## Stop solving because in the last 5000 steps no successful valid
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
253 ## value has been found
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
254 if (ireject >= 5000)
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
255 error (["integrate_adaptive: Solving has not been successful. ",
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
256 " The iterative integration loop exited at time",
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
257 " t = %f before endpoint at tend = %f was reached. ",
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
258 " This happened because the iterative integration loop",
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
259 " does not find a valid solution at this time",
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
260 " stamp. Try to reduce the value of 'InitialStep' and/or",
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
261 " 'MaxStep' with the command 'odeset'.\n"],
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
262 t_old, tspan(end));
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
263 endif
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
264
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
265 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
266
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
267 ## Compute next timestep, formula taken from Hairer
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
268 err += eps; # avoid divisions by zero
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
269 dt *= min (facmax, max (facmin, fac * (1 / err)^(1 / (order + 1))));
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
270 dt = dir * min (abs (dt), options.MaxStep);
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
271
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
272 endwhile
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
273
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
274 ## Check if integration of the ode has been successful
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
275 if (dir * t(end) < dir * tspan(end))
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
276 if (solution.vunhandledtermination == true)
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
277 error ("integrate_adaptive: InvalidArgument",
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
278 ["Solving has not been successful. The iterative",
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
279 " integration loop exited at time t = %f",
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
280 " before endpoint at tend = %f was reached. This may",
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
281 " happen if the stepsize grows too small. ",
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
282 " Try to reduce the value of 'InitialStep'",
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
283 " and/or 'MaxStep' with the command 'odeset'.\n"],
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
284 t(end), tspan(end));
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
285 else
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
286 warning ("integrate_adaptive: InvalidArgument",
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
287 ["Solver has been stopped by a call of 'break' in the main",
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
288 " iteration loop at time t = %f before endpoint at tend = %f ",
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
289 " was reached. This may happen because the @odeplot function",
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
290 " returned 'true' or the @event function returned",
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
291 " 'true'.\n"],
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
292 t(end), tspan(end));
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
293 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
294 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
295
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
296 ## Remove not-requested values of time and solution
20635
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
297 solution.t = t;
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20596
diff changeset
298 solution.x = x.';
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
299
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
300 endfunction
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
301