changeset 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 c1a6c31ac29a
children bc6daa38ff50
files scripts/ode/module.mk scripts/ode/private/hermite_quartic_interpolation.m scripts/ode/private/integrate_adaptive.m scripts/ode/private/ode_rk_interpolate.m
diffstat 4 files changed, 114 insertions(+), 134 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/ode/module.mk	Tue Oct 06 00:20:02 2015 -0400
+++ b/scripts/ode/module.mk	Tue Oct 06 19:28:59 2015 +0200
@@ -5,13 +5,13 @@
 scripts_ode_PRIVATE_FCN_FILES = \
   scripts/ode/private/AbsRel_Norm.m \
   scripts/ode/private/fuzzy_compare.m \
-  scripts/ode/private/hermite_quartic_interpolation.m \
   scripts/ode/private/integrate_adaptive.m \
   scripts/ode/private/integrate_const.m \
   scripts/ode/private/integrate_n_steps.m \
   scripts/ode/private/kahan.m \
   scripts/ode/private/odepkg_event_handle.m \
   scripts/ode/private/odepkg_structure_check.m \
+  scripts/ode/private/ode_rk_interpolate.m \
   scripts/ode/private/ode_struct_value_check.m \
   scripts/ode/private/runge_kutta_45_dorpri.m \
   scripts/ode/private/starting_stepsize.m
--- a/scripts/ode/private/hermite_quartic_interpolation.m	Tue Oct 06 00:20:02 2015 -0400
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,65 +0,0 @@
-## Copyright (C) 2015 Jacopo Corno <jacopo.corno@gmail.com>
-##
-## This file is part of Octave.
-##
-## Octave is free software; you can redistribute it and/or modify it
-## under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 3 of the License, or (at
-## your option) any later version.
-##
-## Octave is distributed in the hope that it will be useful, but
-## WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-## General Public License for more details.
-##
-## You should have received a copy of the GNU General Public License
-## along with Octave; see the file COPYING.  If not, see
-## <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn {Function File} {@var{x_out} =} hermite_quartic_interpolation (@var{t}, @var{x}, @var{der}, @var{t_out})
-##
-## This function file can be called by an ODE solver function in order to
-## interpolate the solution at the time @var{t_out} using a 4th order
-## Hermite interpolation.
-##
-## The first input @var{t} is a vector with two given times.
-##
-## The second input argument is the vector with the values of the function
-## to interpolate at the times specified in @var{t} and at the middle point.
-##
-## The third input argument is the value of the derivatives of the function
-## evaluated at the two extreme points.
-##
-## The output @var{x_out} is the evaluation of the Hermite interpolant at
-## @var{t_out}.
-##
-## @end deftypefn
-##
-## @seealso{linear_interpolation, quadratic_interpolation,
-## hermite_cubic_interpolation, hermite_quintic_interpolation,
-## dorpri_interpolation}
-
-function x_out = hermite_quartic_interpolation (t, x, der, t_out)
-
-  ## Rescale time on [0,1]
-  s = (t_out - t(1)) / (t(2) - t(1));
-
-  ## Hermite basis functions
-  ## H0 = 1   - 11*s.^2 + 18*s.^3 -  8*s.^4;
-  ## H1 =   s -  4*s.^2 +  5*s.^3 -  2*s.^4;
-  ## H2 =       16*s.^2 - 32*s.^3 + 16*s.^4;
-  ## H3 =     -  5*s.^2 + 14*s.^3 -  8*s.^4;
-  ## H4 =          s.^2 -  3*s.^3 +  2*s.^4;
-
-  x_out = zeros (rows (x), length (t_out));
-  for ii = 1:rows (x)
-    x_out(ii,:) = (1   - 11*s.^2 + 18*s.^3 -  8*s.^4)*x(ii,1) ...
-                + (  s -  4*s.^2 +  5*s.^3 -  2*s.^4)*(t(2)-t(1))*der(ii,1) ...
-                + (      16*s.^2 - 32*s.^3 + 16*s.^4)*x(ii,2) ...
-                + (    -  5*s.^2 + 14*s.^3 -  8*s.^4)*x(ii,3) ...
-                + (         s.^2 -  3*s.^3 +  2*s.^4)*(t(2)-t(1))*der(ii,2);
-  endfor
-
-endfunction
-
--- a/scripts/ode/private/integrate_adaptive.m	Tue Oct 06 00:20:02 2015 -0400
+++ b/scripts/ode/private/integrate_adaptive.m	Tue Oct 06 19:28:59 2015 +0200
@@ -1,3 +1,4 @@
+## Copyright (C) 2015, Carlo de Falco
 ## Copyright (C) 2013, Roberto Porcu' <roberto.porcu@polimi.it>
 ##
 ## This file is part of Octave.
