Mercurial > octave

annotate scripts/ode/private/runge_kutta_45_dorpri.m @ 27919:1891570abac8

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 2020.
author John W. Eaton <jwe@octave.org>
date Mon, 06 Jan 2020 22:29:51 -0500
parents b442ec6dda5c
children bd51beb6205e
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
27919
1891570abac8 update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents: 27918
diff changeset
1 ## Copyright (C) 2013-2020 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
2 ##
b442ec6dda5c use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents: 26376
diff changeset
3 ## See the file COPYRIGHT.md in the top-level directory of this distribution
b442ec6dda5c use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents: 26376
diff changeset
4 ## or <https://octave.org/COPYRIGHT.html/>.
b442ec6dda5c use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents: 26376
diff changeset
5 ##
20536
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
6 ##
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
7 ## 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
8 ##
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23565
diff changeset
9 ## 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
10 ## 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: 23565
diff changeset
11 ## the Free Software Foundation, either version 3 of the License, or
22755
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22661
diff changeset
12 ## (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
13 ##
20536
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
14 ## 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
15 ## WITHOUT ANY WARRANTY; without even the implied warranty of
22755
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22661
diff changeset
16 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22661
diff changeset
17 ## 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
18 ##
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
19 ## 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
20 ## 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: 23565
diff changeset
21 ## <https://www.gnu.org/licenses/>.
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
22
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
23 ## -*- texinfo -*-
22768
a75a0900e216 maint: Match documentation names and variable names in more m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
24 ## @deftypefn {} {[@var{t_next}, @var{x_next}] =} runge_kutta_45_dorpri (@var{@@fun}, @var{t}, @var{x}, @var{dt})
23565
3a730821e4a2 doc: Peridoc grammarcheck of documentation.
Rik <rik@octave.org>
parents: 23220
diff changeset
25 ## @deftypefnx {} {[@var{t_next}, @var{x_next}] =} runge_kutta_45_dorpri (@var{@@fun}, @var{t}, @var{x}, @var{dt}, @var{options})
3a730821e4a2 doc: Peridoc grammarcheck of documentation.
Rik <rik@octave.org>
parents: 23220
diff changeset
26 ## @deftypefnx {} {[@var{t_next}, @var{x_next}] =} runge_kutta_45_dorpri (@var{@@fun}, @var{t}, @var{x}, @var{dt}, @var{options}, @var{k_vals})
3a730821e4a2 doc: Peridoc grammarcheck of documentation.
Rik <rik@octave.org>
parents: 23220
diff changeset
27 ## @deftypefnx {} {[@var{t_next}, @var{x_next}] =} runge_kutta_45_dorpri (@var{@@fun}, @var{t}, @var{x}, @var{dt}, @var{options}, @var{k_vals}, @var{t_next})
22768
a75a0900e216 maint: Match documentation names and variable names in more m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
28 ## @deftypefnx {} {[@var{t_next}, @var{x_next}, @var{x_est}] =} runge_kutta_45_dorpri (@dots{})
a75a0900e216 maint: Match documentation names and variable names in more m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
29 ## @deftypefnx {} {[@var{t_next}, @var{x_next}, @var{x_est}, @var{k_vals_out}] =} runge_kutta_45_dorpri (@dots{})
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
30 ##
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
31 ## This function can be used to integrate a system of ODEs with a given initial
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20600
diff changeset
32 ## condition @var{x} from @var{t} to @var{t+dt} with the Dormand-Prince method.
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
33 ## For the definition of this method see
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
34 ## @url{http://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method}.
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
35 ##
20548
25623ef2ff4f doc: Rewrite docstrings for ode* family of functions.
Rik <rik@octave.org>
parents: 20543
diff changeset
36 ## First input argument is the function describing the system of ODEs to be
25623ef2ff4f doc: Rewrite docstrings for ode* family of functions.
Rik <rik@octave.org>
parents: 20543
diff changeset
37 ## integrated.
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
38 ##
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
39 ## Second input parameter is the first extreme of integration interval.
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
40 ##
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
41 ## Third input argument is the initial condition of the system.
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 ## Fourth input argument is the timestep, that is the length of the
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
44 ## integration interval.
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
45 ##
20548
25623ef2ff4f doc: Rewrite docstrings for ode* family of functions.
Rik <rik@octave.org>
parents: 20543
diff changeset
46 ## Fifth input parameter is optional and describes a set of options useful to
25623ef2ff4f doc: Rewrite docstrings for ode* family of functions.
Rik <rik@octave.org>
parents: 20543
diff changeset
47 ## adapt the computation to what is needed.
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 ## Sixth input parameter is optional and describes the Runge-Kutta evaluations
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20600
diff changeset
50 ## of the previous step to use in an FSAL scheme.
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
51 ##
22768
a75a0900e216 maint: Match documentation names and variable names in more m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
52 ## Seventh input parameter is optional and is the time (@var{t_next}) to
a75a0900e216 maint: Match documentation names and variable names in more m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
53 ## integrate to. The default is @code{@var{t} + @var{dt}}.
a75a0900e216 maint: Match documentation names and variable names in more m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
54 ##
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
55 ## First output argument is the final integration time value.
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
56 ##
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
57 ## Second output parameter is the higher order computed solution at time
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
58 ## @var{t_next} (local extrapolation).
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
59 ##
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
60 ## Third output parameter is a lower order solution for the estimation of the
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
61 ## error.
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
62 ##
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
63 ## Fourth output parameter is matrix containing the Runge-Kutta evaluations
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
64 ## to use in an FSAL scheme or for dense output.
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
65 ## @end deftypefn
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
66
22768
a75a0900e216 maint: Match documentation names and variable names in more m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
67 function [t_next, x_next, x_est, k] = runge_kutta_45_dorpri (fun, t, x, dt,
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20600
diff changeset
68 options = [],
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20600
diff changeset
69 k_vals = [],
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20600
diff changeset
70 t_next = t + dt)
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
71
22661
c00578ac8dda Resolve FIXME notes in ode code base.
Rik <rik@octave.org>
parents: 22626
diff changeset
72 ## Reference: Hairer, Ernst; Nørsett, Syvert Paul; Wanner, Gerhard (2008),
c00578ac8dda Resolve FIXME notes in ode code base.
Rik <rik@octave.org>
parents: 22626
diff changeset
73 ## Solving ordinary differential equations I: Nonstiff problems,
c00578ac8dda Resolve FIXME notes in ode code base.
Rik <rik@octave.org>
parents: 22626
diff changeset
74 ## Berlin, New York: Springer-Verlag, ISBN 978-3-540-56670-0
20540
72cd24aa5f7a Clean up and vetorize Dormant&Prince RK timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20536
diff changeset
75 persistent a = [0 0 0 0 0 0;
72cd24aa5f7a Clean up and vetorize Dormant&Prince RK timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20536
diff changeset
76 1/5 0 0 0 0 0;
72cd24aa5f7a Clean up and vetorize Dormant&Prince RK timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20536
diff changeset
77 3/40 9/40 0 0 0 0;
72cd24aa5f7a Clean up and vetorize Dormant&Prince RK timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20536
diff changeset
78 44/45 -56/15 32/9 0 0 0;
72cd24aa5f7a Clean up and vetorize Dormant&Prince RK timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20536
diff changeset
79 19372/6561 -25360/2187 64448/6561 -212/729 0 0;
20903
3d3da166dac5 2015 Code Sprint: finish import of ode23 into core
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20852
diff changeset
80 9017/3168 -355/33 46732/5247 49/176 -5103/18656 0];
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
81 persistent b = [0, 1/5, 3/10, 4/5, 8/9, 1, 1];
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
82 persistent c = [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84];
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
83 persistent c_prime = [5179/57600, 0, 7571/16695, 393/640, ...
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
84 -92097/339200, 187/2100, 1/40];
20540
72cd24aa5f7a Clean up and vetorize Dormant&Prince RK timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20536
diff changeset
85
72cd24aa5f7a Clean up and vetorize Dormant&Prince RK timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20536
diff changeset
86 s = t + dt * b;
72cd24aa5f7a Clean up and vetorize Dormant&Prince RK timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20536
diff changeset
87 cc = dt * c;
72cd24aa5f7a Clean up and vetorize Dormant&Prince RK timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20536
diff changeset
88 aa = dt * a;
72cd24aa5f7a Clean up and vetorize Dormant&Prince RK timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20536
diff changeset
89 k = zeros (rows (x), 7);
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
90
22661
c00578ac8dda Resolve FIXME notes in ode code base.
Rik <rik@octave.org>
parents: 22626
diff changeset
91 if (! isempty (options)) # extra arguments for function evaluator
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
92 args = options.funarguments;
20543
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20540
diff changeset
93 else
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20540
diff changeset
94 args = {};
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
95 endif
20543
3339c9bdfe6a Activate FSAL property in dorpri timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20540
diff changeset
96
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20600
diff changeset
97 if (! isempty (k_vals)) # k values from previous step are passed
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20554
diff changeset
98 k(:,1) = k_vals(:,end); # FSAL property
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
99 else
22768
a75a0900e216 maint: Match documentation names and variable names in more m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
100 k(:,1) = feval (fun, t, x, args{:});
20600
a22d8a2eb0e5 fix adaptive strategy in ode solvers.
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20554
diff changeset
101 endif
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
102
22768
a75a0900e216 maint: Match documentation names and variable names in more m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
103 k(:,2) = feval (fun, s(2), x + k(:,1) * aa(2, 1).' , args{:});
a75a0900e216 maint: Match documentation names and variable names in more m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
104 k(:,3) = feval (fun, s(3), x + k(:,1:2) * aa(3, 1:2).', args{:});
a75a0900e216 maint: Match documentation names and variable names in more m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
105 k(:,4) = feval (fun, s(4), x + k(:,1:3) * aa(4, 1:3).', args{:});
a75a0900e216 maint: Match documentation names and variable names in more m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
106 k(:,5) = feval (fun, s(5), x + k(:,1:4) * aa(5, 1:4).', args{:});
a75a0900e216 maint: Match documentation names and variable names in more m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
107 k(:,6) = feval (fun, s(6), x + k(:,1:5) * aa(6, 1:5).', args{:});
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
108
20552
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20548
diff changeset
109 ## compute new time and new values for the unknowns
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20600
diff changeset
110 ## t_next = t + dt;
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20600
diff changeset
111 x_next = x + k(:,1:6) * cc(:); # 5th order approximation
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
112
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
113 ## if the estimation of the error is required
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
114 if (nargout >= 3)
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
115 ## new solution to be compared with the previous one
22768
a75a0900e216 maint: Match documentation names and variable names in more m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
116 k(:,7) = feval (fun, t_next, x_next, args{:});
20540
72cd24aa5f7a Clean up and vetorize Dormant&Prince RK timestepper
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20536
diff changeset
117 cc_prime = dt * c_prime;
20552
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20548
diff changeset
118 x_est = x + k * cc_prime(:);
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
119 endif
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
120
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
121 endfunction