Mercurial > hg > octave-nkf

annotate scripts/ode/private/ode_rk_interpolate.m @ 20796:e5f36a7854a5

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Remove fuzzy matching from odeset/odeget. * levenshtein.cc: Deleted file. * libinterp/corefcn/module.mk: Remove levenshtein.cc from build system. * fuzzy_compare.m: Deleted file. * scripts/ode/module.mk: Remove fuzzy_compare.m from build system * odeget.m: Reword docstring. Use a persistent cellstr variable to keep track of all options. Replace fuzzy_compare() calls with combination of strcmpi and strncmpi. Report errors relative to function odeget rather than OdePkg. Rewrite and extend BIST tests. Add input validation BIST tests. * odeset.m: Reword docstring. Use a persistent cellstr variable to keep track of all options. Replace fuzzy_compare() calls with combination of strcmpi and strncmpi. Report errors relative to function odeset rather than OdePkg. Use more meaningful variables names and create intermediate variables with logical names to help make code readable. Remove interactive input when multiple property names match and just issue an error. Rewrite BIST tests. * ode_struct_value_check.m: Remove input checking for private function which must always be invoked correctly by caller. Use intermediate variables opt and val to make the code more understandable. Consolidate checks on values into single if statements. Use 'val == fix (val)' to check for integer. * __unimplemented__.m: Removed odeset, odeget, ode45 from list.
author Rik <rik@octave.org>
date Fri, 09 Oct 2015 12:03:23 -0700
parents ea6a1c00763a
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
20771
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
1 ## Copyright (C) 2015 Carlo de Falco
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
2 ## Copyright (C) 2015 Jacopo Corno <jacopo.corno@gmail.com>
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
3 ##
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
4 ## 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
5 ##
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
6 ## Octave is free software; you can redistribute it and/or modify it
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
7 ## under the terms of the GNU General Public License as published by
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
8 ## the Free Software Foundation; either version 3 of the License, or (at
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
9 ## your option) any later version.
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
10 ##
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
11 ## 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
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
14 ## General Public License for more details.
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
15 ##
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
16 ## 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
17 ## along with Octave; see the file COPYING. If not, see
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
18 ## <http://www.gnu.org/licenses/>.
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
19
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
20 function u_interp = ode_rk_interpolate (order, z, u, t, k_vals, dt, args)
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
21
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
22 switch order
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
23
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
24 #{
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
25 case 1
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
26 u_interp = linear_interpolation (z, u, t);
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
27 case 2
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
28 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
29 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
30 else
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
31 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
32 endif
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
33 u_interp = quadratic_interpolation (z, u, der, t);
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
34 case 3
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
35 u_interp = ...
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
36 hermite_cubic_interpolation (z, u, k_vals, t);
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
37 case 4
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
38 ## if ode45 is used without local extrapolation this function
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
39 ## doesn't require a new function evaluation.
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
40 u_interp = dorpri_interpolation ([z(i-1) z(i)],
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
41 [u(:,i-1) u(:,i)],
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
42 k_vals, tspan(counter));
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
43
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
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
46 case 5
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
47 ## ode45 with Dormand-Prince scheme:
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
48 u_interp = ...
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
49 hermite_quartic_interpolation (z, u, k_vals, t);
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
50
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
51 ## it is also possible to do a new function evaluation and use
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
52 ## the quintic hermite interpolator
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
53 ## f_half = feval (func, t+1/2*dt, u_half,
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
54 ## options.vfunarguments{:});
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
55 ## u_interp =
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
56 ## hermite_quintic_interpolation ([z(i-1) z(i)],
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
57 ## [u(:,i-1) u_half u(:,i)],
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
58 ## [k_vals(:,1) f_half ...
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
59 ## k_vals(:,end)],
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
60 ## tspan(counter));
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
61 otherwise
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
62 warning ("High order interpolation not yet implemented: ",
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
63 "using cubic interpolation instead");
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
64 der(:,1) = 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
65 der(:,2) = feval (func, z(2) , u(:,2), args);
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
66 u_interp = hermite_cubic_interpolation (z, u, der, t);
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
67 endswitch
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
68
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
69 endfunction
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
70
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
71
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
72
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
73 ## The function below can be used in an ODE solver to
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
74 ## interpolate the solution at the time t_out using 4th order
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
75 ## hermite interpolation.
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
76 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
77
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
78 persistent coefs_u_half = ...
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
79 [(6025192743/30085553152), 0, (51252292925/65400821598), ...
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
80 (-2691868925/45128329728), (187940372067/1594534317056), ...
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
81 (-1776094331/19743644256), (11237099/235043384)].';
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
82
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
83 ## 4th order approximation of y in t+dt/2 as proposed by
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
84 ## Shampine in Lawrence, Shampine, "Some Practical
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
85 ## Runge-Kutta Formulas", 1986.
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
86 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
87 u_half = x(:,1) + (1/2) * dt * (der(:,1:7) * coefs_u_half);
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
88
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
89 ## Rescale time on [0,1]
20773
ea6a1c00763a fix interpolation bug introduced with 87b557ee8e5d
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20771
diff changeset
90 s = (t_out - t(1)) / dt;
20771
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
91
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
92 ## Hermite basis functions
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
93 ## 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
94 ## 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
95 ## 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
96 ## 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
97 ## 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
98
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
99 x_out = zeros (rows (x), length (t_out));
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
100 x_out = (1 - 11*s.^2 + 18*s.^3 - 8*s.^4) .* x(:,1) + ...
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
101 ( s - 4*s.^2 + 5*s.^3 - 2*s.^4) .* (dt * der(:,1)) + ...
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
102 ( 16*s.^2 - 32*s.^3 + 16*s.^4) .* u_half + ...
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
103 ( - 5*s.^2 + 14*s.^3 - 8*s.^4) .* x(:,2) + ...
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
104 ( s.^2 - 3*s.^3 + 2*s.^4) .* (dt * der(:,end));
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
105
87b557ee8e5d clean up and vectorize code for dense output in ode45
Carlo de Falco <carlo.defalco@polimi.it>
parents:
diff changeset
106 endfunction