Mercurial > octave-nkf

comparison scripts/ode/private/starting_stepsize.m @ 20584:eb9e2d187ed2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Use Octave coding conventions in scripts/ode/private dir. * AbsRel_Norm.m, fuzzy_compare.m, hermite_quartic_interpolation.m, integrate_adaptive.m, integrate_const.m, integrate_n_steps.m, kahan.m, ode_struct_value_check.m, odepkg_event_handle.m, odepkg_structure_check.m, runge_kutta_45_dorpri.m, starting_stepsize.m: Wrap long lines to < 80 chars. Use double quotes rather than single quotes where possible. Use ';' at end of keywords "return;" and "break;" Use '##" for stand-alone comments and '#' for end-of-line comments. Use two spaces after period before starting new sentence. Use '!' instead of '~' for logical negation. Use specific form of end (endif, endfor, etc.). Don't use line continuation marker '...' unless necessary.
author Rik <rik@octave.org>
date Sun, 04 Oct 2015 22:18:54 -0700
parents 25623ef2ff4f
children
comparison
equal deleted inserted replaced
20583:d746695bf494 20584:eb9e2d187ed2
15 ## You should have received a copy of the GNU General Public License 15 ## You should have received a copy of the GNU General Public License
16 ## along with Octave; see the file COPYING. If not, see 16 ## along with Octave; see the file COPYING. If not, see
17 ## <http://www.gnu.org/licenses/>. 17 ## <http://www.gnu.org/licenses/>.
18 18
19 ## -*- texinfo -*- 19 ## -*- texinfo -*-
20 ## @deftypefn {Function File} {[@var{h}] =} starting_stepsize (@var{order}, @var{@@fun}, @var{t0}, @var{x0}) 20 ## @deftypefn {Function File} {@var{h} =} starting_stepsize (@var{order}, @var{@@fun}, @var{t0}, @var{x0})
21 ## 21 ##
22 ## This function file can be used to determine a good initial step for an ODE 22 ## This function file can be used to determine a good initial step for an ODE
23 ## solver of order @var{order}. The algorithm is that one described in [1]. 23 ## solver of order @var{order}. The algorithm is that one described in [1].
24 ## 24 ##
25 ## Second input argument, which is @var{@@fun}, is the function describing 25 ## Second input argument, which is @var{@@fun}, is the function describing
43 43
44 ## compute norm of the function evaluated at initial conditions 44 ## compute norm of the function evaluated at initial conditions
45 y = func (t0, x0); 45 y = func (t0, x0);
46 d1 = AbsRel_Norm (y, y, AbsTol, RelTol, normcontrol); 46 d1 = AbsRel_Norm (y, y, AbsTol, RelTol, normcontrol);
47 47
48 if (d0 < 1.e-5 || d1 < 1.e-5) 48 if (d0 < 1e-5 || d1 < 1e-5)
49 h0 = 1.e-6; 49 h0 = 1e-6;
50 else 50 else
51 h0 = .01 * (d0 / d1); 51 h0 = .01 * (d0 / d1);
52 endif 52 endif
53 53
54 ## compute one step of Explicit-Euler 54 ## compute one step of Explicit-Euler
57 ## approximate the derivative norm 57 ## approximate the derivative norm
58 d2 = (1 / h0) * ... 58 d2 = (1 / h0) * ...
59 AbsRel_Norm (func (t0+h0, x1) - y, 59 AbsRel_Norm (func (t0+h0, x1) - y,
60 func (t0+h0, x1) - y, AbsTol, RelTol, normcontrol); 60 func (t0+h0, x1) - y, AbsTol, RelTol, normcontrol);
61 61
62 if (max(d1, d2) <= 1.e-15) 62 if (max(d1, d2) <= 1e-15)
63 h1 = max (1.e-6, h0*1.e-3); 63 h1 = max (1e-6, h0*1e-3);
64 else 64 else
65 h1 = (1.e-2 / max (d1, d2)) ^(1 / (order+1)); 65 h1 = (1e-2 / max (d1, d2)) ^(1 / (order+1));
66 endif 66 endif
67 67
68 h = min (100*h0, h1); 68 h = min (100*h0, h1);
69 69
70 endfunction 70 endfunction
71 71
72 ## Local Variables: ***
73 ## mode: octave ***
74 ## End: ***