Mercurial > octave-nkf

diff scripts/ode/private/integrate_adaptive.m @ 20584:eb9e2d187ed2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Use Octave coding conventions in scripts/ode/private dir. * AbsRel_Norm.m, fuzzy_compare.m, hermite_quartic_interpolation.m, integrate_adaptive.m, integrate_const.m, integrate_n_steps.m, kahan.m, ode_struct_value_check.m, odepkg_event_handle.m, odepkg_structure_check.m, runge_kutta_45_dorpri.m, starting_stepsize.m: Wrap long lines to < 80 chars. Use double quotes rather than single quotes where possible. Use ';' at end of keywords "return;" and "break;" Use '##" for stand-alone comments and '#' for end-of-line comments. Use two spaces after period before starting new sentence. Use '!' instead of '~' for logical negation. Use specific form of end (endif, endfor, etc.). Don't use line continuation marker '...' unless necessary.
author Rik <rik@octave.org>
date Sun, 04 Oct 2015 22:18:54 -0700
parents 25623ef2ff4f
children b7ac1e94266e
line wrap: on
line diff
--- a/scripts/ode/private/integrate_adaptive.m	Sun Oct 04 16:24:32 2015 +0100
+++ b/scripts/ode/private/integrate_adaptive.m	Sun Oct 04 22:18:54 2015 -0700
@@ -61,7 +61,7 @@
 
 function solution = integrate_adaptive (stepper, order, func, tspan, x0, options)
 
-  solution = struct;
+  solution = struct ();
 
   ## first values for time and solution
   t = tspan(1);
@@ -78,16 +78,16 @@
   endif
   dt = vdirection * min (abs (dt), options.MaxStep);
 
-  ## set parameters
+  ## Set parameters
   k = length (tspan);
   counter = 2;
   comp = 0.0;
   tk = tspan(1);
   options.comp = comp;
   
-  ## factor multiplying the stepsize guess
+  ## Factor multiplying the stepsize guess
   facmin = 0.8;
-  fac = 0.38^(1/(order+1)); ## formula taken from Hairer
+  fac = 0.38^(1/(order+1));  # formula taken from Hairer
   t_caught = false;
 
 
@@ -122,7 +122,7 @@
   while (counter <= k)
     facmax = 1.5;
 
-    ## compute integration step from t to t+dt
+    ## Compute integration step from t to t+dt
     if (isempty (k_vals))
       [s, y, y_est, new_k_vals] = stepper (func, z(end), u(:,end),
                                            dt, options);
@@ -144,7 +144,7 @@
     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
+    ## Solution accepted only if the error is less or equal to 1.0
     if (err <= 1)
       
       [tk, comp] = kahan (tk, comp, dt);
@@ -162,8 +162,9 @@
               (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
@@ -233,11 +234,12 @@
                 ## 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)],
+                ##                                  [k_vals(:,1) f_half ...
+                ##                                   k_vals(:,end)],
                 ##                                  tspan(counter));
               otherwise
                 warning ("High order interpolation not yet implemented: ",
-                         "using cubic iterpolation instead");
+                         "using cubic interpolation instead");
                 der(:,1) = feval (func, z(i-1) , u(:,i-1),
                                   options.vfunarguments{:});
                 der(:,2) = feval (func, z(i) , u(:,i),
@@ -281,7 +283,7 @@
       vcntiter = 0;
       
       ## Call plot only if a valid result has been found, therefore this
-      ## code fragment has moved here. Stop integration if plot function
+      ## code fragment has moved here.  Stop integration if plot function
       ## returns false
       if (options.vhaveoutputfunction)
         for vcnt = 0:options.Refine # Approximation between told and t
@@ -300,27 +302,28 @@
           vpltret = feval (options.OutputFcn, vapproxtime,
                            vapproxvals, [], options.vfunarguments{:});
           if (vpltret) # Leave refinement loop
-            break
+            break;
           endif
         endfor
         if (vpltret) # Leave main loop
           solution.vunhandledtermination = false;
-          break
+          break;
         endif
       endif
       
       ## Call event only if a valid result has been found, therefore this
-      ## code fragment has moved here. Stop integration if veventbreak is
+      ## code fragment has moved here.  Stop integration if veventbreak is
       ## 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,:);
           x(:,solution.vcntloop-1) = solution.vevent{4}(end,:)';
           solution.vunhandledtermination = false; 
-          break
+          break;
         endif
       endif
       
@@ -337,18 +340,18 @@
     dt = vdirection * min (abs (dt), options.MaxStep);
     
     ## Update counters that count the number of iteration cycles
-    solution.vcntcycles = solution.vcntcycles + 1; # Needed for cost statistics
-    vcntiter = vcntiter + 1; # Needed to find iteration problems
+    solution.vcntcycles += 1; # Needed for cost statistics
+    vcntiter += 1; # Needed to find iteration problems
 
     ## Stop solving because in the last 1000 steps no successful valid
     ## value has been found
     if (vcntiter >= 5000)
-      error (["Solving has not been successful. The iterative",
+      error (["Solving has not been successful.  The iterative",
               " integration loop exited at time t = %f before endpoint at",
-              " tend = %f was reached. This happened because the iterative",
+              " tend = %f was reached.  This happened because the iterative",
               " integration loop does not find a valid solution at this time",
-              " stamp. Try to reduce the value of ''InitialStep'' and/or",
-              " ''MaxStep'' with the command ''odeset''.\n"],
+              " stamp.  Try to reduce the value of 'InitialStep' and/or",
+              " 'MaxStep' with the command 'odeset'.\n"],
              s(end), tspan(end));
     endif
 
@@ -362,20 +365,20 @@
   if (vdirection * z(end) < vdirection * tspan(end))
     if (solution.vunhandledtermination == true)
       error ("OdePkg:InvalidArgument",
-             ["Solving has not been successful. The iterative",
+             ["Solving has not been successful.  The iterative",
               " integration loop exited at time t = %f",
-              " before endpoint at tend = %f was reached. This may",
+              " before endpoint at tend = %f was reached.  This may",
               " happen if the stepsize grows smaller than defined in",
-              " vminstepsize. Try to reduce the value of ''InitialStep''",
-              " and/or ''MaxStep'' with the command ''odeset''.\n"],
+              " vminstepsize.  Try to reduce the value of 'InitialStep'",
+              " and/or 'MaxStep' with the command 'odeset'.\n"],
              z(end), tspan(end));
     else
       warning ("OdePkg:InvalidArgument",
-               ["Solver has been stopped by a call of ''break'' in the main",
+               ["Solver has been stopped by a call of 'break' in the main",
                 " iteration loop at time t = %f before endpoint at tend = %f ",
-                " was reached. This may happen because the @odeplot function",
-                " returned ''true'' or the @event function returned",
-                " ''true''.\n"],
+                " was reached.  This may happen because the @odeplot function",
+                " returned 'true' or the @event function returned",
+                " 'true'.\n"],
                z(end), tspan(end));
     endif
   endif
@@ -389,3 +392,4 @@
   solution.x = x(:,1:end-f)';
   
 endfunction
+