Mercurial > octave
annotate scripts/ode/private/integrate_adaptive.m @ 20621:b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
* AbsRel_Norm.m: Use retval as return variable. Use sumsq() rather than explicit
squaring of vector and sum(). Combine multiple lines where possible.
* integrate_adaptive.m: Rewrite docstring. Only call starting_stepsize() if
InitialStep option is empty.
* integrate_const.m: Rewrite docstring. Remove useless commented out code.
Combine multiple lines where possible.
* integrate_n_steps.m: Rewrite docstring. Remove useless commented out code.
Combine multiple lines where possible.
* kahan.m: Remove excessive 4-space indentation, use 2-space indentation.
* ode_rk_interpolate.m: Use parentheses around condition for switch stmt.
Combine multiple lines where possible.
* ode_struct_value_check.m: Remove comma from Copyright statement that Octave
doesn't use.
* odepkg_event_handle.m: Remove comma from Copyright statement that Octave
doesn't use.
* odepkg_structure_check.m: Remove comma from Copyright statement that Octave
doesn't use.
* runge_kutta_45_dorpri.m: Remove comma from Copyright statement that Octave
doesn't use. Improve docstring. Match variable names in documentation to
those in code.
* starting_stepsize.m: Rewrite docstring. Use spaces between function name
and opening parenthesis.
author | Rik <rik@octave.org> |
---|---|
date | Wed, 14 Oct 2015 10:35:53 -0700 |
parents | a260a6acb70f |
children | 80e630b37ba1 |
rev | line source |
---|---|
20621
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
1 ## Copyright (C) 2015 Carlo de Falco |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
2 ## Copyright (C) 2013 Roberto Porcu' <roberto.porcu@polimi.it> |
20536
6256f6e366ac
Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20533
diff
changeset
|
3 ## |
6256f6e366ac
Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20533
diff
changeset
|
4 ## This file is part of Octave. |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
5 ## |
20536
6256f6e366ac
Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20533
diff
changeset
|
6 ## Octave is free software; you can redistribute it and/or modify it |
6256f6e366ac
Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20533
diff
changeset
|
7 ## under the terms of the GNU General Public License as published by |
6256f6e366ac
Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20533
diff
changeset
|
8 ## the Free Software Foundation; either version 3 of the License, or (at |
6256f6e366ac
Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20533
diff
changeset
|
9 ## your option) any later version. |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
10 ## |
20536
6256f6e366ac
Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20533
diff
changeset
|
11 ## Octave is distributed in the hope that it will be useful, but |
6256f6e366ac
Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20533
diff
changeset
|
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
6256f6e366ac
Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20533
diff
changeset
|
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
6256f6e366ac
Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20533
diff
changeset
|
14 ## General Public License for more details. |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
15 ## |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
16 ## You should have received a copy of the GNU General Public License |
20536
6256f6e366ac
Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20533
diff
changeset
|
17 ## along with Octave; see the file COPYING. If not, see |
6256f6e366ac
Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20533
diff
changeset
|
18 ## <http://www.gnu.org/licenses/>. |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
19 |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
20 ## -*- texinfo -*- |
20621
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
21 ## @deftypefn {Function File} {@var{solution} =} integrate_adaptive (@var{@@stepper}, @var{order}, @var{@@func}, @var{tspan}, @var{x0}, @var{options}) |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
22 ## |
20548
25623ef2ff4f
doc: Rewrite docstrings for ode* family of functions.
Rik <rik@octave.org>
parents:
20543
diff
changeset
|
23 ## This function file can be called by an ODE solver function in order to |
20621
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
24 ## integrate the set of ODEs on the interval @var{[t0, t1]} with an |
20533
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 ## |
20621
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
27 ## The function returns a structure @var{solution} with two fieldss: @var{t} |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
28 ## and @var{y}. @var{t} is a column vector and contains the time stamps. |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
29 ## @var{y} is a matrix in which each column refers to a different unknown |
20621
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
30 ## of the problem and the row number is the same as the @var{t} row number. |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
31 ## Thus, each row of the matrix @var{y} contains the values of all unknowns at |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
32 ## the time value contained in the corresponding row in @var{t}. |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
33 ## |
20621
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
34 ## The first input argument must be a function handle or inline function |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
35 ## representing the stepper, i.e., the function responsible for step-by-step |
20548
25623ef2ff4f
doc: Rewrite docstrings for ode* family of functions.
Rik <rik@octave.org>
parents:
20543
diff
changeset
|
36 ## integration. This function discriminates one method from the others. |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
37 ## |
20548
25623ef2ff4f
doc: Rewrite docstrings for ode* family of functions.
Rik <rik@octave.org>
parents:
20543
diff
changeset
|
38 ## The second input argument is the order of the stepper. It is needed |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
39 ## to compute the adaptive timesteps. |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
40 ## |
20621
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
41 ## The third input argument is a function handle or inline function that |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
42 ## defines the ODE: |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
43 ## |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
44 ## @ifhtml |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
45 ## @example |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
46 ## @math{y' = f(t,y)} |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
47 ## @end example |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
48 ## @end ifhtml |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
49 ## @ifnothtml |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
50 ## @math{y' = f(t,y)}. |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
51 ## @end ifnothtml |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
52 ## |
20621
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
53 ## The fourth input argument is the time vector which defines the integration |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
54 ## interval, i.e, @var{[tspan(1), tspan(end)]} and all intermediate elements |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
55 ## are taken as times at which the solution is required. |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
56 ## |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
57 ## 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
|
58 ## 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
|
59 ## |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
60 ## @end deftypefn |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
61 ## |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
62 ## @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
|
63 |
20621
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
64 function solution = integrate_adaptive (stepper, order, func, tspan, x0, |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
65 options) |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
66 |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
67 fixed_times = numel (tspan) > 2; |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
68 |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
69 t_new = t_old = t = tspan(1); |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
70 x_new = x_old = x = x0(:); |
20602
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20601
diff
changeset
|
71 |
20621
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
72 ## Get first initial timestep |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
73 dt = odeget (options, "InitialStep", [], "fast"); |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
74 if (isempty (dt)) |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
75 dt = starting_stepsize (order, func, t, x, options.AbsTol, options.RelTol, |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
76 options.vnormcontrol); |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
77 endif |
20602
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20601
diff
changeset
|
78 |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
79 dir = odeget (options, "vdirection", [], "fast"); |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
80 dt = dir * min (abs (dt), options.MaxStep); |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
81 |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
82 options.comp = 0.0; |
20602
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20601
diff
changeset
|
83 |
20552
eb9e2d187ed2
maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents:
20548
diff
changeset
|
84 ## Factor multiplying the stepsize guess |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
85 facmin = 0.8; |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
86 facmax = 1.5; |
20552
eb9e2d187ed2
maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents:
20548
diff
changeset
|
87 fac = 0.38^(1/(order+1)); # formula taken from Hairer |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
88 |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
89 ## Initialize the OutputFcn |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
90 if (options.vhaveoutputfunction) |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
91 if (options.vhaveoutputselection) |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
92 solution.vretout = x(options.OutputSel,end); |
20602
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20601
diff
changeset
|
93 else |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
94 solution.vretout = x; |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
95 endif |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
96 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
|
97 "init", options.vfunarguments{:}); |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
98 endif |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
99 |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
100 ## Initialize the EventFcn |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
101 if (options.vhaveeventfunction) |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
102 odepkg_event_handle (options.Events, tspan(end), x, |
20533
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 |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
106 if (options.vhavenonnegative) |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
107 nn = options.NonNegative; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
108 endif |
20602
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20601
diff
changeset
|
109 |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
110 solution.vcntloop = 2; |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
111 solution.vcntcycles = 1; |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
112 solution.vcntsave = 2; |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
113 solution.vunhandledtermination = true; |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
114 ireject = 0; |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
115 |
20602
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20601
diff
changeset
|
116 k_vals = []; |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
117 iout = istep = 1; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
118 while (dir * t_old < dir * tspan(end)) |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
119 |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
120 ## Compute integration step from t_old to t_new = t_old + dt |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
121 [t_new, options.comp] = kahan (t_old, options.comp, dt); |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
122 [t_new, x_new, x_est, new_k_vals] = ... |
20621
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
123 stepper (func, t_old, x_old, dt, options, k_vals, t_new); |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
124 |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
125 solution.vcntcycles++; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
126 |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
127 if (options.vhavenonnegative) |
20604
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20602
diff
changeset
|
128 x_new(nn, end) = abs (x_new(nn, end)); |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20602
diff
changeset
|
129 x_est(nn, end) = abs (x_est(nn, end)); |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
130 endif |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
131 |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
132 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:
20565
diff
changeset
|
133 options.vnormcontrol, x_est); |
20602
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20601
diff
changeset
|
134 |
20621
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
135 ## Accept solution only if err <= 1.0 |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
136 if (err <= 1) |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
137 |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
138 solution.vcntloop++; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
139 ireject = 0; |
20602
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20601
diff
changeset
|
140 |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
141 ## if output time steps are fixed |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
142 if (fixed_times) |
20543
3339c9bdfe6a
Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20536
diff
changeset
|
143 |
20604
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20602
diff
changeset
|
144 t_caught = find ((dir * tspan(iout:end) > dir * t_old) |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20602
diff
changeset
|
145 & (dir * tspan(iout:end) <= dir * t_new)); |
20601
43822bda4f65
fix indexing bug introduced with a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20600
diff
changeset
|
146 t_caught = t_caught + iout - 1; |
20602
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20601
diff
changeset
|
147 |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
148 if (! isempty (t_caught)) |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
149 t(t_caught) = tspan(t_caught); |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
150 iout = max (t_caught); |
20601
43822bda4f65
fix indexing bug introduced with a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20600
diff
changeset
|
151 x(:, t_caught) = ... |
20621
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
152 ode_rk_interpolate (order, [t_old t_new], [x_old x_new], |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
153 tspan(t_caught), new_k_vals, dt, |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
154 options.vfunarguments{:}); |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
155 |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
156 istep++; |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
157 |
20621
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
158 ## Call Events function only if a valid result has been found. |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
159 ## Stop integration if veventbreak is true. |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
160 if (options.vhaveeventfunction) |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
161 break_loop = false; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
162 for idenseout = 1:numel (t_caught) |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
163 id = t_caught(idenseout); |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
164 td = t(id); |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
165 solution.vevent = ... |
20621
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
166 odepkg_event_handle (options.Events, t(id), x(:, id), [], |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
167 options.vfunarguments{:}); |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
168 if (! isempty (solution.vevent{1}) && solution.vevent{1} == 1) |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
169 t(id) = solution.vevent{3}(end); |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
170 t = t(1:id); |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
171 x(:, id) = solution.vevent{4}(end, :).'; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
172 x = x(:,1:id); |
20602
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20601
diff
changeset
|
173 solution.vunhandledtermination = false; |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
174 break_loop = true; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
175 break; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
176 endif |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
177 endfor |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
178 if (break_loop) |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
179 break; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
180 endif |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
181 endif |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
182 |
20621
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
183 ## Call OutputFcn only if a valid result has been found. |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
184 ## Stop integration if function returns false. |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
185 if (options.vhaveoutputfunction) |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
186 vcnt = options.Refine + 1; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
187 vapproxtime = linspace (t_old, t_new, vcnt); |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
188 vapproxvals = interp1 ([t_old, t(t_caught), t_new], |
20604
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20602
diff
changeset
|
189 [x_old, x(:, t_caught), x_new] .', |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20602
diff
changeset
|
190 vapproxtime, 'linear') .'; |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
191 if (options.vhaveoutputselection) |
20604
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20602
diff
changeset
|
192 vapproxvals = vapproxvals(options.OutputSel, :); |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
193 endif |
20604
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20602
diff
changeset
|
194 for ii = 1:numel (vapproxtime) |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20602
diff
changeset
|
195 vpltret = feval (options.OutputFcn, vapproxtime(ii), |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20602
diff
changeset
|
196 vapproxvals(:, ii), [], |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20602
diff
changeset
|
197 options.vfunarguments{:}); |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20602
diff
changeset
|
198 endfor |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
199 if (vpltret) # Leave main loop |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
200 solution.vunhandledtermination = false; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
201 break; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
202 endif |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
203 endif |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
204 |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
205 endif |
20602
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20601
diff
changeset
|
206 |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
207 else |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
208 |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
209 t(++istep) = t_new; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
210 x(:, istep) = x_new; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
211 iout = istep; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
212 |
20621
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
213 ## Call Events function only if a valid result has been found. |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
214 ## Stop integration if veventbreak is true. |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
215 if (options.vhaveeventfunction) |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
216 solution.vevent = ... |
20621
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
217 odepkg_event_handle (options.Events, t(istep), x(:, istep), [], |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
218 options.vfunarguments{:}); |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
219 if (! isempty (solution.vevent{1}) && solution.vevent{1} == 1) |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
220 t(istep) = solution.vevent{3}(end); |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
221 x(:, istep) = solution.vevent{4}(end, :).'; |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
222 solution.vunhandledtermination = false; |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
223 break; |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
224 endif |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
225 endif |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
226 |
20621
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
227 ## Call OutputFcn only if a valid result has been found. |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
228 ## Stop integration if function returns false. |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
229 if (options.vhaveoutputfunction) |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
230 vcnt = options.Refine + 1; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
231 vapproxtime = linspace (t_old, t_new, vcnt); |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
232 vapproxvals = interp1 ([t_old, t_new], |
20604
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20602
diff
changeset
|
233 [x_old, x_new] .', |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20602
diff
changeset
|
234 vapproxtime, 'linear') .'; |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
235 if (options.vhaveoutputselection) |
20604
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20602
diff
changeset
|
236 vapproxvals = vapproxvals(options.OutputSel, :); |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
237 endif |
20604
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20602
diff
changeset
|
238 for ii = 1:numel (vapproxtime) |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20602
diff
changeset
|
239 vpltret = feval (options.OutputFcn, vapproxtime(ii), |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20602
diff
changeset
|
240 vapproxvals(:, ii), [], options.vfunarguments{:}); |
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20602
diff
changeset
|
241 endfor |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
242 if (vpltret) # Leave main loop |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
243 solution.vunhandledtermination = false; |
20552
eb9e2d187ed2
maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents:
20548
diff
changeset
|
244 break; |
20533
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 endif |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
247 |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
248 endif |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
249 |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
250 ## move to next time-step |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
251 t_old = t_new; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
252 x_old = x_new; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
253 k_vals = new_k_vals; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
254 |
20621
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
255 solution.vcntloop += 1; |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
256 vcntiter = 0; |
20602
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20601
diff
changeset
|
257 |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
258 else |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
259 |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
260 ireject++; |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
261 |
20621
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
262 ## Stop solving because, in the last 5,000 steps, no successful valid |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
263 ## value has been found |
20621
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
264 if (ireject >= 5_000) |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
265 error (["integrate_adaptive: Solving was not successful. ", ... |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
266 " The iterative integration loop exited at time", ... |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
267 " t = %f before the endpoint at tend = %f was reached. ", ... |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
268 " This happened because the iterative integration loop", ... |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
269 " did not find a valid solution at this time stamp. ", ... |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
270 " Try to reduce the value of 'InitialStep' and/or", ... |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
271 " 'MaxStep' with the command 'odeset'.\n"], |
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
272 t_old, tspan(end)); |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
273 endif |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
274 |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
275 endif |
20602
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20601
diff
changeset
|
276 |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
277 ## Compute next timestep, formula taken from Hairer |
20621
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
278 err += eps; # avoid divisions by zero |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
279 dt *= min (facmax, max (facmin, fac * (1 / err)^(1 / (order + 1)))); |
20602
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20601
diff
changeset
|
280 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:
20601
diff
changeset
|
281 |
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20601
diff
changeset
|
282 ## 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:
20601
diff
changeset
|
283 dt = dir * min (abs (dt), abs (tspan(end) - t_old)); |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
284 |
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
285 endwhile |
20602
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20601
diff
changeset
|
286 |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
287 ## Check if integration of the ode has been successful |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
288 if (dir * t(end) < dir * tspan(end)) |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
289 if (solution.vunhandledtermination == true) |
20601
43822bda4f65
fix indexing bug introduced with a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20600
diff
changeset
|
290 error ("integrate_adaptive:unexpected_termination", |
20621
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
291 [" Solving was not successful. ", ... |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
292 " The iterative integration loop exited at time", ... |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
293 " t = %f before the endpoint at tend = %f was reached. ", ... |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
294 " This may happen if the stepsize becomes too small. ", ... |
20604
a260a6acb70f
fix test failures introduced by a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20602
diff
changeset
|
295 " Try to reduce the value of 'InitialStep'", ... |
20552
eb9e2d187ed2
maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents:
20548
diff
changeset
|
296 " and/or 'MaxStep' with the command 'odeset'.\n"], |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
297 t(end), tspan(end)); |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
298 else |
20601
43822bda4f65
fix indexing bug introduced with a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20600
diff
changeset
|
299 warning ("integrate_adaptive:unexpected_termination", |
20621
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
300 ["Solver was stopped by a call of 'break'", ... |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
301 " in the main iteration loop at time ", ... |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
302 " t = %f before the endpoint at tend = %f was reached. ", ... |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
303 " This may happen because the @odeplot function ", ... |
b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents:
20604
diff
changeset
|
304 " returned 'true' or the @event function returned 'true'.\n"], |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
305 t(end), tspan(end)); |
20533
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 |
20601
43822bda4f65
fix indexing bug introduced with a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20600
diff
changeset
|
309 ## Set up return structure |
43822bda4f65
fix indexing bug introduced with a22d8a2eb0e5
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20600
diff
changeset
|
310 solution.t = t(:); |
20600
a22d8a2eb0e5
fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20565
diff
changeset
|
311 solution.x = x.'; |
20602
756b052037fb
avoid stepping beyond end of thspan in ode solvers
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20601
diff
changeset
|
312 |
20533
fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff
changeset
|
313 endfunction |
20552
eb9e2d187ed2
maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents:
20548
diff
changeset
|
314 |