[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Review ode csets
From: |
c. |
Subject: |
Re: Review ode csets |
Date: |
Thu, 20 Oct 2016 02:37:10 +0200 |
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.
- Review ode csets, Rik, 2016/10/18
- Re: Review ode csets, c., 2016/10/19
- Re: Review ode csets, c., 2016/10/19
- Re: Review ode csets, c., 2016/10/19
- Re: Review ode csets,
c. <=
- Re: Review ode csets, Rik, 2016/10/21
- Message not available
- Message not available
- Re: Review ode csets, Carlo de Falco, 2016/10/21
- Re: ode output orientation, Rik, 2016/10/21
- Re: ode output orientation, c., 2016/10/21
- Re: ode output orientation, Sebastian Schöps, 2016/10/21
- Re: ode output orientation, Rik, 2016/10/21
- Re: Review ode csets, c., 2016/10/21
- Re: ode Refine option missing, Rik, 2016/10/21