Re: Review ode csets

[c.]

On 20 Oct 2016, at 02:20, c. <address@hidden> wrote:

> 
> On 18 Oct 2016, at 20:40, Rik <address@hidden> wrote:
> 
>> 10/18/16
>> 
>> Carlo,
>> 
>> Can you review this cset that I just committed? 
>> http://hg.savannah.gnu.org/hgweb/octave/rev/7efa2d0e22c9.
>> 
>> I changed OutputFcn to behave in a more Matlab-compatible way.  However,
>> part of the change was subtle.  When using the Refine option to generate
>> intermediate timesteps, OutputFcn was called for each of the timesteps. 
>> However, only the final call in the last timestep had the ability to halt
>> the solver by returning a value of TRUE.  This seemed wrong, but maybe
>> there was a reason for it.  The code is:
>> 
>> +            stop_solve = false;
>>            for ii = 1:numel (approxtime)
>> -              pltret = feval (options.OutputFcn, approxtime(ii),
>> -                              approxvals(:, ii), [],
>> -                              options.funarguments{:});
>> +              stop_solve = feval (options.OutputFcn,
>> +                                  approxtime(ii), approxvals(:, ii), [],
>> +                                  options.funarguments{:});
>> +              if (stop_solve)
>> +                break;  # break from inner loop
>> +              endif
>>            endfor
>> -            if (pltret)  # Leave main loop
>> +            if (stop_solve)  # Leave main loop
>> 
>> There is also a patch I submitted to https://savannah.gnu.org/bugs/?49364
>> that needs review.
>> 
>> Thanks,
>> Rik
> 
> 
> 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?

c.