Mercurial > octave

annotate scripts/ode/private/runge_kutta_23s.m @ 29359:7854d5752dd2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: merge stable to default.
author John W. Eaton <jwe@octave.org>
date Wed, 10 Feb 2021 10:10:40 -0500
parents 967cfcde2e35 0a5b15007766
children 796f54d4ddbf
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
28308
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
1 ########################################################################
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
2 ##
29358
0a5b15007766 update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents: 28308
diff changeset
3 ## Copyright (C) 2013-2021 The Octave Project Developers
28308
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
4 ##
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
5 ## See the file COPYRIGHT.md in the top-level directory of this
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
6 ## distribution or <https://octave.org/copyright/>.
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
7 ##
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
8 ## This file is part of Octave.
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
9 ##
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
10 ## Octave is free software: you can redistribute it and/or modify it
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
11 ## under the terms of the GNU General Public License as published by
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
12 ## the Free Software Foundation, either version 3 of the License, or
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
13 ## (at your option) any later version.
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
14 ##
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
15 ## Octave is distributed in the hope that it will be useful, but
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
16 ## WITHOUT ANY WARRANTY; without even the implied warranty of
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
17 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
18 ## GNU General Public License for more details.
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
19 ##
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
20 ## You should have received a copy of the GNU General Public License
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
21 ## along with Octave; see the file COPYING. If not, see
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
22 ## <https://www.gnu.org/licenses/>.
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
23 ##
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
24 ########################################################################
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
25
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
26 ## -*- texinfo -*-
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
27 ## @deftypefn {} {[@var{t_next}, @var{x_next}] =} runge_kutta_23s (@var{fun}, @var{t}, @var{x}, @var{dt})
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
28 ## @deftypefnx {} {[@var{t_next}, @var{x_next}] =} runge_kutta_23s (@var{fun}, @var{t}, @var{x}, @var{dt}, @var{options})
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
29 ## @deftypefnx {} {[@var{t_next}, @var{x_next}] =} runge_kutta_23s (@var{fun}, @var{t}, @var{x}, @var{dt}, @var{options}, @var{k_vals})
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
30 ## @deftypefnx {} {[@var{t_next}, @var{x_next}] =} runge_kutta_23s (@var{fun}, @var{t}, @var{x}, @var{dt}, @var{options}, @var{k_vals}, @var{t_next})
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
31 ## @deftypefnx {} {[@var{t_next}, @var{x_next}, @var{x_est}] =} runge_kutta_23s (@dots{})
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
32 ## @deftypefnx {} {[@var{t_next}, @var{x_next}, @var{x_est}, @var{k_vals_out}] =} runge_kutta_23s (@dots{})
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
33 ##
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
34 ## This function can be used to integrate a system of ODEs with a given initial
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
35 ## condition @var{x} from @var{t} to @var{t+dt}, with a Rosenbrock method of
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
36 ## order (2,3). All the mathematical formulas are from Shampine, L. F. and
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
37 ## M. W. Reichelt, "The MATLAB ODE Suite", SIAM Journal on Scientific
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
38 ## Computing, Vol. 18, 1997, pp. 1–22.
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
39 ##
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
40 ## @var{f} is a function handle that defines the ODE: @code{y' = f(tau,y)}.
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
41 ## The function must accept two inputs where the first is time @var{tau} and
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
42 ## the second is a column vector of unknowns @var{y}.
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
43 ##
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
44 ## @var{t} is the first extreme of integration interval.
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
45 ##
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
46 ## @var{x} is the initial condition of the system..
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
47 ##
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
48 ## @var{dt} is the timestep, that is the length of the integration interval.
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
49 ##
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
50 ## The optional fourth argument @var{options} specifies options for the ODE
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
51 ## solver. It is a structure generated by @code{odeset}. In particular it
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
52 ## contains the field @var{funarguments} with the optional arguments to be used
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
53 ## in the evaluation of @var{fun}.
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
54 ##
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
55 ## The optional fifth argument @var{k_vals_in} contains the Runge-Kutta
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
56 ## evaluations of the previous step to use in a FSAL scheme.
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
57 ##
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
58 ## The optional sixth argument @var{t_next} (@code{t_next = t + dt}) specifies
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
59 ## the end of the integration interval. The output @var{x_next} s the higher
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
60 ## order computed solution at time @var{t_next} (local extrapolation is
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
61 ## performed).
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
62 ##
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
63 ## Optionally the functions can also return @var{x_est}, a lower order solution
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
64 ## for the estimation of the error, and @var{k_vals_out}, a matrix containing
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
65 ## the Runge-Kutta evaluations to use in a FSAL scheme or for dense output.
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
66 ##
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
67 ## @seealso{runge_kutta_23}
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
68 ## @end deftypefn
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
69
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
70 function [t_next, x_next, x_est, k] = runge_kutta_23s (fun, t, x, dt,
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
71 options = [],
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
72 k_vals = [],
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
73 t_next = t + dt)
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
74
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
75 persistent d = 1 / (2 + sqrt (2));
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
76 persistent a = 1 / 2;
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
77 persistent e32 = 6 + sqrt (2);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
78
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
79 ## extra arguments for function evaluator
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
80 if (! isempty (options))
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
81 args = options.funarguments;
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
82 else
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
83 args = {};
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
84 endif
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
85
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
86 jacfun = false;
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
87 jacmat = false;
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
88 if (! isempty (options.Jacobian))
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
89 if (ischar (options.Jacobian))
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
90 jacfun = true;
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
91 jac = str2fun (options.Jacobian);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
92 elseif (is_function_handle (options.Jacobian))
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
93 jacfun = true;
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
94 jac = options.Jacobian;
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
95 elseif (ismatrix (options.Jacobian))
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
96 jacmat = true;
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
97 jac = options.Jacobian;
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
98 else
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
99 error (["ode23s: the jacobian should be passed as a matrix, ", ...
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
100 "a string or a function handle"])
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
101 endif
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
102 endif
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
103
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
104 jacpat = false;
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
105 if (! isempty (options.JPattern))
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
106 jacpat = true;
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
107 pattern = logical (options.JPattern);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
108 endif
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
109
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
110 ## Jacobian matrix, dfxpdp
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
111 if (jacmat)
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
112 J = jac;
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
113 elseif (jacfun)
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
114 J = jac (t, x);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
115 elseif (! jacpat)
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
116 J = __dfxpdp__ (x, @(a) feval (fun, t, a, args{:}), options.RelTol);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
117 elseif (jacpat)
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
118 J = __dfxpdp__ (x, @(a) feval (fun, t, a, args{:}), options.RelTol, pattern);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
119 endif
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
120
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
121 T = (feval (fun, t + .1 * dt, x) - feval (fun, t, x, args{:})) / (.1 * dt);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
122
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
123 ## Wolfbrandt coefficient
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
124 if (isempty (options.Mass))
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
125 M = speye (length (x));
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
126 else
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
127 M = options.Mass;
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
128 endif
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
129 W = M - dt*d*J;
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
130
28921
967cfcde2e35 maint: Use parentheses around conditional expressions.
Rik <rik@octave.org>
parents: 28308
diff changeset
131 if (issparse (W))
28308
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
132 [Lw, Uw, Pw, Qw, Rw] = lu (W);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
133 else
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
134 [Lw, Uw, Pw] = lu (W);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
135 endif
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
136
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
137 ## compute the slopes
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
138 F(:,1) = feval (fun, t, x, args{:});
28921
967cfcde2e35 maint: Use parentheses around conditional expressions.
Rik <rik@octave.org>
parents: 28308
diff changeset
139 if (issparse (W))
28308
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
140 k(:,1) = Qw * (Uw \ (Lw \ (Pw * (Rw \ (F(:,1) + dt*d*T)))));
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
141 else
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
142 k(:,1) = Uw \ (Lw \ (Pw * (F(:,1) + dt*d*T)));
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
143 endif
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
144 F(:,2) = feval (fun, t+a*dt, x+a*dt*k(:,1), args{:});
28921
967cfcde2e35 maint: Use parentheses around conditional expressions.
Rik <rik@octave.org>
parents: 28308
diff changeset
145 if (issparse (W))
28308
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
146 k(:,2) = Uw * (Uw \ (Lw \ (Pw * (Rw \ (F(:,2) - M*k(:,1)))))) + k(:,1);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
147 else
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
148 k(:,2) = Uw \ (Lw \ (Pw * (F(:,2) - M*k(:,1)))) + k(:,1);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
149 endif
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
150
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
151 ## compute the 2nd order estimate
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
152 x_next = x + dt*k(:,2);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
153
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
154 if (nargout >= 3)
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
155 ## 3rd order, needed in error formula
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
156 F(:,3) = feval (fun, t+dt, x_next, args{:});
28921
967cfcde2e35 maint: Use parentheses around conditional expressions.
Rik <rik@octave.org>
parents: 28308
diff changeset
157 if (issparse (W))
28308
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
158 k(:,3) = Qw * (Uw \ (Lw \ (Pw * (Rw \ (F(:,3) - e32 * (M*k(:,2) - F(:,2)) - 2 * (M*k(:,1) - F(:,1)) + dt*d*T)))));
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
159 else
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
160 k(:,3) = Uw \ (Lw \ (Pw * (F(:,3) - e32 * (M*k(:,2) - F(:,2)) - 2 * (M*k(:,1) - F(:,1)) + dt*d*T)));
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
161 endif
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
162
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
163 ## estimate the error
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
164 err_est = (dt/6) * (k(:,1) - 2*k(:,2) + k(:,3));
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
165
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
166 ## FIXME: to use in AbsRel_Norm function I need x_est and not err directly
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
167 x_est = x_next + err_est;
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
168 endif
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
169
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
170 endfunction
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
171
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
172
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
173 function prt = __dfxpdp__ (p, func, rtol, pattern)
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
174
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
175 ## This subfunction was copied 2011 from the OF "optim" package
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
176 ## "inst/private/__dfdp__.m".
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
177
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
178 f = func (p)(:);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
179 m = numel (f);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
180 n = numel (p);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
181
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
182 diffp = rtol .* ones (n, 1);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
183
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
184 del = ifelse (p == 0, diffp, diffp .* p);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
185 absdel = abs (del);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
186
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
187 ## double sided interval
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
188 p1 = p + absdel / 2;
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
189 p2 = p - absdel / 2;
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
190
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
191 ps = p;
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
192 if (nargin > 3 && issparse (pattern))
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
193 prt = pattern; # initialize Jacobian
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
194 for j = find (any (pattern, 1))
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
195 ps(j) = p1(j);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
196 tp1 = func (ps);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
197 ps(j) = p2(j);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
198 tp2 = func (ps);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
199 pattern_nnz = find (pattern(:, j));
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
200 prt(pattern_nnz, j) = (tp1(pattern_nnz) - tp2(pattern_nnz)) / absdel(j);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
201 ps(j) = p(j);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
202 endfor
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
203 else
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
204 prt = zeros (m, n); # initialize Jacobian
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
205 for j = 1:n
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
206 ps(j) = p1(j);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
207 tp1 = func (ps);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
208 ps(j) = p2(j);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
209 tp2 = func (ps);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
210 prt(:, j) = (tp1(:) - tp2(:)) / absdel(j);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
211 ps(j) = p(j);
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
212 endfor
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
213 endif
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
214
8085ae13cc4a ode23s.m: new function from former odepkg (bug #57309)
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
diff changeset
215 endfunction