Mercurial > octave-nkf
annotate scripts/ode/private/integrate_adaptive.m @ 20654:b65888ec820e draft default tip gccjit
dmalcom gcc jit import
author | Stefan Mahr <dac922@gmx.de> |
---|---|
date | Fri, 27 Feb 2015 16:59:36 +0100 |
parents | a260a6acb70f |
children |
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(:); |
20637
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20636
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"); |
20637
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20636
diff
changeset
|
74 |
20635
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; |
20637
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20636
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); |
20637
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20636
diff
changeset
|
89 else |
20568
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 |
20637
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20636
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 |
20637
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20636
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) |
20639
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
124 x_new(nn, end) = abs (x_new(nn, end)); |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
125 x_est(nn, end) = abs (x_est(nn, end)); |
20568
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
126 endif |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
127 |
20635
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
128 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
|
129 options.vnormcontrol, x_est); |
20637
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20636
diff
changeset
|
130 |
20635
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
131 ## 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
|
132 if (err <= 1) |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
133 |
20635
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
134 solution.vcntloop++; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
135 ireject = 0; |
20637
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20636
diff
changeset
|
136 |
20635
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
137 ## if output time steps are fixed |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
138 if (fixed_times) |
20575
3339c9bdfe6a
Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20570
diff
changeset
|
139 |
20639
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
140 t_caught = find ((dir * tspan(iout:end) > dir * t_old) |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
141 & (dir * tspan(iout:end) <= dir * t_new)); |
20636
43822bda4f65
fix indexing bug introduced with a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20635
diff
changeset
|
142 t_caught = t_caught + iout - 1; |
20637
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20636
diff
changeset
|
143 |
20635
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
144 if (! isempty (t_caught)) |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
145 t(t_caught) = tspan(t_caught); |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
146 iout = max (t_caught); |
20636
43822bda4f65
fix indexing bug introduced with a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20635
diff
changeset
|
147 x(:, t_caught) = ... |
43822bda4f65
fix indexing bug introduced with a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20635
diff
changeset
|
148 ode_rk_interpolate (order, [t_old t_new], [x_old x_new], |
43822bda4f65
fix indexing bug introduced with a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20635
diff
changeset
|
149 tspan(t_caught), new_k_vals, dt, |
43822bda4f65
fix indexing bug introduced with a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20635
diff
changeset
|
150 options.vfunarguments{:}); |
20635
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
151 |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
152 istep++; |
20568
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
153 |
20635
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
154 if (options.vhaveeventfunction) |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
155 ## 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
|
156 ## Stop integration if veventbreak is true |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
157 break_loop = false; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
158 for idenseout = 1:numel (t_caught) |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
159 id = t_caught(idenseout); |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
160 td = t(id); |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
161 solution.vevent = ... |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
162 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
|
163 options.vfunarguments{:}); |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
164 if (! isempty (solution.vevent{1}) |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
165 && solution.vevent{1} == 1) |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
166 t(id) = solution.vevent{3}(end); |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
167 t = t(1:id); |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
168 x(:, id) = solution.vevent{4}(end, :).'; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
169 x = x(:,1:id); |
20637
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20636
diff
changeset
|
170 solution.vunhandledtermination = false; |
20635
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
171 break_loop = true; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
172 break; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
173 endif |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
174 endfor |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
175 if (break_loop) |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
176 break; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
177 endif |
20568
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
178 endif |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
179 |
20635
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
180 ## 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
|
181 ## returns true. |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
182 if (options.vhaveoutputfunction) |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
183 vcnt = options.Refine + 1; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
184 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
|
185 vapproxvals = interp1 ([t_old, t(t_caught), t_new], |
20639
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
186 [x_old, x(:, t_caught), x_new] .', |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
187 vapproxtime, 'linear') .'; |
20635
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
188 if (options.vhaveoutputselection) |
20639
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
189 vapproxvals = vapproxvals(options.OutputSel, :); |
20635
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
190 endif |
20639
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
191 for ii = 1:numel (vapproxtime) |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
192 vpltret = feval (options.OutputFcn, vapproxtime(ii), |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
193 vapproxvals(:, ii), [], |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
194 options.vfunarguments{:}); |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
195 endfor |
20635
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
196 if (vpltret) # Leave main loop |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
197 solution.vunhandledtermination = false; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
198 break; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
199 endif |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
200 endif |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
201 |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
202 endif |
20637
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20636
diff
changeset
|
203 |
20635
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
204 else |
20568
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
205 |
20635
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
206 t(++istep) = t_new; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
207 x(:, istep) = x_new; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
208 iout = istep; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
209 |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
210 ## Call event handler on new timestep. |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
211 ## Stop integration if veventbreak is true |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
212 if (options.vhaveeventfunction) |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
213 solution.vevent = ... |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
214 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
|
215 options.vfunarguments{:}); |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
216 if (! isempty (solution.vevent{1}) |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
217 && solution.vevent{1} == 1) |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
218 t(istep) = solution.vevent{3}(end); |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
219 x(:, istep) = solution.vevent{4}(end, :).'; |
20637
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20636
diff
changeset
|
220 solution.vunhandledtermination = false; |
20635
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
221 break; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
222 endif |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
223 endif |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
224 |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
225 ## 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
|
226 ## returns true. |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
227 if (options.vhaveoutputfunction) |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
228 vcnt = options.Refine + 1; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
229 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
|
230 vapproxvals = interp1 ([t_old, t_new], |
20639
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
231 [x_old, x_new] .', |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
232 vapproxtime, 'linear') .'; |
20568
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
233 if (options.vhaveoutputselection) |
20639
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
234 vapproxvals = vapproxvals(options.OutputSel, :); |
20568
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
235 endif |
20639
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
236 for ii = 1:numel (vapproxtime) |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
237 vpltret = feval (options.OutputFcn, vapproxtime(ii), |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
238 vapproxvals(:, ii), [], options.vfunarguments{:}); |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
239 endfor |
20635
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
240 if (vpltret) # Leave main loop |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
241 solution.vunhandledtermination = false; |
20584
eb9e2d187ed2
maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents:
20580
diff
changeset
|
242 break; |
20568
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
243 endif |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
244 endif |
20635
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
245 |
20568
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
246 endif |
20635
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 ## move to next time-step |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
249 t_old = t_new; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
250 x_old = x_new; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
251 k_vals = new_k_vals; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
252 |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
253 solution.vcntloop = solution.vcntloop + 1; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
254 vcntiter = 0; |
20637
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20636
diff
changeset
|
255 |
20635
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
256 else |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
257 |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
258 ireject++; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
259 |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
260 ## 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
|
261 ## value has been found |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
262 if (ireject >= 5000) |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
263 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
|
264 " 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
|
265 " 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
|
266 " 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
|
267 " 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
|
268 " 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
|
269 " 'MaxStep' with the command 'odeset'.\n"], |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
270 t_old, tspan(end)); |
20568
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
271 endif |
20635
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
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 |
20637
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20636
diff
changeset
|
274 |
20568
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 |
20635
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
276 err += eps; # avoid divisions by zero |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
277 dt *= min (facmax, max (facmin, fac * (1 / err)^(1 / (order + 1)))); |
20637
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20636
diff
changeset
|
278 dt = dir * min (abs (dt), options.MaxStep); |
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20636
diff
changeset
|
279 |
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20636
diff
changeset
|
280 ## make sure we don't go past tpan(end) |
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20636
diff
changeset
|
281 dt = dir * min (abs (dt), abs (tspan(end) - t_old)); |
20568
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
282 |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
283 endwhile |
20637
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20636
diff
changeset
|
284 |
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20636
diff
changeset
|
285 |
20568
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
286 ## 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
|
287 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
|
288 if (solution.vunhandledtermination == true) |
20636
43822bda4f65
fix indexing bug introduced with a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20635
diff
changeset
|
289 error ("integrate_adaptive:unexpected_termination", |
20639
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
290 [" Solving has not been successful. The iterative", ... |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
291 " integration loop exited at time t = %f ", ... |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
292 " before endpoint at tend = %f was reached. This may", ... |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
293 " happen if the stepsize grows too small. ", ... |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
294 " 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
|
295 " 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
|
296 t(end), tspan(end)); |
20568
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
297 else |
20636
43822bda4f65
fix indexing bug introduced with a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20635
diff
changeset
|
298 warning ("integrate_adaptive:unexpected_termination", |
20639
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
299 ["Solver has been stopped by a call of 'break' ", ... |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
300 " in the main iteration loop at time t = %f before", ... |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
301 " endpoint at tend = %f was reached. This may", ... |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
302 " happen because the @odeplot function", ... |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20637
diff
changeset
|
303 " returned 'true' or the @event function returned", ... |
20584
eb9e2d187ed2
maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents:
20580
diff
changeset
|
304 " 'true'.\n"], |
20635
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
305 t(end), tspan(end)); |
20568
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
306 endif |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
307 endif |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
308 |
20636
43822bda4f65
fix indexing bug introduced with a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20635
diff
changeset
|
309 ## Set up return structure |
43822bda4f65
fix indexing bug introduced with a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20635
diff
changeset
|
310 solution.t = t(:); |
20635
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
311 solution.x = x.'; |
20637
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20636
diff
changeset
|
312 |
20568
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
313 endfunction |
20584
eb9e2d187ed2
maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents:
20580
diff
changeset
|
314 |