Mercurial > octave-nkf
comparison scripts/ode/private/ode_rk_interpolate.m @ 20596:87b557ee8e5d
clean up and vectorize code for dense output in ode45
* scripts/ode/private/ode_rk_interpolate.m: new file
* scripts/ode/private/ode_rk_interpolate.m(hermite_quartic_interpolation):
move to internal function, use vectorization and broadcasting.
* scripts/ode/private/hermite_quartic_interpolation.m: remove file
* scripts/ode/module.mk: list added and removed files
* scripts/ode/private/integrate_adaptive.m: use new interpolation code.
author | Carlo de Falco <carlo.defalco@polimi.it> |
---|---|
date | Tue, 06 Oct 2015 19:28:59 +0200 |
parents | |
children | ea6a1c00763a |
comparison
equal
deleted
inserted
replaced
20595:c1a6c31ac29a | 20596:87b557ee8e5d |
---|---|
1 ## Copyright (C) 2015 Carlo de Falco | |
2 ## Copyright (C) 2015 Jacopo Corno <jacopo.corno@gmail.com> | |
3 ## | |
4 ## This file is part of Octave. | |
5 ## | |
6 ## Octave is free software; you can redistribute it and/or modify it | |
7 ## under the terms of the GNU General Public License as published by | |
8 ## the Free Software Foundation; either version 3 of the License, or (at | |
9 ## your option) any later version. | |
10 ## | |
11 ## Octave is distributed in the hope that it will be useful, but | |
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
14 ## General Public License for more details. | |
15 ## | |
16 ## You should have received a copy of the GNU General Public License | |
17 ## along with Octave; see the file COPYING. If not, see | |
18 ## <http://www.gnu.org/licenses/>. | |
19 | |
20 function u_interp = ode_rk_interpolate (order, z, u, t, k_vals, dt, args) | |
21 | |
22 switch order | |
23 | |
24 #{ | |
25 case 1 | |
26 u_interp = linear_interpolation (z, u, t); | |
27 case 2 | |
28 if (! isempty (k_vals)) | |
29 der = k_vals(:,1); | |
30 else | |
31 der = feval (func, z(1) , u(:,1), args); | |
32 endif | |
33 u_interp = quadratic_interpolation (z, u, der, t); | |
34 case 3 | |
35 u_interp = ... | |
36 hermite_cubic_interpolation (z, u, k_vals, t); | |
37 case 4 | |
38 ## if ode45 is used without local extrapolation this function | |
39 ## doesn't require a new function evaluation. | |
40 u_interp = dorpri_interpolation ([z(i-1) z(i)], | |
41 [u(:,i-1) u(:,i)], | |
42 k_vals, tspan(counter)); | |
43 | |
44 #} | |
45 | |
46 case 5 | |
47 ## ode45 with Dormand-Prince scheme: | |
48 u_interp = ... | |
49 hermite_quartic_interpolation (z, u, k_vals, t); | |
50 | |
51 ## it is also possible to do a new function evaluation and use | |
52 ## the quintic hermite interpolator | |
53 ## f_half = feval (func, t+1/2*dt, u_half, | |
54 ## options.vfunarguments{:}); | |
55 ## u_interp = | |
56 ## hermite_quintic_interpolation ([z(i-1) z(i)], | |
57 ## [u(:,i-1) u_half u(:,i)], | |
58 ## [k_vals(:,1) f_half ... | |
59 ## k_vals(:,end)], | |
60 ## tspan(counter)); | |
61 otherwise | |
62 warning ("High order interpolation not yet implemented: ", | |
63 "using cubic interpolation instead"); | |
64 der(:,1) = feval (func, z(1) , u(:,1), args); | |
65 der(:,2) = feval (func, z(2) , u(:,2), args); | |
66 u_interp = hermite_cubic_interpolation (z, u, der, t); | |
67 endswitch | |
68 | |
69 endfunction | |
70 | |
71 | |
72 | |
73 ## The function below can be used in an ODE solver to | |
74 ## interpolate the solution at the time t_out using 4th order | |
75 ## hermite interpolation. | |
76 function x_out = hermite_quartic_interpolation (t, x, der, t_out) | |
77 | |
78 persistent coefs_u_half = ... | |
79 [(6025192743/30085553152), 0, (51252292925/65400821598), ... | |
80 (-2691868925/45128329728), (187940372067/1594534317056), ... | |
81 (-1776094331/19743644256), (11237099/235043384)].'; | |
82 | |
83 ## 4th order approximation of y in t+dt/2 as proposed by | |
84 ## Shampine in Lawrence, Shampine, "Some Practical | |
85 ## Runge-Kutta Formulas", 1986. | |
86 dt = t(2) - t(1); | |
87 u_half = x(:,1) + (1/2) * dt * (der(:,1:7) * coefs_u_half); | |
88 | |
89 ## Rescale time on [0,1] | |
90 s = (t_out - t) / dt; | |
91 | |
92 ## Hermite basis functions | |
93 ## H0 = 1 - 11*s.^2 + 18*s.^3 - 8*s.^4; | |
94 ## H1 = s - 4*s.^2 + 5*s.^3 - 2*s.^4; | |
95 ## H2 = 16*s.^2 - 32*s.^3 + 16*s.^4; | |
96 ## H3 = - 5*s.^2 + 14*s.^3 - 8*s.^4; | |
97 ## H4 = s.^2 - 3*s.^3 + 2*s.^4; | |
98 | |
99 x_out = zeros (rows (x), length (t_out)); | |
100 x_out = (1 - 11*s.^2 + 18*s.^3 - 8*s.^4) .* x(:,1) + ... | |
101 ( s - 4*s.^2 + 5*s.^3 - 2*s.^4) .* (dt * der(:,1)) + ... | |
102 ( 16*s.^2 - 32*s.^3 + 16*s.^4) .* u_half + ... | |
103 ( - 5*s.^2 + 14*s.^3 - 8*s.^4) .* x(:,2) + ... | |
104 ( s.^2 - 3*s.^3 + 2*s.^4) .* (dt * der(:,end)); | |
105 | |
106 endfunction |