diff scripts/ode/private/runge_kutta_45_dorpri.m @ 20543:3339c9bdfe6a

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.
author Carlo de Falco <carlo.defalco@polimi.it>
date Sat, 03 Oct 2015 07:32:50 +0200
parents 72cd24aa5f7a
children 25623ef2ff4f
line wrap: on
line diff
--- a/scripts/ode/private/runge_kutta_45_dorpri.m	Fri Oct 02 15:07:37 2015 -0400
+++ b/scripts/ode/private/runge_kutta_45_dorpri.m	Sat Oct 03 07:32:50 2015 +0200
@@ -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