diff scripts/ode/private/runge_kutta_45_dorpri.m @ 20635:a22d8a2eb0e5

fix adaptive strategy in ode solvers. * script/ode/ode45.m: remove unused option OutputSave * script/ode/private/integrate_adaptive.m: rewrite algorithm to be more compatible. * script/ode/private/runge_kutta_45_dorpri.m: use kahan summation for time increment.
author Carlo de Falco <carlo.defalco@polimi.it>
date Sun, 11 Oct 2015 18:44:58 +0200
parents b7ac1e94266e
children
line wrap: on
line diff
--- a/scripts/ode/private/runge_kutta_45_dorpri.m	Sat Oct 10 16:52:59 2015 -0700
+++ b/scripts/ode/private/runge_kutta_45_dorpri.m	Sun Oct 11 18:44:58 2015 +0200
@@ -58,7 +58,8 @@
 ## @seealso{odepkg}
 
 function [t_out, x_out, x_est, k] = ...
-         runge_kutta_45_dorpri (f, t, x, dt, varargin)
+         runge_kutta_45_dorpri (f, t, x, dt, opts = [], k_vals = [],
+                                t_out = t + dt)
 
   persistent a = [0           0          0           0        0          0;
                   1/5         0          0           0        0          0;
@@ -80,17 +81,18 @@
   aa = dt * a;
   k = zeros (rows (x), 7);
 
-  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
+  if (! isempty (opts)) # extra arguments for function evaluator
+    args = opts.vfunarguments;
   else
     args = {};
   endif
 
+  if (! isempty (k_vals))  # k values from previous step are passed
+    k(:,1) = k_vals(:,end);  # FSAL property
+  else      
+    k(:,1) = feval (f, t, x, args{:});
+  endif
+    
   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{:});
@@ -98,13 +100,13 @@
   k(:,6) = feval (f, s(6), x + k(:,1:5) * aa(6, 1:5).', args{:});
 
   ## compute new time and new values for the unknowns
-  t_out = t + dt;
+  ## 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, x_out, args{:});
+    k(:,7) = feval (f, t_out, x_out, args{:});
     cc_prime = dt * c_prime;
     x_est = x + k * cc_prime(:);
   endif