Dense Output

Dense Output

Specifying specific output times for the solution, should not affect the internal time steps that the solver uses. The basic idea of the Dense Output concept is to provide the solution at a given time \(s \in [t, t+dt]\) with the same order of accuracy as the solutions computed at the internal time points by using suitable interpolation methods.
Up to now only linear interpolation was performed and this significantly lowered the accuracy if a higher order solver was used.
I then implemented a series of interpolation function:

• linear_interpolation:
  

  x_out = linear_interpolation (t, x, t_out)

  Given the time span \(t=[t_0, t_1]\) and the function values \(x=[x_0, x_1]\), it returns the linear interpolation value \(x_{out}\) at the point \(t_{out}\).
• quadratic_interpolation:
  

  x_out = quadratic_interpolation (t, x, der, t_out)

  Given the time span \(t=[t_0, t_1]\), the function values \(x=[x_0, x_1]\) and the derivative of the function at the point \(x_0\), it returns the quadratic interpolation value \(x_{out}\) at the point \(t_{out}\).
• hermite_cubic_interpolation:
  

  x_out = hermite_cubic_interpolation (t, x, der, t_out)

  Given the time span \(t=[t_0, t_1]\), the function values \(x=[x_0, x_1]\) and the derivatives of the function at both points \(x_0\) and \(x_1\), it returns the 3rd order approximation \(x_{out}\) at the point \(t_{out}\) by performing Hermite interpolation.
• hermite_quartic_interpolation:
  

  x_out = hermite_quartic_interpolation (t, x, der, t_out)

  Given the time span \(t=[t_0, t_1]\), the function values \(x=[x_0, x_{1/2}, x_1]\) (where \(x_{1/2}\) is the value of the function at the time \(t_0+dt/2\)) and the derivatives of the function at the extremes \(x0\) and \(x1\), it returns the 4th order approximation \(x_{out}\) at the point \(t_{out}\) by performing Hermite interpolation.
• dorpri_interpolation:
  

  x_out = dorpri_interpolation (t, x, k, t_out)

  This interpolation method is specific for the Dormand-Prince Runge-Kutta scheme. Given the time span \(t=[t_0, t_1]\), the function value \(x=x_0\) and the vector \(k\) with the function evaluations required in the Dormand-Prince method, it returns the 4th order approximation \(x_{out}\) at the point \(t_{out}\). For more information on the method have a look at Hairer, Noersett, Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems (pag. 191-193).
• hermite_quintic_interpolation:
  

  x_out = hermite_quintic_interpolation (t, x, der, t_out)

  Given the time span \(t=[t_0, t_1]\), the function values \(x=[x_0, x_{1/2}, x_1]\) and the derivatives of the function at each point, it returns the 5th order approximation \(x_{out}\) at the point \(t_{out}\) by performing Hermite interpolation.

These methods are then used to perform the Dense Output according to the order of the solver chosen. This is the piece of code in integrate_adaptive.m that performs the interpolation:

% if next tspan value is caught, update counter
if( (z(end) == tspan(counter)) || (abs (z(end) - tspan(counter)) / ...
    (max (abs (z(end)), abs (tspan(counter)))) < 8*eps) )
  counter++;

