Re: Review ode csets

[Rik]

On 10/19/2016 05:37 PM, c. wrote:
> > > > Rik, > > I see the following problem with this changeset: > > Let myofun be the following: > > function val = myofun (t, y, s) >  disp ('t = ') >  disp (t) >  val = any ((t > 3.2e-2) & (t < .1)); >  disp ('val =') >  disp (val) > end > > and consider the following example: > >>>> myfun  = @(t, y) 1; >>>> [t, y] = ode45 (myfun, [0 1], 1, odeset ('maxstep', 1, 'initialstep', .1, 'refine', 3, 'outputfcn', @myofun)); t >> >> In Matlab this prints: >> ------------------------------------------ >> t = >>     0     1 >> >> val = >>     0 >> >> t = >>    0.0333    0.0667    0.1000 >> >> val = >>     1 >> >> t = >> val = >>     0 >> >> >> t = >> >>         0 >>    0.0333 >>    0.0667 >>    0.1000 >> ------------------------------------------ >> >> So integration stops after the call to myofun with t = [1/30 2/30 3/30] but all the results are returned >> >> In Octave this produces an error: >> >> ------------------------------------------ >> t = >>   0 >>   1 >> val = >> 0 >> t = >> 0 >> val = >> 0 >> error: approxvals(_,2): but approxvals has size 4x1 >> error: called from >>    integrate_adaptive at line 249 column 24 >>    ode45 at line 222 column 12 >> ------------------------------------------ >> >> c. > > > If I change integrate_adaptive to avoid the error: > > ------------------------------------------ > diff --git a/scripts/ode/private/integrate_adaptive.m b/scripts/ode/private/integrate_adaptive.m > --- a/scripts/ode/private/integrate_adaptive.m > +++ b/scripts/ode/private/integrate_adaptive.m > @@ -199,7 +199,7 @@ >              stop_solve = false; >              for ii = 1:numel (approxtime) >                stop_solve = feval (options.OutputFcn, > -                                  approxtime(ii), approxvals(:, ii), [], > +                                  approxtime(ii), approxvals(ii), [], >                                    options.funarguments{:}); >                if (stop_solve) >                  break;  # break from inner loop > @@ -247,7 +247,7 @@ >            stop_solve = false; >            for ii = 1:numel (approxtime) >              stop_solve = feval (options.OutputFcn, > -                                approxtime(ii), approxvals(:, ii), [], > +                                approxtime(ii), approxvals(ii), [], >                                  options.funarguments{:}); >              if (stop_solve) >                break;  # break from inner loop > ------------------------------------------ > > the behaviour of OutpuFcn still looks incompatible > > ------------------------------------------ > t = >    0 >    1 > val = > 0 > t = > 0 > val = > 0 > t = >  0.033333 > val = > 1 > warning: Solver was stopped by a call of 'break' in the main iteration loop at time t = 0.100000 before the endpoint at tend = 1.000000 was reached.  This may happen because the @odeplot function returned 'true' or the @event function returned 'true'. > warning: called from >     integrate_adaptive at line 313 column 7 >     ode45 at line 222 column 12 > t = > [] > val = > 0 > t = > >    0.00000 >    0.10000 > ------------------------------------------ > > data at intermediate times has been thrown away and only the end of the timestep has been returned. > does the above change to integrate_adaptive.m look correct to you?

This turns out to be slightly more complicated.  The problem is that the call to interp1 immediately before the feval code will return different orientations depending on whether the input was a vector (only one variable in y) or a matrix (multiple variables in y).  This is a bit of a hack, but I checked in this cset

changeset:   22649:7458bacc631a
branch:      stable
user:        Rik <address@hidden>
date:        Thu Oct 20 21:14:17 2016 -0700
summary:     integrate_adaptive.m: Fix orientation of approxvals when y is a vector.

diff -r 9bc03a3f7a34 -r 7458bacc631a scripts/ode/private/integrate_adaptive.m
--- a/scripts/ode/private/integrate_adaptive.m  Wed Oct 19 15:17:17 2016 -0700
+++ b/scripts/ode/private/integrate_adaptive.m  Thu Oct 20 21:14:17 2016 -0700
@@ -193,6 +193,9 @@ function solution = integrate_adaptive (
             approxvals = interp1 ([t_old, t(t_caught), t_new],
                                   [x_old, x(:, t_caught), x_new] .',
                                   approxtime, "linear") .';
+            if (isvector (approxvals))
+              approxvals = approxvals.';
+            endif
             if (! isempty (options.OutputSel))
               approxvals = approxvals(options.OutputSel, :);
             endif
@@ -241,6 +244,9 @@ function solution = integrate_adaptive (
           approxvals = interp1 ([t_old, t_new],
                                 [x_old, x_new] .',
                                 approxtime, "linear") .';
+          if (isvector (approxvals))
+            approxvals = approxvals.';
+          endif
           if (! isempty (options.OutputSel))
             approxvals = approxvals(options.OutputSel, :);
           endif

As for the difference in the output variable t, this is because the code which detects a stop_solve event breaks out of the main loop without ever assigning to t.  If you look at the code for the Event function, it does assign to t and x.

              if (! isempty (solution.event{1}) && solution.event{1} == 1)
                t(id) = solution.event{3}(end);
                t = t(1:id);
                x(:, id) = solution.event{4}(end, :).';
                x = x(:,1:id);
                solution.unhandledtermination = false;
                break_loop = true;
                break;
              endif

If you want to get this to work, I think the place to investigate is the exit code for the OutputFcn

          if (stop_solve)  # Leave main loop
            solution.unhandledtermination = false;
            break;
          endif

It seems to me that you would need something like

          if (stop_solve)  # Leave main loop
            t = approxtime;
            x = approxvals;
            solution.unhandledtermination = false;
            break;
          endif

But, you might actually need to append to variables t and x instead.

I fixed the issue I identified earlier where odedefaults did not depend on the input variables to the function due to persistent variables (http://hg.savannah.gnu.org/hgweb/octave/rev/9476c7ddf584).

Should we we make the 4.2.0 release now, or do you want to fix that last compatibility issue?

--Rik