Mercurial > octave

annotate scripts/ode/decic.m @ 31225:3eab70385569

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
sparse-xpow.cc: Use faster multiplication technique, this time for complex
author Arun Giridhar <arungiridhar@gmail.com>
date Sun, 11 Sep 2022 13:53:38 -0400
parents e1788b1a315f
children 597f3ee61a48
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: 29359
diff changeset
3 ## Copyright (C) 2016-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/>.
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
7 ##
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
8 ## This file is part of Octave.
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
9 ##
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 24007
diff changeset
10 ## Octave is free software: you can redistribute it and/or modify it
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
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: 24007
diff changeset
12 ## the Free Software Foundation, either version 3 of the License, or
24007
e8a74d95b4f3 maint: Use same format for Copyright statement throught code base.
Rik <rik@octave.org>
parents: 23086
diff changeset
13 ## (at your option) any later version.
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
14 ##
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
15 ## Octave is distributed in the hope that it will be useful, but
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
16 ## WITHOUT ANY WARRANTY; without even the implied warranty of
24007
e8a74d95b4f3 maint: Use same format for Copyright statement throught code base.
Rik <rik@octave.org>
parents: 23086
diff changeset
17 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
e8a74d95b4f3 maint: Use same format for Copyright statement throught code base.
Rik <rik@octave.org>
parents: 23086
diff changeset
18 ## GNU General Public License for more details.
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
19 ##
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
20 ## You should have received a copy of the GNU General Public License
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
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: 24007
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 ########################################################################
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
25
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
26 ## -*- texinfo -*-
30893
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files.
Rik <rik@octave.org>
parents: 30564
diff changeset
27 ## @deftypefn {} {[@var{y0_new}, @var{yp0_new}] =} decic (@var{fcn}, @var{t0}, @var{y0}, @var{fixed_y0}, @var{yp0}, @var{fixed_yp0})
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files.
Rik <rik@octave.org>
parents: 30564
diff changeset
28 ## @deftypefnx {} {[@var{y0_new}, @var{yp0_new}] =} decic (@var{fcn}, @var{t0}, @var{y0}, @var{fixed_y0}, @var{yp0}, @var{fixed_yp0}, @var{options})
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
29 ## @deftypefnx {} {[@var{y0_new}, @var{yp0_new}, @var{resnorm}] =} decic (@dots{})
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
30 ##
25026
f886561f9696 doc: improve differential eqtn docs and mention ode15i/ode15s (bug #51965).
Colin Macdonald <cbm@m.fsf.org>
parents: 24534
diff changeset
31 ## Compute consistent implicit ODE initial conditions @var{y0_new} and
f886561f9696 doc: improve differential eqtn docs and mention ode15i/ode15s (bug #51965).
Colin Macdonald <cbm@m.fsf.org>
parents: 24534
diff changeset
32 ## @var{yp0_new} given initial guesses @var{y0} and @var{yp0}.
f886561f9696 doc: improve differential eqtn docs and mention ode15i/ode15s (bug #51965).
Colin Macdonald <cbm@m.fsf.org>
parents: 24534
diff changeset
33 ##
f886561f9696 doc: improve differential eqtn docs and mention ode15i/ode15s (bug #51965).
Colin Macdonald <cbm@m.fsf.org>
parents: 24534
diff changeset
34 ## A maximum of @code{length (@var{y0})} components between @var{fixed_y0} and
f886561f9696 doc: improve differential eqtn docs and mention ode15i/ode15s (bug #51965).
Colin Macdonald <cbm@m.fsf.org>
parents: 24534
diff changeset
35 ## @var{fixed_yp0} may be chosen as fixed values.
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
36 ##
30893
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files.
Rik <rik@octave.org>
parents: 30564
diff changeset
37 ## @var{fcn} is a function handle. The function must accept three inputs where
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
38 ## the first is time @var{t}, the second is a column vector of unknowns
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
39 ## @var{y}, and the third is a column vector of unknowns @var{yp}.
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
40 ##
25026
f886561f9696 doc: improve differential eqtn docs and mention ode15i/ode15s (bug #51965).
Colin Macdonald <cbm@m.fsf.org>
parents: 24534
diff changeset
41 ## @var{t0} is the initial time such that
30893
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files.
Rik <rik@octave.org>
parents: 30564
diff changeset
42 ## @code{@var{fcn}(@var{t0}, @var{y0_new}, @var{yp0_new}) = 0}, specified as a
25026
f886561f9696 doc: improve differential eqtn docs and mention ode15i/ode15s (bug #51965).
Colin Macdonald <cbm@m.fsf.org>
parents: 24534
diff changeset
43 ## scalar.
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
44 ##
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
45 ## @var{y0} is a vector used as the initial guess for @var{y}.
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
46 ##
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
47 ## @var{fixed_y0} is a vector which specifies the components of @var{y0} to
25026
f886561f9696 doc: improve differential eqtn docs and mention ode15i/ode15s (bug #51965).
Colin Macdonald <cbm@m.fsf.org>
parents: 24534
diff changeset
48 ## hold fixed. Choose a maximum of @code{length (@var{y0})} components between
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
49 ## @var{fixed_y0} and @var{fixed_yp0} as fixed values.
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
50 ## Set @var{fixed_y0}(i) component to 1 if you want to fix the value of
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
51 ## @var{y0}(i).
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
52 ## Set @var{fixed_y0}(i) component to 0 if you want to allow the value of
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
53 ## @var{y0}(i) to change.
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
54 ##
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
55 ## @var{yp0} is a vector used as the initial guess for @var{yp}.
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
56 ##
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
57 ## @var{fixed_yp0} is a vector which specifies the components of @var{yp0} to
25026
f886561f9696 doc: improve differential eqtn docs and mention ode15i/ode15s (bug #51965).
Colin Macdonald <cbm@m.fsf.org>
parents: 24534
diff changeset
58 ## hold fixed. Choose a maximum of @code{length (@var{yp0})} components
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
59 ## between @var{fixed_y0} and @var{fixed_yp0} as fixed values.
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
60 ## Set @var{fixed_yp0}(i) component to 1 if you want to fix the value of
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
61 ## @var{yp0}(i).
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
62 ## Set @var{fixed_yp0}(i) component to 0 if you want to allow the value of
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
63 ## @var{yp0}(i) to change.
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
64 ##
25026
f886561f9696 doc: improve differential eqtn docs and mention ode15i/ode15s (bug #51965).
Colin Macdonald <cbm@m.fsf.org>
parents: 24534
diff changeset
65 ## The optional seventh argument @var{options} is a structure array. Use
f886561f9696 doc: improve differential eqtn docs and mention ode15i/ode15s (bug #51965).
Colin Macdonald <cbm@m.fsf.org>
parents: 24534
diff changeset
66 ## @code{odeset} to generate this structure. The relevant options are
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
67 ## @code{RelTol} and @code{AbsTol} which specify the error thresholds used to
22910
23847979b91e maint: Strip trailing whitespace from source files.
John W. Eaton <jwe@octave.org>
parents: 22901
diff changeset
68 ## compute the initial conditions.
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
69 ##
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
70 ## The function typically returns two outputs. Variable @var{y0_new} is a
25579
07c2c42f457e doc: Miscellaneous documentation fixes all over the manual (bug #54288).
Rik <rik@octave.org>
parents: 25054
diff changeset
71 ## column vector and contains the consistent initial value of @var{y}. The
07c2c42f457e doc: Miscellaneous documentation fixes all over the manual (bug #54288).
Rik <rik@octave.org>
parents: 25054
diff changeset
72 ## output @var{yp0_new} is a column vector and contains the consistent initial
07c2c42f457e doc: Miscellaneous documentation fixes all over the manual (bug #54288).
Rik <rik@octave.org>
parents: 25054
diff changeset
73 ## value of @var{yp}.
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
74 ##
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
75 ## The optional third output @var{resnorm} is the norm of the vector of
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
76 ## residuals. If @var{resnorm} is small, @code{decic} has successfully
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
77 ## computed the initial conditions. If the value of @var{resnorm} is large,
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
78 ## use @code{RelTol} and @code{AbsTol} to adjust it.
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
79 ##
25026
f886561f9696 doc: improve differential eqtn docs and mention ode15i/ode15s (bug #51965).
Colin Macdonald <cbm@m.fsf.org>
parents: 24534
diff changeset
80 ## Example: Compute initial conditions for @nospell{Robertson's} equations:
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
81 ##
25026
f886561f9696 doc: improve differential eqtn docs and mention ode15i/ode15s (bug #51965).
Colin Macdonald <cbm@m.fsf.org>
parents: 24534
diff changeset
82 ## @smallexample
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
83 ## @group
25026
f886561f9696 doc: improve differential eqtn docs and mention ode15i/ode15s (bug #51965).
Colin Macdonald <cbm@m.fsf.org>
parents: 24534
diff changeset
84 ## function r = robertson_dae (@var{t}, @var{y}, @var{yp})
f886561f9696 doc: improve differential eqtn docs and mention ode15i/ode15s (bug #51965).
Colin Macdonald <cbm@m.fsf.org>
parents: 24534
diff changeset
85 ## r = [ -(@var{yp}(1) + 0.04*@var{y}(1) - 1e4*@var{y}(2)*@var{y}(3))
f886561f9696 doc: improve differential eqtn docs and mention ode15i/ode15s (bug #51965).
Colin Macdonald <cbm@m.fsf.org>
parents: 24534
diff changeset
86 ## -(@var{yp}(2) - 0.04*@var{y}(1) + 1e4*@var{y}(2)*@var{y}(3) + 3e7*@var{y}(2)^2)
f886561f9696 doc: improve differential eqtn docs and mention ode15i/ode15s (bug #51965).
Colin Macdonald <cbm@m.fsf.org>
parents: 24534
diff changeset
87 ## @var{y}(1) + @var{y}(2) + @var{y}(3) - 1 ];
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
88 ## endfunction
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
89 ## @end group
25026
f886561f9696 doc: improve differential eqtn docs and mention ode15i/ode15s (bug #51965).
Colin Macdonald <cbm@m.fsf.org>
parents: 24534
diff changeset
90 ## [@var{y0_new},@var{yp0_new}] = decic (@@robertson_dae, 0, [1; 0; 0], [1; 1; 0],
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
91 ## [-1e-4; 1; 0], [0; 0; 0]);
25026
f886561f9696 doc: improve differential eqtn docs and mention ode15i/ode15s (bug #51965).
Colin Macdonald <cbm@m.fsf.org>
parents: 24534
diff changeset
92 ## @end smallexample
22956
f9e6259c5fc1 decic.m: Fix bad @seealso macro from cset 80ac3e38b03d.
Rik <rik@octave.org>
parents: 22943
diff changeset
93 ## @seealso{ode15i, odeset}
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
94 ## @end deftypefn
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
95
30893
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files.
Rik <rik@octave.org>
parents: 30564
diff changeset
96 function [y0_new, yp0_new, resnorm] = decic (fcn, t0,
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
97 y0, fixed_y0, yp0, fixed_yp0,
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
98 options)
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
99
28789
28de41192f3c Eliminate unneeded verification of nargin, nargout in m-files.
Rik <rik@octave.org>
parents: 27923
diff changeset
100 if (nargin < 6)
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
101 print_usage ();
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
102 endif
22910
23847979b91e maint: Strip trailing whitespace from source files.
John W. Eaton <jwe@octave.org>
parents: 22901
diff changeset
103
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
104 ## Validate inputs
30893
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files.
Rik <rik@octave.org>
parents: 30564
diff changeset
105 if (! is_function_handle (fcn))
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
106 error ("Octave:invalid-input-arg",
30893
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files.
Rik <rik@octave.org>
parents: 30564
diff changeset
107 "decic: FCN must be a valid function handle");
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
108 endif
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
109
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
110 if (! (isnumeric (t0) && isscalar (t0)))
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
111 error ("Octave:invalid-input-arg",
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
112 "decic: T0 must be a numeric scalar");
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
113 endif
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
114
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
115 if ( ! (isnumeric (y0) && isvector (y0))
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
116 || ! (isnumeric (fixed_y0) && isvector (fixed_y0))
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
117 || ! (isnumeric (yp0) && isvector (yp0))
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
118 || ! (isnumeric (fixed_yp0) && isvector (fixed_yp0)))
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
119 error ("Octave:invalid-input-arg",
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
120 "decic: Y0, FIXED_Y0, YP0, and FIXED_YP0 must be numeric vectors");
22910
23847979b91e maint: Strip trailing whitespace from source files.
John W. Eaton <jwe@octave.org>
parents: 22901
diff changeset
121
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
122 elseif (! isequal (numel (y0), numel (fixed_y0), numel (yp0),
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
123 numel (fixed_yp0)))
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
124 error ("Octave:invalid-input-arg",
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
125 "decic: length of Y0, FIXED_Y0, YP0, and FIXED_YP0 must be equal");
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
126 endif
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
127
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
128 ## FIXME: This level of checking isn't necessary
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
129 for i = 1:numel (y0)
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
130 if (! (fixed_y0(i) == 0 || fixed_y0(i) == 1) || ! (fixed_yp0(i) == 0
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
131 || fixed_yp0(i) == 1))
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
132 error ("Octave:invalid-input-arg",
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
133 "decic: FIXED_Y0 and FIXED_YP0 must be boolean vectors");
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
134 endif
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
135 endfor
22910
23847979b91e maint: Strip trailing whitespace from source files.
John W. Eaton <jwe@octave.org>
parents: 22901
diff changeset
136
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
137 n = numel (y0);
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
138 nl = sum (! fixed_y0);
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
139 nu = sum (! fixed_yp0);
22910
23847979b91e maint: Strip trailing whitespace from source files.
John W. Eaton <jwe@octave.org>
parents: 22901
diff changeset
140
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
141 if (n - nl - nu > 0)
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
142 error ("Octave:invalid-input-arg",
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
143 "decic: cannot fix more than length (Y0) components");
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
144 endif
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
145
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
146 ## Set default values
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
147 TolFun = 0;
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
148 TolX = eps;
22910
23847979b91e maint: Strip trailing whitespace from source files.
John W. Eaton <jwe@octave.org>
parents: 22901
diff changeset
149
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
150 ## Check AbsTol and RelTol
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
151 if (nargin == 7)
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
152 if (! isempty (options.AbsTol))
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
153 if (! isscalar (options.AbsTol))
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
154 error ("Octave:invalid-input-arg",
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
155 "decic: AbsTol must be a scalar value");
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
156 else
22910
23847979b91e maint: Strip trailing whitespace from source files.
John W. Eaton <jwe@octave.org>
parents: 22901
diff changeset
157 TolFun = options.AbsTol;
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
158 endif
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
159 endif
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
160
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
161 if (! isempty (options.RelTol))
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
162 if (! isscalar (options.RelTol))
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
163 error ("Octave:invalid-input-arg",
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
164 "decic: RelTol must be a scalar value");
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
165 else
22910
23847979b91e maint: Strip trailing whitespace from source files.
John W. Eaton <jwe@octave.org>
parents: 22901
diff changeset
166 TolX = options.RelTol;
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
167 endif
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
168 endif
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
169 endif
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
170
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
171 x0 = [y0(! fixed_y0); yp0(! fixed_yp0)];
26153
7bd67e786e5f decic.m: Fix failing BIST test.
Rik <rik@octave.org>
parents: 25803
diff changeset
172 opt = optimset ("tolfun", TolFun, "tolx", TolX, "FinDiffType", "central");
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
173 x = ...
30893
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files.
Rik <rik@octave.org>
parents: 30564
diff changeset
174 fminunc (@(x) objective (x, t0, y0, fixed_y0, yp0, fixed_yp0, nl, nu, fcn),
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
175 x0, opt);
22910
23847979b91e maint: Strip trailing whitespace from source files.
John W. Eaton <jwe@octave.org>
parents: 22901
diff changeset
176
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
177 y0_new = y0;
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
178 yp0_new = yp0;
22910
23847979b91e maint: Strip trailing whitespace from source files.
John W. Eaton <jwe@octave.org>
parents: 22901
diff changeset
179
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
180 y0_new(! fixed_y0) = x(1:nl);
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
181 yp0_new(! fixed_yp0) = x(nl+1:nl+nu);
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
182 if (isargout (3))
30893
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files.
Rik <rik@octave.org>
parents: 30564
diff changeset
183 resnorm = fcn (t0, y0_new, yp0_new);
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
184 endif
22910
23847979b91e maint: Strip trailing whitespace from source files.
John W. Eaton <jwe@octave.org>
parents: 22901
diff changeset
185
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
186 endfunction
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
187
30893
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files.
Rik <rik@octave.org>
parents: 30564
diff changeset
188 function res = objective (x, t0, y0, fixed_y0, yp0, fixed_yp0, nl, nu, fcn)
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
189
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
190 y = y0;
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
191 y(! fixed_y0) = x(1:nl);
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
192 yp = yp0;
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
193 yp(! fixed_yp0) = x(nl+1:nl+nu);
30893
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files.
Rik <rik@octave.org>
parents: 30564
diff changeset
194 res = sqrt (sum (fcn (t0, y, yp) .^ 2));
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
195
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
196 endfunction
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
197
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
198
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
199 %!function res = rob (t, y, yp)
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
200 %! res =[-(yp(1) + 0.04*y(1) - 1e4*y(2)*y(3));
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
201 %! -(yp(2) - 0.04*y(1) + 1e4*y(2)*y(3) + 3e7*y(2)^2);
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
202 %! y(1) + y(2) + y(3) - 1];
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
203 %!endfunction
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
204
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
205 %!test # Without options
22910
23847979b91e maint: Strip trailing whitespace from source files.
John W. Eaton <jwe@octave.org>
parents: 22901
diff changeset
206 %! ref1 = [1;0;0];
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
207 %! ref2 = [-4e-2; 4e-2; 0];
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
208 %! [ynew, ypnew] = decic (@rob,0,[1;0;0],[1;1;0],[23;110;0],[0;0;1]);
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
209 %! assert ([ynew(1:end), ypnew(1:end)], [ref1(1:end), ref2(1:end)], 1e-10);
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
210 %!test # With options
22910
23847979b91e maint: Strip trailing whitespace from source files.
John W. Eaton <jwe@octave.org>
parents: 22901
diff changeset
211 %! ref1 = [1;0;0];
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
212 %! ref2 = [-4e-2; 4e-2; 0];
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
213 %! opt = odeset ("AbsTol", 1e-8, "RelTol", 1e-4);
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
214 %! [ynew, ypnew] = decic (@rob,0,[1;0;0],[1;1;0],[23;110;0],[0;0;1],opt);
22900
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
215 %! assert ([ynew(1:end), ypnew(1:end)], [ref1(1:end), ref2(1:end)], 1e-5);
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
216
ee1c77705fcd Refactor code for ode solvers and private functions
Francesco Faccio <francesco.faccio@mail.polimi.it>
parents:
diff changeset
217 ## Test input validation
28896
90fea9cc9caa test: Add expected error message <Invalid call> to BIST tests for nargin.
Rik <rik@octave.org>
parents: 28789
diff changeset
218 %!error <Invalid call> decic ()
90fea9cc9caa test: Add expected error message <Invalid call> to BIST tests for nargin.
Rik <rik@octave.org>
parents: 28789
diff changeset
219 %!error <Invalid call> decic (1)
90fea9cc9caa test: Add expected error message <Invalid call> to BIST tests for nargin.
Rik <rik@octave.org>
parents: 28789
diff changeset
220 %!error <Invalid call> decic (1,2)
90fea9cc9caa test: Add expected error message <Invalid call> to BIST tests for nargin.
Rik <rik@octave.org>
parents: 28789
diff changeset
221 %!error <Invalid call> decic (1,2,3)
90fea9cc9caa test: Add expected error message <Invalid call> to BIST tests for nargin.
Rik <rik@octave.org>
parents: 28789
diff changeset
222 %!error <Invalid call> decic (1,2,3,4)
90fea9cc9caa test: Add expected error message <Invalid call> to BIST tests for nargin.
Rik <rik@octave.org>
parents: 28789
diff changeset
223 %!error <Invalid call> decic (1,2,3,4,5)
30893
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files.
Rik <rik@octave.org>
parents: 30564
diff changeset
224 %!error <FCN must be a valid function handle>
22943
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
225 %! decic (1, 0, [1; 0; 0], [1; 1; 0], [-1e-4; 1; 0], [0; 0; 0]);
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
226 %!error <T0 must be a numeric scalar>
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
227 %! decic (@rob, [1, 1], [1; 0; 0], [1; 1; 0], [-1e-4; 1; 0], [0; 0; 0]);
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
228 %!error <length of Y0, FIXED_Y0, YP0, and FIXED_YP0 must be equal>
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
229 %! decic (@rob, 0, [0; 0], [1; 1; 0], [-1e-4; 1; 0], [0; 0; 0]);
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
230 %!error <Y0, FIXED_Y0, YP0, and FIXED_YP0 must be numeric vectors>
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
231 %! decic (@rob, 0, [1; 0; 0], [1; 0],"", [0; 0; 0]);
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
232 %!error <Y0, FIXED_Y0, YP0, and FIXED_YP0 must be numeric vectors>
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
233 %! decic (@rob, 0, [1; 0; 0], [1; 1; 0], [-1e-4; 1; 0], [0; 0; "1"]);
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
234 %!error <FIXED_Y0 and FIXED_YP0 must be boolean vectors>
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
235 %! decic (@rob, 0, [1; 0; 0], [5; 5; 0], [-1e-4; 1; 0], [0; 0; 0]);
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
236 %!error <FIXED_Y0 and FIXED_YP0 must be boolean vectors>
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
237 %! decic (@rob, 0, [1; 0; 0], [1; 1; 0], [-1e-4; 1; 0], [0; 4; 0]);
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
238 %!error <cannot fix more than length \(Y0\) components>
80ac3e38b03d decic.m: Use Octave coding conventions.
Rik <rik@octave.org>
parents: 22910
diff changeset
239 %! decic (@rob, 0, [1; 0; 0], [1; 1; 0], [-1e-4; 1; 0], [0; 1; 1]);