# HG changeset patch # User Carlo de Falco # Date 1443849665 -7200 # Sat Oct 03 07:21:05 2015 +0200 # Node ID 2fad9d8cae73801fa3be1f7daf76d2c518c6f137 # Parent dd6345fd8a976c2487ee1fbc62c63353d45f1132 Activate FSAL property in dorpri timestepper * scripts/ode/private/runge_kutta_45_dorpri.m: don't compute first stage if values from previous iteration are passed. * scripts/ode/private/integrate_adaptive.m: do not update cmputed stages if timestep is rejected. diff --git a/scripts/ode/ode45.m b/scripts/ode/ode45.m --- a/scripts/ode/ode45.m +++ b/scripts/ode/ode45.m @@ -106,8 +106,8 @@ endif if (length (vslot) < 2 && ... - (isempty (vodeoptions.TimeStepSize) ... - || isempty (vodeoptions.TimeStepNumber))) + (isempty (vodeoptions.TimeStepSize) ... + || isempty (vodeoptions.TimeStepNumber))) error ("OdePkg:InvalidArgument", ... "second input argument must be a valid vector"); elseif (vslot(2) == vslot(1)) @@ -251,8 +251,8 @@ ## option can be set by the user to another value than default value. if (isempty (vodeoptions.InitialStep) && strcmp (integrate_func, "adaptive")) vodeoptions.InitialStep = vodeoptions.vdirection* ... - starting_stepsize (vorder, vfun, vslot(1), vinit, vodeoptions.AbsTol, ... - vodeoptions.RelTol, vodeoptions.vnormcontrol); + starting_stepsize (vorder, vfun, vslot(1), vinit, vodeoptions.AbsTol, ... + vodeoptions.RelTol, vodeoptions.vnormcontrol); warning ("OdePkg:InvalidArgument", ... "option ''InitialStep'' not set, estimated value %f is used", ... vodeoptions.InitialStep); @@ -359,11 +359,11 @@ if (vmassdependence) ## constant mass matrices have already vmass = @(t,x) vodeoptions.Mass (t, x, vodeoptions.vfunarguments{:}); vfun = @(t,x) vmass (t, x, vodeoptions.vfunarguments{:}) ... - \ vfun (t, x, vodeoptions.vfunarguments{:}); + \ vfun (t, x, vodeoptions.vfunarguments{:}); else ## if (vmassdependence == false) vmass = @(t) vodeoptions.Mass (t, vodeoptions.vfunarguments{:}); vfun = @(t,x) vmass (t, vodeoptions.vfunarguments{:}) ... - \ vfun (t, x, vodeoptions.vfunarguments{:}); + \ vfun (t, x, vodeoptions.vfunarguments{:}); endif endif @@ -385,11 +385,11 @@ ## Postprocessing, do whatever when terminating integration algorithm if (vodeoptions.vhaveoutputfunction) ## Cleanup plotter feval (vodeoptions.OutputFcn, solution.t(end), ... - solution.x(end,:)', "done", vodeoptions.vfunarguments{:}); + solution.x(end,:)', "done", vodeoptions.vfunarguments{:}); endif if (vodeoptions.vhaveeventfunction) ## Cleanup event function handling odepkg_event_handle (vodeoptions.Events, solution.t(end), ... - solution.x(end,:)', "done", vodeoptions.vfunarguments{:}); + solution.x(end,:)', "done", vodeoptions.vfunarguments{:}); endif ## Print additional information if option Stats is set diff --git a/scripts/ode/private/integrate_adaptive.m b/scripts/ode/private/integrate_adaptive.m --- a/scripts/ode/private/integrate_adaptive.m +++ b/scripts/ode/private/integrate_adaptive.m @@ -69,7 +69,7 @@ ## first values for time and solution t = tspan(1); x = x0(:); - + ## get first initial timestep dt = odeget (options, "InitialStep", starting_stepsize (order, func, t, x, options.AbsTol, @@ -120,15 +120,20 @@ z = t; u = x; - k_vals = feval (func, t , x, options.vfunarguments{:}); - + k_vals = []; + while (counter <= k) facmax = 1.5; ## compute integration step from t to t+dt - [s, y, y_est, k_vals] = stepper (func, z(end), u(:,end), - dt, options, k_vals); - + if (isempty (k_vals)) + [s, y, y_est, new_k_vals] = stepper (func, z(end), u(:,end), + dt, options); + else + [s, y, y_est, new_k_vals] = stepper (func, z(end), u(:,end), + dt, options, k_vals); + endif + if (options.vhavenonnegative) x(options.NonNegative,end) = abs (x(options.NonNegative,end)); y(options.NonNegative,end) = abs (y(options.NonNegative,end)); @@ -141,27 +146,29 @@ err = AbsRel_Norm (y(:,end), u(:,end), options.AbsTol, options.RelTol, options.vnormcontrol, y_est(:,end)); - + ## 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; s(end) = tk; ## values on this interval for time and solution - z = [z(end);s]; - u = [u(:,end),y]; - + z = [z(end); s]; + u = [u(:,end), y]; + k_vals = new_k_vals; + ## if next tspan value is caught, update counter if ((z(end) == tspan(counter)) || (abs (z(end) - tspan(counter)) / (max (abs (z(end)), abs (tspan(counter)))) < 8*eps) ) counter++; - - ## if there is an element in time vector at which the solution is required - ## the program must compute this solution before going on with next steps + + ## if there is an element in time vector at which the solution is required + ## the program must compute this solution before going on with next steps elseif (vdirection * z(end) > vdirection * tspan(counter)) + ## initialize counter for the following cycle i = 2; while (i <= length (z)) @@ -179,9 +186,9 @@ ## 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)); + 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); @@ -194,10 +201,10 @@ 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)); + 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. @@ -217,10 +224,10 @@ - (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)); + 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 @@ -239,16 +246,16 @@ 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)); + hermite_cubic_interpolation ([z(i-1) z(i)], + [u(:,i-1) u(:,i)], + der, tspan(counter)); endswitch ## add the interpolated value of the solution u = [u(:,1:i-1), u_interp, u(:,i:end)]; ## add the time requested - z = [z(1:i-1);tspan(counter);z(i:end)]; + z = [z(1:i-1); tspan(counter); z(i:end)]; ## update counters counter++; @@ -310,7 +317,7 @@ ## true if (options.vhaveeventfunction) solution.vevent = odepkg_event_handle (options.Events, t(end), - x(:,end), [], options.vfunarguments{:}); + x(:,end), [], options.vfunarguments{:}); if (! isempty (solution.vevent{1}) && solution.vevent{1} == 1) t(solution.vcntloop-1,:) = solution.vevent{3}(end,:); @@ -321,13 +328,15 @@ endif else + facmax = 1.0; + endif ## Compute next timestep, formula taken from Hairer err += eps; # adding an eps to avoid divisions by zero - dt = dt * min (facmax, max (facmin, - fac * (1 / err)^(1 / (order + 1)))); + dt = dt * min (facmax, + max (facmin, fac * (1 / err)^(1 / (order + 1)))); dt = vdirection * min (abs (dt), options.MaxStep); ## Update counters that count the number of iteration cycles @@ -383,7 +392,3 @@ solution.x = x(:,1:end-f)'; endfunction - -## Local Variables: *** -## mode: octave *** -## End: *** diff --git a/scripts/ode/private/runge_kutta_45_dorpri.m b/scripts/ode/private/runge_kutta_45_dorpri.m --- a/scripts/ode/private/runge_kutta_45_dorpri.m +++ b/scripts/ode/private/runge_kutta_45_dorpri.m @@ -63,7 +63,8 @@ ## ## @seealso{odepkg} -function varargout = runge_kutta_45_dorpri (f, t, x, dt, varargin) +function [t_out, x_out, x_est, k] = ... + runge_kutta_45_dorpri (f, t, x, dt, varargin) persistent a = [0 0 0 0 0 0; 1/5 0 0 0 0 0; @@ -83,17 +84,19 @@ s = t + dt * b; cc = dt * c; aa = dt * a; - - args = varargin{1}.vfunarguments; k = zeros (rows (x), 7); - if (nargin == 5) # only the options are passed - k(:,1) = feval (f, t, x, args{:}); - elseif (nargin == 6) # both the options and the k values are passed - k(:,1) = varargin{2}(:,end); # FSAL property + if (nargin >= 5) # options are passed + args = varargin{1}.vfunarguments; + if (nargin >= 6) # both the options and the k values are passed + k(:,1) = varargin{2}(:,end); # FSAL property + else + k(:,1) = feval (f, t, x, args{:}); + endif + else + args = {}; endif - - k(:,1) = feval (f, s(1), x, args{:}); + k(:,2) = feval (f, s(2), x + k(:,1) * aa(2, 1).', args{:}); k(:,3) = feval (f, s(3), x + k(:,1:2) * aa(3, 1:2).', args{:}); k(:,4) = feval (f, s(4), x + k(:,1:3) * aa(4, 1:3).', args{:}); @@ -101,16 +104,15 @@ k(:,6) = feval (f, s(6), x + k(:,1:5) * aa(6, 1:5).', args{:}); ## compute new time and new values for the unkwnowns - varargout{1} = t + dt; - varargout{2} = x + k(:,1:6) * cc(:); # 5th order approximation + t_out = t + dt; + x_out = x + k(:,1:6) * cc(:); # 5th order approximation ## if the estimation of the error is required if (nargout >= 3) ## new solution to be compared with the previous one - k(:,7) = feval (f, t + dt, varargout{2}, args{:}); + k(:,7) = feval (f, t + dt, x_out, args{:}); cc_prime = dt * c_prime; - varargout{3} = x + k * cc_prime(:); # x_est - varargout{4} = k; + x_est = x + k * cc_prime(:); # x_est endif endfunction