Starting Stepsize

Starting Step Size

The implementation of integrate_adaptive function takes (for now) the first timestep as last input argument. Obviously the choice of the first timestep made by the user may be not the best, even though the adaptive technique will adjust the step at the second iteration. So that there exist many techniques to find out a good timestep for the first iteration.

Here I present the one I implemented, which was proposed by Gladwell, Shampine and Brankin in 1987. It is based on the hypothesis that:
$$ \text{local error}\approx Ch^{p+1}y^{(p+1)}(x_0) $$ The algorithm is the following, and usually gives a good guess for the initial timestep.

function h = starting_stepsize(f,x0,t0,p)
  d0 = norm(x0,2);
  d1 = norm(f(t0,x0),2);
  if( d0<1.e-5 || d1<1.e-5 )
    h0 = 1.e-6;
  else
    h0 = .01*(d0/d1);
  end
  x1 = x0 + h0.*f(t0,x0);
  d2 = (1.0/h0)*norm(f(t0+h0,x1)-f(t0,x0),2);
  if( max(d1,d2)<= 1.e-15 )
    h1 = max(1.e-6,h0*1.e-3);
  else
    h1 = (1.e-2 / max(d1,d2))^(1/(p+1));
  end
  h = min(100*h0,h1);
end

So that, what I will do is include this algorithm into integrate_adaptive function which will also have a new signature, because the initial timestep is no more needed:

function [x,t] = integrate_adaptive(stepper,order,f,x0,t0,t1)