@@ -147,8 +148,7 @@
     ## Solution accepted only if the error is less or equal to 1.0
     if (err <= 1)
       
-      [tk, comp] = kahan (tk, comp, dt);
-      options.comp = comp;
+      [tk, options.comp] = kahan (tk, options.comp, dt);
       s(end) = tk;
 
       ## values on this interval for time and solution
@@ -182,73 +182,12 @@
           while ((counter <= k)
                  && (vdirection * z(i) > vdirection * tspan(counter)))
             ## choose interpolation scheme according to order of the solver
-            switch (order)
-              case 1
-                u_interp = linear_interpolation ([z(i-1) z(i)],
-                                                 [u(:,i-1) u(:,i)],
-                                                 tspan(counter));
-              case 2
-                if (! isempty (k_vals))
-                  der = k_vals(:,1);
-                else
-                  der = feval (func, z(i-1) , u(:,i-1),
-                               options.vfunarguments{:});
-                endif
-                u_interp = quadratic_interpolation ([z(i-1) z(i)],
-                                                    [u(:,i-1) u(:,i)],
-                                                    der, tspan(counter));
-              case 3
-                u_interp = ...
-                hermite_cubic_interpolation ([z(i-1) z(i)],
-                                             [u(:,i-1) u(:,i)],
-                                             [k_vals(:,1) k_vals(:,end)],
-                                             tspan(counter));
-              case 4
-                ## if ode45 is used without local extrapolation this function
-                ## doesn't require a new function evaluation.
-                u_interp = dorpri_interpolation ([z(i-1) z(i)],
-                                                 [u(:,i-1) u(:,i)],
-                                                 k_vals, tspan(counter));
-              case 5
-                ## ode45 with Dormand-Prince scheme:
-                ## 4th order approximation of y in t+dt/2 as proposed by
-                ## Shampine in Lawrence, Shampine, "Some Practical
-                ## Runge-Kutta Formulas", 1986.
-                u_half = u(:,i-1) ...
-                         + 1/2*dt*((6025192743/30085553152) * k_vals(:,1)
-                                   + (51252292925/65400821598) * k_vals(:,3)
-                                   - (2691868925/45128329728) * k_vals(:,4)
-                                   + (187940372067/1594534317056) * k_vals(:,5)
-                                   - (1776094331/19743644256) * k_vals(:,6)
-                                   + (11237099/235043384) * k_vals(:,7));
-                u_interp = ...
-                hermite_quartic_interpolation ([z(i-1) z(i)],
-                                               [u(:,i-1) u_half u(:,i)],
-                                               [k_vals(:,1) k_vals(:,end)],
-                                               tspan(counter));
 
