# HG changeset patch # User Carlo de Falco # Date 1444152539 -7200 # Node ID 87b557ee8e5d2956db0417bbeb4e4fc2ce32c3e8 # Parent c1a6c31ac29a30f2693ca73b964819fc0944187d 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. diff -r c1a6c31ac29a -r 87b557ee8e5d scripts/ode/module.mk --- 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 diff -r c1a6c31ac29a -r 87b557ee8e5d scripts/ode/private/hermite_quartic_interpolation.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 -## -## 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 -## . - -## -*- 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 - diff -r c1a6c31ac29a -r 87b557ee8e5d scripts/ode/private/integrate_adaptive.m --- 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' ## ## 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)]; diff -r c1a6c31ac29a -r 87b557ee8e5d scripts/ode/private/ode_rk_interpolate.m --- /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 +## +## 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 +## . + +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