Mercurial > octave-nkf
annotate scripts/ode/private/ode_rk_interpolate.m @ 20598:ea6a1c00763a
fix interpolation bug introduced with 87b557ee8e5d
* ode_rk_interpolate.m(hermite_quartic_interpolation):
fix typo that lead to incorrect size for interpolated values.
author | Carlo de Falco <carlo.defalco@polimi.it> |
---|---|
date | Tue, 06 Oct 2015 22:14:41 +0200 |
parents | 87b557ee8e5d |
children |
rev | line source |
---|---|
20596
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] |
20598
ea6a1c00763a
fix interpolation bug introduced with 87b557ee8e5d
Carlo de Falco <carlo.defalco@polimi.it>
parents:
20596
diff
changeset
|
90 s = (t_out - t(1)) / dt; |
20596
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 |