Mercurial > octave

annotate scripts/ode/private/runge_kutta_interpolate.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) 2015-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 ##
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
6 ##
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
7 ## This file is part of Octave.
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
8 ##
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23407
diff changeset
9 ## Octave is free software: you can redistribute it and/or modify it
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
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: 23407
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.
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
13 ##
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
14 ## Octave is distributed in the hope that it will be useful, but
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
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.
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
18 ##
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
19 ## You should have received a copy of the GNU General Public License
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
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: 23407
diff changeset
21 ## <https://www.gnu.org/licenses/>.
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
22
20847
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
23 function u_interp = runge_kutta_interpolate (order, z, u, t, k_vals, dt, func, args)
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
24
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
25 switch (order)
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
26
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
27 case 1
22655
6b134d294d61 ode solvers: use ordinary transpose instead of Hermitian conjugate (bug #49410).
Carlo de Falco <carlo.defalco@polimi.it>
parents: 22626
diff changeset
28 u_interp = interp1 (z, u.', t, "linear");
20847
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
29
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
30 case 2
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
31 if (! isempty (k_vals))
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
32 der = k_vals(:,1);
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
33 else
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
34 der = feval (func, z(1) , u(:,1), args);
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
35 endif
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
36 u_interp = quadratic_interpolation (z, u, der, t);
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
37
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
38 case 3
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
39 u_interp = hermite_cubic_interpolation (z, u, k_vals, t);
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
40
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
41 case 5
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
42 ## ode45 with Dormand-Prince scheme:
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
43 u_interp = hermite_quartic_interpolation (z, u, k_vals, t);
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
44
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
45 otherwise
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
46 warning (["High order interpolation not yet implemented: ", ...
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
47 "using cubic interpolation instead"]);
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
48 der(:,1) = feval (func, z(1), u(:,1), args);
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
49 der(:,2) = feval (func, z(2), u(:,2), args);
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
50 u_interp = hermite_cubic_interpolation (z, u, der, t);
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
51
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
52 endswitch
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
53
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
54 endfunction
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
55
20847
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
56 ## The function below can be used in an ODE solver to interpolate the solution
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
57 ## at the time t_out using 2nd order Hermite interpolation.
20847
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
58 function x_out = quadratic_interpolation (t, x, der, t_out)
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
59
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
60 # coefficients of the parabola
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
61 a = -(x(:,1) - x(:,2) - der(:).*(t(1)-t(2))) / (t(1) - t(2))^2;
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
62 b = der(:) - 2*t(1).*a;
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
63 c = x(:,1) - a*t(1)^2 - b*t(1);
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
64
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
65 # evaluate in t_out
20847
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
66 x_out = a*t_out.^2 + b*t_out + c;
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
67
ddc18b909ec7 codesprint: ode: first and second order interpolator for dense_output
jcorno <jacopo.corno@gmail.com>
parents: 20634
diff changeset
68 endfunction
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
69
20903
3d3da166dac5 2015 Code Sprint: finish import of ode23 into core
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20847
diff changeset
70 ## The function below can be used in an ODE solver to interpolate the
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
71 ## solution at the time t_out using 3rd order Hermite interpolation.
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
72 function x_out = hermite_cubic_interpolation (t, x, der, t_out)
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
73
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
74 dt = (t(2) - t(1));
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
75 s = (t_out - t(1)) / dt;
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
76 x_out = ((1 + 2*s) .* (1-s).^2) .* x(:,1) + ...
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
77 (s .* (1-s).^2 * dt ) .* der(:,1) + ...
23407
e265ae9e7a6c Fix hermite cubic interpolation in ode23
Andreas Stahel <Andreas.Stahel@bfh.ch>
parents: 23220
diff changeset
78 ((3-2*s) .* s.^2 ) .* x(:,end) + ...
e265ae9e7a6c Fix hermite cubic interpolation in ode23
Andreas Stahel <Andreas.Stahel@bfh.ch>
parents: 23220
diff changeset
79 ((s-1) .* s.^2 * dt ) .* der(:,end);
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
80
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
81 endfunction
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
82
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
83 ## The function below can be used in an ODE solver to interpolate the
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
84 ## solution at the time t_out using 4th order Hermite interpolation.
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
85 function x_out = hermite_quartic_interpolation (t, x, der, t_out)
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
86
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
87 persistent coefs_u_half = ...
22626
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
88 [6025192743/30085553152; 0; 51252292925/65400821598;
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
89 -2691868925/45128329728; 187940372067/1594534317056;
869c02fde46c Further clean-up of ode functions.
Rik <rik@octave.org>
parents: 22323
diff changeset
90 -1776094331/19743644256; 11237099/235043384];
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
91
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
92 ## 4th order approximation of y in t+dt/2 as proposed by Shampine in
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
93 ## Lawrence, Shampine, "Some Practical Runge-Kutta Formulas", 1986.
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
94 dt = t(2) - t(1);
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
95 u_half = x(:,1) + (1/2) * dt * (der(:,1:7) * coefs_u_half);
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20621
diff changeset
96
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
97 ## Rescale time on [0,1]
20567
ea6a1c00763a fix interpolation bug introduced with 87b557ee8e5d
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20565
diff changeset
98 s = (t_out - t(1)) / dt;
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
99
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
100 ## Hermite basis functions
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
101 ## H0 = 1 - 11*s.^2 + 18*s.^3 - 8*s.^4;
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
102 ## H1 = s - 4*s.^2 + 5*s.^3 - 2*s.^4;
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
103 ## H2 = 16*s.^2 - 32*s.^3 + 16*s.^4;
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
104 ## H3 = - 5*s.^2 + 14*s.^3 - 8*s.^4;
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
105 ## H4 = s.^2 - 3*s.^3 + 2*s.^4;
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
106
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
107 x_out = (1 - 11*s.^2 + 18*s.^3 - 8*s.^4) .* x(:,1) + ...
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
108 ( s - 4*s.^2 + 5*s.^3 - 2*s.^4) .* (dt * der(:,1)) + ...
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
109 ( 16*s.^2 - 32*s.^3 + 16*s.^4) .* u_half + ...
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
110 ( - 5*s.^2 + 14*s.^3 - 8*s.^4) .* x(:,2) + ...
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20567
diff changeset
111 ( s.^2 - 3*s.^3 + 2*s.^4) .* (dt * der(:,end));
20565
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
112
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
113 endfunction