Steppers and Integrate functions

In this first post I want to explain the organization of the code that I'm going to implement. The main idea is to have a structured organization and to subdivide the code according to the most important operations which must be executed. So that should be easier to optimize the bottlenecks and to extend the code with new functionalities.

We will have two main actors: the steppers and the integrate functions.

Steppers

A stepper will be a function and will represent the numerical method used for the integration. Its job is to execute just one integration step and its signature will be:

[x_next,err] = stepper(f,x,t,dt)

x_next is the solution at the next step, err is an estimation of the error (obtainable for example with Richardson Extrapolation), f is a function handle representing the equations to be integrated, x is the solution at the current time t and finally dt is the time step.

Inside this function I'll check if err is really requested by means of nargout keyword. The estimation of the error will be useful to determine the optimal dt in adaptive integrators.

To set the parameters for both the numerical method and the solver I'll use odeset and odeget which are not embedded into Octave but are functions of Odepkg. My intent is to edit these functions in order to make them suitable for my parameters setting.

Integrate Functions

An integrate function will be the function executing the algorithm of integration on more steps. We will have different integrate functions depending on how we want to obtain the solution. For example we can have:

• integrate_const(stepper,f,x0,t0,t1,dt) integrating from t0 to t<=t1 with a fixed dt;
• integrate_n_steps(stepper,f,x0,t0,dt,n) integrating from t0 for n integration steps with fixed dt;
• integrate_adaptive(stepper,p,f,x0,t0,t1,dt) integrating from t0 to t1 with an adaptive timestep, where p is the order of the stepper.

A simple example: Forward Euler

In this part I'll show you an example of two simple steppers (fe_richardson and fe_heun). Both the steppers make use of the Forward Euler method to find the new solution: $$x_{k+1} = x_{k} + h f(t_{k},x_{k})\ .$$ They differ in the way the error is estimated: fe_richardson makes use of Richardson Extrapolation while fe_heun makes use of Heun-Euler method.

This is the implementation of the fe_richardson stepper:

function [t_next,x_next,err] = fwe_richardson(f,t,x,dt,options)
  x_next = x + dt.*f(t,x);
  t_next = t+dt;
  if(nargout == 3)
    x1 = x + (.5*dt).*f(t,x); 
    x_est = x1 + (.5*dt).*f(t+ .5*dt,x1);
    err = norm(x_next - x_est,2);
  end
end

The Heun-Euler method combines the Heun method, which is of order 2, with the Euler method, which is of order 1. This combination allows to get an estimation of the local error. We can write the Heun method as $$\begin{aligned} \tilde{x}_{k+1} &= x_{k} + h f(t,x_{k}) \\ x_{k+1}^{*} &= x_{k} + \frac{h}{2} (f(t,x_{k}) + f(t+h,\tilde{x}_{k+1})) \end{aligned}$$ So that the error estimation is given by: $$ e_{k+1} = \tilde{x}_{k+1} - x^*_{k+1}=\frac{h}{2}(f(t+h,\tilde{x}_{k+1}) - f(t,x_{k}))$$ This is the implementation of the fe_heun method:

function [t_next,x_next,err]=fwe_heun(f,t,x,dt,options)
  x_next = x + dt.*f(t,x);
  t_next = t+dt;
  if(nargout == 3)
    err = .5 * dt * norm(f(t+dt,x_next) - f(t,x),2);
  end
end

The adaptive integrator

Finally I'll describe how integrate_adaptive is implemented. As i said this function does the integration on the total interval [t0,t1] adapting the timestep at each iteration, by means of the error estimation returned by the stepper. This function takes as input:

• stepper: the stepper;
• p: the order of the stepper;
• f: the function to be integrated;
• x0: the initial condition;
• t0, t1: the extremes of the time interval;
• dt: the first timestep.

While the actual time is less or equal to the final instant, this function recursively calls the stepper and then updates the solution, the time and the next timestep. Finally, since the very last instant may be greater than t1, the function substitutes the last element of the solution with a linear interpolation calculated in t1 between the last two elements of the solution.

function [t,x] = integrate_adaptive(stepper,order, ...
                func,tspan,x0,options)
  t = tspan(1);
  x = x0(:);
  dt = odeget(options,'InitialStep',[],'fast');
  tau = 1.e-6;
  counter = 2;
  while( t(end) < (tspan(end)-1.e-13) )
    [s,y,err] = stepper(func,t(end),x(:,end),dt,options);
    if( s(end) > (tspan(counter)+1.e-13) )
      z = [t(end);s];
      u = [x(:,end),y];
      i = 1;
      while( i<=length(z) )
        if( abs(z(i)-tspan(counter)) < 1.e-13 )
          x = [x,y];
          t = [t;s];
          i_old = i;
          i = length(z)+1;
        elseif( z(i) > (tspan(counter)+1.e-13) )
          x = [x(:,1:end-1),u(:,1:i-1),u(:,i-1)+ ...
             ((tspan(counter)-z(i-1))/(z(i)-z(i-1)))* ...
             (u(:,i)-u(:,i-1)),u(:,i:end)];
          t = [t(1:end-1);z(1:i-1);tspan(counter);z(i:end)];
          i_old = i;
          i = length(z)+1;
        end
        i++;
      end
      counter++;
    else
      x = [x,y];
      t = [t;s];
    end
    err += eps;
    dt = dt.*(tau/abs(err))^(1.0/(order+1));
  end
  if( counter > max(size(tspan)) )
    k = (length(z)-i_old)+1;
    x = x(:,1:end-k);
    t = t(1:end-k);
  end
  x = x';
end

The method used to calculate the next timestep is taken from the book of Shampine, Gladwell and Thomson. tau is a parameter which now is fixed but in the future will be setted by means of odeset and odeget functions.

How it will be at the end

At the end, functions like ode45 will be just aliases. For example, we define a function called odefwe which does the integration by means of Forward Euler method, adapting the timestep with Richardson Extrapolation. The call to this function will just hide a call to integrate_adaptive function in this way:

odefwe(f,x0,t0,t1,dt) = ...
    integrate_adaptive(fe_richardson,1,f,x0,t0,t1,dt)

In the same way we'll define a function called ode12 which does the integration in an adaptive way using the Heun-Euler method. So that, now we will have:

ode12(f,x0,t0,t1,dt) = integrate_adaptive(fe_heun,1,f,x0,t0,t1,dt)