Mercurial > octave

annotate scripts/ode/private/starting_stepsize.m @ 30564:796f54d4ddbf stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
update Octave Project Developers copyright for the new year In files that have the "Octave Project Developers" copyright notice, update for 2021. In all .txi and .texi files except gpl.txi and gpl.texi in the doc/liboctave and doc/interpreter directories, change the copyright to "Octave Project Developers", the same as used for other source files. Update copyright notices for 2022 (not done since 2019). For gpl.txi and gpl.texi, change the copyright notice to be "Free Software Foundation, Inc." and leave the date at 2007 only because this file only contains the text of the GPL, not anything created by the Octave Project Developers. Add Paul Thomas to contributors.in.
author John W. Eaton <jwe@octave.org>
date Tue, 28 Dec 2021 18:22:40 -0500
parents 0a5b15007766
children e1788b1a315f
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
27923
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
1 ########################################################################
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
2 ##
30564
796f54d4ddbf update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents: 29358
diff changeset
3 ## Copyright (C) 2013-2022 The Octave Project Developers
27918
b442ec6dda5c use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents: 26376
diff changeset
4 ##
27923
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
5 ## See the file COPYRIGHT.md in the top-level directory of this
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
6 ## distribution or <https://octave.org/copyright/>.
20536
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
7 ##
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
8 ## 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
9 ##
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23220
diff changeset
10 ## Octave is free software: you can redistribute it and/or modify it
20536
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
11 ## under the terms of the GNU General Public License as published by
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23220
diff changeset
12 ## the Free Software Foundation, either version 3 of the License, or
22755
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22626
diff changeset
13 ## (at 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
14 ##
20536
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
15 ## 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
16 ## WITHOUT ANY WARRANTY; without even the implied warranty of
22755
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22626
diff changeset
17 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22626
diff changeset
18 ## GNU 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
19 ##
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
20 ## 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
21 ## along with Octave; see the file COPYING. If not, see
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23220
diff changeset
22 ## <https://www.gnu.org/licenses/>.
27923
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
23 ##
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
24 ########################################################################
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
25
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
26 ## -*- texinfo -*-
22768
a75a0900e216 maint: Match documentation names and variable names in more m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
27 ## @deftypefn {} {@var{h} =} starting_stepsize (@var{order}, @var{func}, @var{t0}, @var{x0}, @var{AbsTol}, @var{RelTol}, @var{normcontrol}, @var{args})
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
28 ##
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20552
diff changeset
29 ## Determine a good initial timestep for an ODE solver of order @var{order}
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20552
diff changeset
30 ## using the algorithm described in reference [1].
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
31 ##
22625
081a201b77c7 Clean up ode options implementation to follow Octave coding standards.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 22605
diff changeset
32 ## The input argument @var{func}, is the function describing the differential
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20552
diff changeset
33 ## equations, @var{t0} is the initial time, and @var{x0} is the initial
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20552
diff changeset
34 ## condition. @var{AbsTol} and @var{RelTol} are the absolute and relative
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20552
diff changeset
35 ## tolerance on the ODE integration taken from an ode options structure.
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
36 ##
27984
b09432b20a84 maint: Remove special cases of old version control keywords in code base.
Rik <rik@octave.org>
parents: 27923
diff changeset
37 ## Reference:
b09432b20a84 maint: Remove special cases of old version control keywords in code base.
Rik <rik@octave.org>
parents: 27923
diff changeset
38 ## E. Hairer, S.P. Norsett and G. Wanner,
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20552
diff changeset
39 ## @cite{Solving Ordinary Differential Equations I: Nonstiff Problems},
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20552
diff changeset
40 ## Springer.
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
41 ## @end deftypefn
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
42 ##
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
43 ## @seealso{odepkg}
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
44
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
45 function h = starting_stepsize (order, func, t0, x0,
22605
177e0c71bcc0 make sure the additional function arguments are always passed to odefun.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 22596
diff changeset
46 AbsTol, RelTol, normcontrol,
177e0c71bcc0 make sure the additional function arguments are always passed to odefun.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 22596
diff changeset
47 args = {})
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
48
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
49 ## compute norm of initial conditions
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22625
diff changeset
50 d0 = AbsRel_norm (x0, x0, AbsTol, RelTol, normcontrol);
20533
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 ## compute norm of the function evaluated at initial conditions
22596
8d3a2d1af389 fix missing function arguments in call to starting_stepsize
jcorno <jacopo.corno@gmail.com>
parents: 22323
diff changeset
53 y = func (t0, x0, args{:});
20847
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
54 if (iscell (y))
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
55 y = y{1};
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
56 endif
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22625
diff changeset
57 d1 = AbsRel_norm (y, y, AbsTol, RelTol, normcontrol);
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
58
20552
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20548
diff changeset
59 if (d0 < 1e-5 || d1 < 1e-5)
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20548
diff changeset
60 h0 = 1e-6;
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
61 else
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
62 h0 = .01 * (d0 / d1);
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
63 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
64
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
65 ## compute one step of Explicit-Euler
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
66 x1 = x0 + h0 * y;
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
67
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
68 ## approximate the derivative norm
22596
8d3a2d1af389 fix missing function arguments in call to starting_stepsize
jcorno <jacopo.corno@gmail.com>
parents: 22323
diff changeset
69 yh = func (t0+h0, x1, args{:});
20847
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
70 if (iscell (yh))
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
71 yh = yh{1};
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
72 endif
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
73 d2 = (1 / h0) * ...
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22625
diff changeset
74 AbsRel_norm (yh - y, yh - y, AbsTol, RelTol, normcontrol);
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
75
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20552
diff changeset
76 if (max (d1, d2) <= 1e-15)
22625
081a201b77c7 Clean up ode options implementation to follow Octave coding standards.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 22605
diff changeset
77 h1 = max (1e-6, h0 * 1e-3);
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
78 else
20552
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20548
diff changeset
79 h1 = (1e-2 / max (d1, d2)) ^(1 / (order+1));
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
80 endif
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
81
22625
081a201b77c7 Clean up ode options implementation to follow Octave coding standards.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 22605
diff changeset
82 h = min (100 * h0, h1);
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
83
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
84 endfunction