Mercurial > octave

changeset 22655:6b134d294d61 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
ode solvers: use ordinary transpose instead of Hermitian conjugate (bug #49410). * ode23.m, ode45.m, ode_event_handler.m, runge_kutta_interpolate.m: Use ordinary transpose instead of Hermitian conjugate.
author Carlo de Falco <carlo.defalco@polimi.it>
date Sat, 22 Oct 2016 22:11:42 +0200
parents bc61ed076549
children 56d7d423aff9
files scripts/ode/ode23.m scripts/ode/ode45.m scripts/ode/private/ode_event_handler.m scripts/ode/private/runge_kutta_interpolate.m
diffstat 4 files changed, 22 insertions(+), 10 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/ode/ode23.m	Fri Oct 21 14:36:38 2016 -0700
+++ b/scripts/ode/ode23.m	Sat Oct 22 22:11:42 2016 +0200
@@ -233,7 +233,7 @@
   endif
   if (! isempty (odeopts.Events))   # Cleanup event function handling
     ode_event_handler (odeopts.Events, solution.t(end),
-                       solution.x(end,:)', "done", odeopts.funarguments{:});
+                       solution.x(end,:).', "done", odeopts.funarguments{:});
   endif
 
   ## Print additional information if option Stats is set
@@ -256,8 +256,8 @@
     varargout{1} = solution.t;      # Time stamps are first output argument
     varargout{2} = solution.x;      # Results are second output argument
   elseif (nargout == 1)
-    varargout{1}.x = solution.t';   # Time stamps are saved in field x (row vector)
-    varargout{1}.y = solution.x';   # Results are saved in field y (row vector)
+    varargout{1}.x = solution.t.';   # Time stamps are saved in field x (row vector)
+    varargout{1}.y = solution.x.';   # Results are saved in field y (row vector)
     varargout{1}.solver = solver;   # Solver name is saved in field solver
     if (! isempty (odeopts.Events))
       varargout{1}.ie = solution.event{2};  # Index info which event occurred
@@ -479,6 +479,12 @@
 ## test for MaxOrder option is missing
 ## test for MvPattern option is missing
 
+%!test # Check that imaginary part of solution does not get inverted
+%! sol = ode23 (@(x,y) 1, [0 1], 1i);
+%! assert (imag (sol.y), ones (size (sol.y)))
+%! [x, y] = ode23 (@(x,y) 1, [0 1], 1i);
+%! assert (imag (y), ones (size (y)))
+
 ## Test input validation
 %!error ode23 ()
 %!error ode23 (1)
--- a/scripts/ode/ode45.m	Fri Oct 21 14:36:38 2016 -0700
+++ b/scripts/ode/ode45.m	Sat Oct 22 22:11:42 2016 +0200
@@ -224,7 +224,7 @@
   endif
   if (! isempty (odeopts.Events))   # Cleanup event function handling
     ode_event_handler (odeopts.Events, solution.t(end),
-                       solution.x(end,:)', "done", odeopts.funarguments{:});
+                       solution.x(end,:).', "done", odeopts.funarguments{:});
   endif
 
   ## Print additional information if option Stats is set
@@ -247,8 +247,8 @@
     varargout{1} = solution.t;      # Time stamps are first output argument
     varargout{2} = solution.x;      # Results are second output argument
   elseif (nargout == 1)
-    varargout{1}.x = solution.t';   # Time stamps are saved in field x (row vector)
-    varargout{1}.y = solution.x';   # Results are saved in field y (row vector)
+    varargout{1}.x = solution.t.';   # Time stamps are saved in field x (row vector)
+    varargout{1}.y = solution.x.';   # Results are saved in field y (row vector)
     varargout{1}.solver = solver;   # Solver name is saved in field solver
     if (! isempty (odeopts.Events))
       varargout{1}.ie = solution.event{2};  # Index info which event occurred
@@ -487,6 +487,12 @@
 ## test for MaxOrder option is missing
 ## test for MvPattern option is missing
 
+%!test # Check that imaginary part of solution does not get inverted
+%! sol = ode45 (@(x,y) 1, [0 1], 1i);
+%! assert (imag (sol.y), ones (size (sol.y)))
+%! [x, y] = ode45 (@(x,y) 1, [0 1], 1i);
+%! assert (imag (y), ones (size (y)))
+
 %!error ode45 ()
 %!error ode45 (1)
 %!error ode45 (1,2)
--- a/scripts/ode/private/ode_event_handler.m	Fri Oct 21 14:36:38 2016 -0700
+++ b/scripts/ode/private/ode_event_handler.m	Sat Oct 22 22:11:42 2016 +0200
@@ -88,7 +88,7 @@
     [evt, term, dir] = feval (inpargs{:});
 
     ## We require that all return values be row vectors
-    evt = evt(:)'; term = term(:)'; dir = dir(:)';
+    evt = evt(:).'; term = term(:).'; dir = dir(:).';
 
     ## Check if one or more signs of the event has changed
     signum = (sign (evtold) != sign (evt));
@@ -115,7 +115,7 @@
         ## calculate new values for the integration results, we do both by
         ## a linear interpolation
         tnew = t - evt(1,idx) * (t - told) / (evt(1,idx) - evtold(1,idx));
-        ynew = (y - (t - tnew) * (y - yold) / (t - told))';
+        ynew = (y - (t - tnew) * (y - yold) / (t - told)).';
         retcell{3}(evtcnt,1) = tnew;
         retcell{4}(evtcnt,:) = ynew;
         evtcnt += 1;
@@ -142,7 +142,7 @@
     [evtold, ~, ~] = feval (inpargs{:});
 
     ## We require that all return values be row vectors
-    evtold = evtold(:)'; told = t; yold = y;
+    evtold = evtold(:).'; told = t; yold = y;
     evtcnt = 1;
     retval = retcell = cell (1,4);
 
--- a/scripts/ode/private/runge_kutta_interpolate.m	Fri Oct 21 14:36:38 2016 -0700
+++ b/scripts/ode/private/runge_kutta_interpolate.m	Sat Oct 22 22:11:42 2016 +0200
@@ -22,7 +22,7 @@
   switch (order)
 
     case 1
-      u_interp = interp1 (z, u', t, "linear");
+      u_interp = interp1 (z, u.', t, "linear");
 
     case 2
       if (! isempty (k_vals))