-                ## it is also possible to do a new function evaluation and use
-                ## the quintic hermite interpolator
-                ## f_half = feval (func, t+1/2*dt, u_half,
-                ##                 options.vfunarguments{:});
-                ## u_interp =
-                ##   hermite_quintic_interpolation ([z(i-1) z(i)],
-                ##                                  [u(:,i-1) u_half u(:,i)],
-                ##                                  [k_vals(:,1) f_half ...
-                ##                                   k_vals(:,end)],
-                ##                                  tspan(counter));
-              otherwise
-                warning ("High order interpolation not yet implemented: ",
-                         "using cubic interpolation instead");
-                der(:,1) = feval (func, z(i-1) , u(:,i-1),
-                                  options.vfunarguments{:});
-                der(:,2) = feval (func, z(i) , u(:,i),
-                                  options.vfunarguments{:});
-                u_interp = ...
-                hermite_cubic_interpolation ([z(i-1) z(i)],
-                                             [u(:,i-1) u(:,i)],
-                                             der, tspan(counter));
-            endswitch
+            u_interp = ...
+            ode_rk_interpolate (order, [z(i-1) z(i)], [u(:,i-1) u(:,i)],
+                                tspan(counter), k_vals, dt,
+                                options.vfunarguments{:});
+            
 
             ## add the interpolated value of the solution
             u = [u(:,1:i-1), u_interp, u(:,i:end)];
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/ode/private/ode_rk_interpolate.m	Tue Oct 06 19:28:59 2015 +0200
@@ -0,0 +1,106 @@
+## Copyright (C) 2015 Carlo de Falco
+## Copyright (C) 2015 Jacopo Corno <jacopo.corno@gmail.com>
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+
+function u_interp = ode_rk_interpolate (order, z, u, t, k_vals, dt, args)
+
+  switch order
+
+    #{
+    case 1
+      u_interp = linear_interpolation (z, u, t);
+    case 2
+      if (! isempty (k_vals))
+        der = k_vals(:,1);
+      else
+        der = feval (func, z(1) , u(:,1), args);
+      endif
+      u_interp = quadratic_interpolation (z, u, der, t);
+    case 3
+      u_interp = ...
+      hermite_cubic_interpolation (z, u, k_vals, t);
+    case 4
+      ## if ode45 is used without local extrapolation this function
+      ## doesn't require a new function evaluation.
+      u_interp = dorpri_interpolation ([z(i-1) z(i)],
+                                       [u(:,i-1) u(:,i)],
+                                       k_vals, tspan(counter));
+
+    #}
+         
+    case 5
+      ## ode45 with Dormand-Prince scheme:     
+      u_interp = ...
+      hermite_quartic_interpolation (z, u, k_vals, t);
+
+      ## it is also possible to do a new function evaluation and use
+      ## the quintic hermite interpolator
+      ## f_half = feval (func, t+1/2*dt, u_half,
+      ##                 options.vfunarguments{:});
+      ## u_interp =
+      ##   hermite_quintic_interpolation ([z(i-1) z(i)],
+      ##                                  [u(:,i-1) u_half u(:,i)],
+      ##                                  [k_vals(:,1) f_half ...
+      ##                                   k_vals(:,end)],
+      ##                                  tspan(counter));
+    otherwise
+      warning ("High order interpolation not yet implemented: ",
+               "using cubic interpolation instead");
+      der(:,1) = feval (func, z(1) , u(:,1), args);
+      der(:,2) = feval (func, z(2) , u(:,2), args);
+      u_interp = hermite_cubic_interpolation (z, u, der, t);
+  endswitch
+
+endfunction
+
+
+
+## The function below can be used in an ODE solver to
+## interpolate the solution at the time t_out using 4th order
+## hermite interpolation.
+function x_out = hermite_quartic_interpolation (t, x, der, t_out)
+
+  persistent coefs_u_half = ...
+  [(6025192743/30085553152), 0, (51252292925/65400821598), ...
+   (-2691868925/45128329728), (187940372067/1594534317056), ...
+   (-1776094331/19743644256), (11237099/235043384)].';
+
+  ## 4th order approximation of y in t+dt/2 as proposed by
+  ## Shampine in Lawrence, Shampine, "Some Practical
+  ## Runge-Kutta Formulas", 1986.
+  dt = t(2) - t(1);
+  u_half = x(:,1) + (1/2) * dt * (der(:,1:7) * coefs_u_half);
+  
+  ## Rescale time on [0,1]
+  s = (t_out - t) / dt;
+
+  ## Hermite basis functions
+  ## H0 = 1   - 11*s.^2 + 18*s.^3 -  8*s.^4;
+  ## H1 =   s -  4*s.^2 +  5*s.^3 -  2*s.^4;
+  ## H2 =       16*s.^2 - 32*s.^3 + 16*s.^4;
+  ## H3 =     -  5*s.^2 + 14*s.^3 -  8*s.^4;
+  ## H4 =          s.^2 -  3*s.^3 +  2*s.^4;
+
+  x_out = zeros (rows (x), length (t_out));
+  x_out = (1   - 11*s.^2 + 18*s.^3 -  8*s.^4)   .* x(:,1) + ...
+          (  s -  4*s.^2 +  5*s.^3 -  2*s.^4)   .* (dt * der(:,1)) + ...
+          (      16*s.^2 - 32*s.^3 + 16*s.^4)   .* u_half + ...
+          (    -  5*s.^2 + 14*s.^3 -  8*s.^4)   .* x(:,2) + ...
+          (         s.^2 -  3*s.^3 +  2*s.^4)   .* (dt * der(:,end));
+
+endfunction