Mercurial > octave

diff scripts/ode/private/runge_kutta_45_dorpri.m @ 20621:b92f8e148936

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Continued clean-up of functions in ode/private dir. * AbsRel_Norm.m: Use retval as return variable. Use sumsq() rather than explicit squaring of vector and sum(). Combine multiple lines where possible. * integrate_adaptive.m: Rewrite docstring. Only call starting_stepsize() if InitialStep option is empty. * integrate_const.m: Rewrite docstring. Remove useless commented out code. Combine multiple lines where possible. * integrate_n_steps.m: Rewrite docstring. Remove useless commented out code. Combine multiple lines where possible. * kahan.m: Remove excessive 4-space indentation, use 2-space indentation. * ode_rk_interpolate.m: Use parentheses around condition for switch stmt. Combine multiple lines where possible. * ode_struct_value_check.m: Remove comma from Copyright statement that Octave doesn't use. * odepkg_event_handle.m: Remove comma from Copyright statement that Octave doesn't use. * odepkg_structure_check.m: Remove comma from Copyright statement that Octave doesn't use. * runge_kutta_45_dorpri.m: Remove comma from Copyright statement that Octave doesn't use. Improve docstring. Match variable names in documentation to those in code. * starting_stepsize.m: Rewrite docstring. Use spaces between function name and opening parenthesis.
author Rik <rik@octave.org>
date Wed, 14 Oct 2015 10:35:53 -0700
parents a22d8a2eb0e5
children 80e630b37ba1
line wrap: on
line diff
--- a/scripts/ode/private/runge_kutta_45_dorpri.m	Wed Oct 14 09:25:04 2015 -0700
+++ b/scripts/ode/private/runge_kutta_45_dorpri.m	Wed Oct 14 10:35:53 2015 -0700
@@ -1,5 +1,5 @@
-## Copyright (C) 2015, Carlo de Falco
-## Copyright (C) 2013, Roberto Porcu' <roberto.porcu@polimi.it>
+## Copyright (C) 2015 Carlo de Falco
+## Copyright (C) 2013 Roberto Porcu' <roberto.porcu@polimi.it>
 ##
 ## This file is part of Octave.
 ##
@@ -23,7 +23,7 @@
 ## @deftypefnx {Function File} {[@var{t_next}, @var{x_next}, @var{x_est}, @var{k_vals_out}] =} runge_kutta_45_dorpri (@var{@@fun}, @var{t}, @var{x}, @var{dt}, @var{options}, @var{k_vals_in})
 ##
 ## This function can be used to integrate a system of ODEs with a given initial
-## condition @var{x} from @var{t} to @var{t+dt}, with the Dormand-Prince method.
+## condition @var{x} from @var{t} to @var{t+dt} with the Dormand-Prince method.
 ## For the definition of this method see
 ## @url{http://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method}.
 ##
@@ -36,7 +36,7 @@
 ## error.
 ##
 ## Fourth output parameter is matrix containing the Runge-Kutta evaluations
-## to use in a FSAL scheme or for dense output.
+## to use in an FSAL scheme or for dense output.
 ##
 ## First input argument is the function describing the system of ODEs to be
 ## integrated.
@@ -52,14 +52,15 @@
 ## adapt the computation to what is needed.
 ##
 ## Sixth input parameter is optional and describes the Runge-Kutta evaluations
-## of the previous step to use in a FSAL scheme.
+## of the previous step to use in an FSAL scheme.
 ## @end deftypefn
 ##
 ## @seealso{odepkg}
 
-function [t_out, x_out, x_est, k] = ...
-         runge_kutta_45_dorpri (f, t, x, dt, opts = [], k_vals = [],
-                                t_out = t + dt)
+function [t_next, x_next, x_est, k] = runge_kutta_45_dorpri (f, t, x, dt,
+                                                             options = [],
+                                                             k_vals = [],
+                                                             t_next = t + dt)
 
   persistent a = [0           0          0           0        0          0;
                   1/5         0          0           0        0          0;
@@ -70,43 +71,43 @@
                   35/384      0          500/1113    125/192 -2187/6784  11/84];
   persistent b = [0 1/5 3/10 4/5 8/9 1 1];
   persistent c = [(35/384) 0 (500/1113) (125/192) (-2187/6784) (11/84)];
-  ## x_est according to Shampine 1986:
+  persistent c_prime = [(5179/57600) 0 (7571/16695) (393/640), ...
+                        (-92097/339200) (187/2100)  (1/40)];
+  ## According to Shampine 1986:
   ## persistent c_prime = [(1951/21600) 0 (22642/50085) (451/720), ...
   ##                       (-12231/42400) (649/6300) (1/60)];
-  persistent c_prime = [(5179/57600) 0 (7571/16695) (393/640), ...
-                        (-92097/339200) (187/2100)  (1/40)];
 
   s = t + dt * b;
   cc = dt * c;
   aa = dt * a;
   k = zeros (rows (x), 7);
 
-  if (! isempty (opts)) # extra arguments for function evaluator
-    args = opts.vfunarguments;
+  if (! isempty (options))  # extra arguments for function evaluator
+    args = options.vfunarguments;
   else
     args = {};
   endif
 
-  if (! isempty (k_vals))  # k values from previous step are passed
+  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(:,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{:});
   k(:,5) = feval (f, s(5), x + k(:,1:4) * aa(5, 1:4).', args{:});
   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;
-  x_out = x + k(:,1:6) * cc(:);  # 5th order approximation
+  ## t_next = t + dt;
+  x_next = 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_out, x_out, args{:});
+    k(:,7) = feval (f, t_next, x_next, args{:});
     cc_prime = dt * c_prime;
     x_est = x + k * cc_prime(:);
   endif