% if there is an element in time vector at which the solution is required
% the program must compute this solution before going on with next steps
elseif( vdirection*z(end) > vdirection*tspan(counter) )
% initializing counter for the following cycle
  i = 2;
  while ( i <= length (z) )

    % if next tspan value is caught, update counter
    if( (counter <= k) && ...
        ( (z(i) == tspan(counter)) || (abs (z(i) - tspan(counter)) / ...
        (max (abs (z(i)), abs (tspan(counter)))) < 8*eps)) )
      counter++;
    endif
    
    % else, loop until there are requested values inside this subinterval
    while((counter <= k) && (vdirection*z(i) > vdirection*tspan(counter)))
      % choose interpolation scheme according to order of the solver
      switch order
      case 1
        u_interp = linear_interpolation ([z(i-1) z(i)], [u(:,i-1) u(:,i)], tspan(counter));
      case 2
        if (~isempty (k_vals))
          der = k_vals(1);
        else
          der = feval (func, z(i-1) , u(:,i-1), options.vfunarguments{:});
        endif
        u_interp = quadratic_interpolation ([z(i-1) z(i)], [u(:,i-1) u(:,i)], der, tspan(counter));
      case 3
        % only ode23 - use k_vals
        u_interp = hermite_cubic_interpolation ([z(i-1) z(i)], [u(:,i-1) u(:,i)], [k_vals(:,1) k_vals(:,end)], tspan(counter));
      case 4
        % if ode45 used without local extrapolation this function doesn't require a new function evaluation.
        u_interp = dorpri_interpolation ([z(i-1) z(i)], [u(:,i-1) u(:,i)], k_vals, tspan(counter));
      case 5
        % ode45 with Dormand-Prince scheme:
        % 4th order approximation of y in t+dt/2 as proposed by Shampine in Lawrence, Shampine, "Some Practical Runge-Kutta Formulas", 1986.
        u_half = u(:,i-1) + 1/2*dt*((6025192743/30085553152)*k_vals(:,1) + (51252292925/65400821598)*k_vals(:,3) - (2691868925/45128329728)*k_vals(:,4) + (187940372067/1594534317056)*k_vals(:,5) - (1776094331/19743644256)*k_vals(:,6) + (11237099/235043384)*k_vals(:,7));
        u_interp = hermite_quartic_interpolation ([z(i-1) z(i)], [u(:,i-1) u_half u(:,i)], [k_vals(:,1) k_vals(:,end)], tspan(counter));

        % it is also possible to do a new function evaluation and the quintic hermite interpolator
        %f_half = feval (func, t+1/2*dt, u_half, options.vfunarguments{:});
        %u_interp = hermite_quintic_interpolation ([z(i-1) z(i)], [u(:,i-1) u_half u(:,i)], [k_vals(:,1) f_half k_vals(:,end)], tspan(counter));
      otherwise
        warning ('high order interpolation not yet implemented: using cubic iterpolation instead');
        der(:,1) = feval (func, z(i-1) , u(:,i-1), options.vfunarguments{:});
        der(:,2) = feval (func, z(i) , u(:,i), options.vfunarguments{:});
        u_interp = hermite_cubic_interpolation ([z(i-1) z(i)], [u(:,i-1) u(:,i)], der, tspan(counter));
      end

      % add the interpolated value of the solution
      u = [u(:,1:i-1), u_interp, u(:,i:end)];
            
      % add the time requested
      z = [z(1:i-1);tspan(counter);z(i:end)];

      % update counters
      counter++;
      i++;
    endwhile

    % if new time requested is not out of this interval
    if ((counter <= k) && (vdirection*z(end) > vdirection*tspan(counter)))
      % update the counter
      i++;
    else
      % else, stop the cycle and go on with the next iteration
      i = length (z) + 1;
    endif

  endwhile

endif

It is important to notice that:

• The 1st order approximation doesn't require any additional function evaluation.
• The 2nd order approximation may require the evaluation of the function at the current time. This can be avoided if the stepper already returns that value.
• The only 3rd order solver implemented is ode23. The 3rd order approximation exploits the Runge-Kutta \(k\) values to avoid further function evaluations.
• There are no 4th order schemes as yet implemented. However if ones were to use ode45 without local extrapolation then the dorpri_interpolation function can be used to obtain a 4th order approximation without any additional function evaluation. For any other 4th order scheme the hermite_quartic_interpolation function can be used.
• For the 5th order method ode45, Shampine proposes to obtain a 4th order approximation at the middle point and to use quartic interpolation. It is however possible to directly do quintic interpolation but this require an additional function evaluation without (according to Shampine) a significant improvement.
• For the higher order solvers (ode78), a suitable interpolator has not yet been implemented.

Finally, since I wrote the interpolation functions in such a way that they are independent of the number of output points requested, a possible improvement would be to compute all the values of the dense output between \(t\) and \(t+dt\) all at once instead of one value at a time.