Re: can ode23s be used with fixed time step?

[c.]

On 23 Feb 2013, at 00:11, Terry Duell <address@hidden> wrote:

> Hello All,
> I would like to use 'ode23s' (latest odepkg), but need to get my solution 
> with constant time steps, to facilitate subsequent analysis of results.
> It appears that 'odeset' only accepts a 'MaxStep' argument.
> Is there a way of using 'ode23s' with fixed time step?
> 
> Cheers,

ode23s uses two embedded rosenbrock methods one of order 2 and the other of 
order 3.
The order 2 method is used for computing the solution while the order 3 method 
is used
to compute an error estimate and adapt the time step.

therefore disabling step-size adaptation defies the whole purpose of using 
ode23s in the 
first place.

Rather than using a fixed stepsize I'd recomend solving with adaptive step and 
then interpolatang on a uniform grid, e.g.:

-----
  vtout = linspace (0, 200, 1e4).';
  fun = @(t,y) [y(2); 100*(1-y(1)^2)*y(2)-y(1)];
  options = odeset ("Jacobian", @(t,y) [0 1; -200*y(1)*y(2)-1, 100*(1-y(1)^2)], 
"InitialStep", 1e-4, "MaxStep", 1);  tic ()
  [vt, vy] = ode23s (fun, [0 200], [2 0], options);
  toc ()
  vyout = interp1 (vt, vy, vtout);
  plot(vt, vy(:, 1), 'ko-', vtout, vyout(:, 1), 'r-');
-----


that said, if you really want to, then you can change  the following lines in 
ode23s: 

-----
    %% 3rd order, needed in error forumula
    F(:,3) = feval (FUN, t+h, x2);
    k(:,3) = W \ (F(:,3) - e32 * (k(:,2)-F(:,2)) - 2 * (k(:,1)-F(:,1)) + h*d*T);

    %% estimate the error
    err = (h/6) * (k(:,1) - 2*k(:,2) + k(:,3));

    %% Estimate the acceptable error
    tau = max (rtol .* abs (x), atol);

    %% Update the solution only if the error is acceptable
    if all (err <= tau)
-----

to

-----
    %% 3rd order, needed in error forumula
    %F(:,3) = feval (FUN, t+h, x2);
    %k(:,3) = W \ (F(:,3) - e32 * (k(:,2)-F(:,2)) - 2 * (k(:,1)-F(:,1)) + 
h*d*T);

    %% estimate the error
    err = 0;%(h/6) * (k(:,1) - 2*k(:,2) + k(:,3));

    %% Estimate the acceptable error
    tau = max (rtol .* abs (x), atol);

    %% Update the solution only if the error is acceptable
    if all (err <= tau)
-----

and then set Initialtep and MaxStep to the same value to maintain a canstant 
step size, e.g.

-----
  fun = @(t,y) [y(2); 100*(1-y(1)^2)*y(2)-y(1)];
  options = odeset ("Jacobian", @(t,y) [0 1; -200*y(1)*y(2)-1, 100*(1-y(1)^2)], 
                  "InitialStep", 1e-2, "MaxStep", 1e-2);
  tic ()
  [vt, vy] = ode23s (fun, [0 200], [2 0], options);
  toc ()
  plot(vt,vy(:,1),'-o');
-----


HTH,
c.