Mercurial > octave-nkf

annotate scripts/ode/private/integrate_adaptive.m @ 20596:87b557ee8e5d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
clean up and vectorize code for dense output in ode45 * scripts/ode/private/ode_rk_interpolate.m: new file * scripts/ode/private/ode_rk_interpolate.m(hermite_quartic_interpolation): move to internal function, use vectorization and broadcasting. * scripts/ode/private/hermite_quartic_interpolation.m: remove file * scripts/ode/module.mk: list added and removed files * scripts/ode/private/integrate_adaptive.m: use new interpolation code.
author Carlo de Falco <carlo.defalco@polimi.it>
date Tue, 06 Oct 2015 19:28:59 +0200
parents b7ac1e94266e
children a22d8a2eb0e5
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
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
65 solution = struct ();
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
66
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
67 ## first values for time and solution
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
68 t = tspan(1);
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
69 x = x0(:);
20575
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20570
diff changeset
70
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
71 ## get first initial timestep
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
72 dt = odeget (options, "InitialStep",
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
73 starting_stepsize (order, func, t, x, options.AbsTol,
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
74 options.RelTol, options.vnormcontrol),
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
75 "fast_not_empty");
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
76 vdirection = odeget (options, "vdirection", [], "fast");
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
77 if (sign (dt) != vdirection)
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
78 dt = -dt;
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
79 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
80 dt = vdirection * min (abs (dt), options.MaxStep);
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
81
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
82 ## Set parameters
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
83 k = length (tspan);
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
84 counter = 2;
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
85 comp = 0.0;
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
86 tk = tspan(1);
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
87 options.comp = comp;
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
88
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
89 ## 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
90 facmin = 0.8;
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
91 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
92 t_caught = false;
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
93
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
94
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
95 ## Initialize the OutputFcn
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
96 if (options.vhaveoutputfunction)
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
97 if (options.vhaveoutputselection)
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
98 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
99 else
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
100 solution.vretout = x;
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 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
103 "init", options.vfunarguments{:});
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
104 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
105
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
106 ## Initialize the EventFcn
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
107 if (options.vhaveeventfunction)
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
108 odepkg_event_handle (options.Events, t(end), x,
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
109 "init", options.vfunarguments{:});
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
110 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
111
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
112 solution.vcntloop = 2;
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
113 solution.vcntcycles = 1;
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
114 vcntiter = 0;
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
115 solution.vunhandledtermination = true;
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
116 solution.vcntsave = 2;
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
117
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
118 z = t;
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
119 u = x;
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
120
20575
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20570
diff changeset
121 k_vals = [];
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20570
diff changeset
122
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
123 while (counter <= k)
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
124 facmax = 1.5;
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
125
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
126 ## Compute integration step from t to t+dt
20575
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20570
diff changeset
127 if (isempty (k_vals))
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20570
diff changeset
128 [s, y, y_est, new_k_vals] = stepper (func, z(end), u(:,end),
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20570
diff changeset
129 dt, options);
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20570
diff changeset
130 else
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20570
diff changeset
131 [s, y, y_est, new_k_vals] = stepper (func, z(end), u(:,end),
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20570
diff changeset
132 dt, options, k_vals);
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20570
diff changeset
133 endif
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20570
diff changeset
134
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
135 if (options.vhavenonnegative)
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
136 x(options.NonNegative,end) = abs (x(options.NonNegative,end));
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
137 y(options.NonNegative,end) = abs (y(options.NonNegative,end));
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
138 y_est(options.NonNegative,end) = abs (y_est(options.NonNegative,end));
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
139 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
140
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
141 if (options.vhaveoutputfunction && options.vhaverefine)
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
142 vSaveVUForRefine = u(:,end);
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
143 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
144
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
145 err = AbsRel_Norm (y(:,end), u(:,end), options.AbsTol, options.RelTol,
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
146 options.vnormcontrol, y_est(:,end));
20575
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20570
diff changeset
147
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
148 ## Solution accepted only if the error is less or equal to 1.0
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
149 if (err <= 1)
20575
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20570
diff changeset
150
20596
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20586
diff changeset
151 [tk, options.comp] = kahan (tk, options.comp, dt);
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
152 s(end) = tk;
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
153
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
154 ## values on this interval for time and solution
20575
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20570
diff changeset
155 z = [z(end); s];
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20570
diff changeset
156 u = [u(:,end), y];
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20570
diff changeset
157 k_vals = new_k_vals;
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20570
diff changeset
158
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
159 ## if next tspan value is caught, update counter
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
160 if ((z(end) == tspan(counter))
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
161 || (abs (z(end) - tspan(counter)) /
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
162 (max (abs (z(end)), abs (tspan(counter)))) < 8*eps) )
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
163 counter++;
20575
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20570
diff changeset
164
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
165 ## if there is an element in time vector at which the solution is
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
166 ## required the program must compute this solution before going on with
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
167 ## next steps
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
168 elseif (vdirection * z(end) > vdirection * tspan(counter))
20575
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20570
diff changeset
169
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
170 ## initialize counter for the following cycle
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
171 i = 2;
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
172 while (i <= length (z))
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
173
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
174 ## if next tspan value is caught, update counter
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
175 if ((counter <= k)
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
176 && ((z(i) == tspan(counter))
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
177 || (abs (z(i) - tspan(counter)) /
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
178 (max (abs (z(i)), abs (tspan(counter)))) < 8*eps)) )
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
179 counter++;
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
180 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
181 ## else, loop until there are requested values inside this subinterval
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
182 while ((counter <= k)
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
183 && (vdirection * z(i) > vdirection * tspan(counter)))
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
184 ## choose interpolation scheme according to order of the solver
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
185
20596
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20586
diff changeset
186 u_interp = ...
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20586
diff changeset
187 ode_rk_interpolate (order, [z(i-1) z(i)], [u(:,i-1) u(:,i)],
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20586
diff changeset
188 tspan(counter), k_vals, dt,
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20586
diff changeset
189 options.vfunarguments{:});
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20586
diff changeset
190
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
191
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
192 ## add the interpolated value of the solution
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
193 u = [u(:,1:i-1), u_interp, u(:,i:end)];
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
194
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
195 ## add the time requested
20575
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20570
diff changeset
196 z = [z(1:i-1); tspan(counter); z(i:end)];
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
197
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
198 ## update counters
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
199 counter++;
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
200 i++;
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
201 endwhile
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
202
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
203 ## if new time requested is not out of this interval
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
204 if ((counter <= k)
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
205 && (vdirection * z(end) > vdirection * tspan(counter)))
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
206 ## update the counter
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
207 i++;
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
208 else
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
209 ## stop the cycle and go on with the next iteration
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
210 i = length (z) + 1;
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
211 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
212
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
213 endwhile
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
214 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
215
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
216 if (mod (solution.vcntloop-1, options.OutputSave) == 0)
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
217 x = [x,u(:,2:end)];
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
218 t = [t;z(2:end)];
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
219 solution.vcntsave = solution.vcntsave + 1;
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
220 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
221 solution.vcntloop = solution.vcntloop + 1;
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
222 vcntiter = 0;
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
223
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
224 ## Call plot only if a valid result has been found, therefore this
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
225 ## code fragment has moved here. Stop integration if plot function
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
226 ## returns false
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
227 if (options.vhaveoutputfunction)
20586
b7ac1e94266e maint: Further clean up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20584
diff changeset
228 for vcnt = 0:options.Refine # Approximation between told and t
b7ac1e94266e maint: Further clean up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20584
diff changeset
229 if (options.vhaverefine) # Do interpolation
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
230 vapproxtime = (vcnt + 1) / (options.Refine + 2);
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
231 vapproxvals = (1 - vapproxtime) * vSaveVUForRefine ...
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
232 + (vapproxtime) * y(:,end);
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
233 vapproxtime = s(end) + vapproxtime * dt;
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
234 else
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
235 vapproxvals = x(:,end);
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
236 vapproxtime = t(end);
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
237 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
238 if (options.vhaveoutputselection)
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
239 vapproxvals = vapproxvals(options.OutputSel);
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
240 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
241 vpltret = feval (options.OutputFcn, vapproxtime,
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
242 vapproxvals, [], options.vfunarguments{:});
20586
b7ac1e94266e maint: Further clean up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20584
diff changeset
243 if (vpltret) # Leave refinement loop
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
244 break;
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
245 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
246 endfor
20586
b7ac1e94266e maint: Further clean up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20584
diff changeset
247 if (vpltret) # Leave main loop
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
248 solution.vunhandledtermination = false;
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
249 break;
20568
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
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
252
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
253 ## Call event only if a valid result has been found, therefore this
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
254 ## code fragment has moved here. Stop integration if veventbreak is
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
255 ## true
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
256 if (options.vhaveeventfunction)
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
257 solution.vevent = odepkg_event_handle (options.Events, t(end),
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
258 x(:,end), [],
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
259 options.vfunarguments{:});
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
260 if (! isempty (solution.vevent{1})
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
261 && solution.vevent{1} == 1)
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
262 t(solution.vcntloop-1,:) = solution.vevent{3}(end,:);
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
263 x(:,solution.vcntloop-1) = solution.vevent{4}(end,:)';
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
264 solution.vunhandledtermination = false;
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
265 break;
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
266 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
267 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
268
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
269 else
20575
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20570
diff changeset
270
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
271 facmax = 1.0;
20575
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20570
diff changeset
272
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
273 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
274
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
275 ## Compute next timestep, formula taken from Hairer
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
276 err += eps; # adding an eps to avoid divisions by zero
20575
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20570
diff changeset
277 dt = dt * min (facmax,
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20570
diff changeset
278 max (facmin, fac * (1 / err)^(1 / (order + 1))));
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
279 dt = vdirection * min (abs (dt), options.MaxStep);
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
280
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
281 ## Update counters that count the number of iteration cycles
20586
b7ac1e94266e maint: Further clean up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20584
diff changeset
282 solution.vcntcycles += 1; # Needed for cost statistics
b7ac1e94266e maint: Further clean up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20584
diff changeset
283 vcntiter += 1; # Needed to find iteration problems
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
284
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
285 ## Stop solving because in the last 1000 steps no successful valid
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
286 ## value has been found
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
287 if (vcntiter >= 5000)
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
288 error (["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
289 " integration loop exited at time t = %f before endpoint at",
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
290 " tend = %f was reached. This happened because the iterative",
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
291 " integration loop does not find a valid solution at this time",
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
292 " stamp. Try to reduce the value of 'InitialStep' and/or",
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
293 " 'MaxStep' with the command 'odeset'.\n"],
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
294 s(end), tspan(end));
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
295 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
296
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
297 ## if this is the last iteration, save the length of last interval
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
298 if (counter > k)
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
299 j = length (z);
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
300 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
301 endwhile
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
302
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
303 ## Check if integration of the ode has been successful
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
304 if (vdirection * z(end) < vdirection * tspan(end))
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
305 if (solution.vunhandledtermination == true)
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
306 error ("OdePkg:InvalidArgument",
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
307 ["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
308 " 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
309 " before endpoint at tend = %f was reached. This may",
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
310 " happen if the stepsize grows smaller than defined in",
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
311 " vminstepsize. Try to reduce the value of 'InitialStep'",
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
312 " and/or 'MaxStep' with the command 'odeset'.\n"],
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
313 z(end), tspan(end));
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
314 else
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
315 warning ("OdePkg:InvalidArgument",
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
316 ["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
317 " 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
318 " 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
319 " 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
320 " 'true'.\n"],
20568
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
321 z(end), tspan(end));
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
322 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
323 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
324
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
325 ## Compute how many values are out of time inerval
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
326 d = vdirection * t((end-(j-1)):end) > vdirection * tspan(end)*ones (j, 1);
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
327 f = sum (d);
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
328
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
329 ## Remove not-requested values of time and solution
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
330 solution.t = t(1:end-f);
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
331 solution.x = x(:,1:end-f)';
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
332
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
333 endfunction
20584
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20580
diff changeset
334