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.

FSAL - new stepper implementation

As stated in the previous post, the implementation of the steppers as it was did not allow the possibility to exploit the FSAL (First Same As Last) property of the Bogacki-Shampine algorithm (ode23) and of the Dormand-Prince algorithm (ode45).
The input and output arguments of the steppers have then be modified. As an example here is the runge_kutta_23 stepper:

function varargout = runge_kutta_23 (f, t, x, dt, varargin)
  options = varargin{1};
  k = zeros (size (x, 1), 4);

  if (nargin == 5) % only the options are passed
    k(:,1) = feval (f, t , x, options.vfunarguments{:});
  elseif (nargin == 6) % both the options and the k values are passed
    k(:,1) = varargin{2}(:,end); % FSAL property
  endif
  k(:,2) = feval (f, t + (1/2)*dt, x + dt*(1/2)*k(:,1), options.vfunarguments{:});
  k(:,3) = feval (f, t + (3/4)*dt, x + dt*(3/4)*k(:,2), options.vfunarguments{:});

  %# computing new time and new values for the unkwnowns
  varargout{1} = t + dt; %t_next
  varargout{2} = x + dt.*((2/9)*k(:,1) + (1/3)*k(:,2) + (4/9)*k(:,3)); % return the 3rd order approximation x_next

  %# if the estimation of the error is required
  if (nargout >= 3)
    %# new solution to be compared with the previous one
    k(:,4) = feval (f, t + dt, varargout{2}, options.vfunarguments{:});
    varargout{3} = x + dt.*((7/24)*k(:,1) + (1/4)*k(:,2) + (1/3)*k(:,3) + (1/8)*k(:,4)); %x_est
    varargout{4} = k;
  endif

endfunction

And the call within the solver becomes:

[s, y, y_est, k_vals] = stepper (func, z(end), u(:,end), dt, options, k_vals);

where k_vals has to be initialized for the first iteration as f(t, x).
This implementation will reduce the number of function evaluation for each step.

Furthermore, after some tests in MATLAB, the return values for the solution and the estimate in the  runge_kutta_23 and runge_kutta_45 steppers have been swapped to automatically perform local extrapolation. The MATLAB functions are in fact of order 3 and 5 respectively.

Status of the code: bugfixes and new issues

ODESET and ODEGET

• odeset and odeget functions have been slightly modified to be compliant with MATLAB. Each MATLAB option is present and all the options are tested. The coding style has been adapted to the GNU-Octave standard.
• ode_struct_value_check: this function has been introduced by Roberto in addition to odepkg_structue_check. The relation between the two functions has to be clarified: in particular it is necessary to understand if it is really necessary to have two different functions or one is sufficient.

CHANGES TO THE STEPPERS

• The runge_kutta_78 stepper has been implemented.
• Two 4th order steppers have been implemented: runge_kutta_45_dopri (Dormand-Prince coefficients) and runge_kutta_45_fehlberg (Fehlberg coefficients).

CHANGES TO THE SOLVERS

• ode78 solver has been updated to the new structure. It now exploits the runge_kutta_78 stepper.
• A series of tests has been added to each solver to check all the functionalities and the all options. This has made me possible to detect some bugs that have been corrected. In particular the adaptive timestep evaluation had some issues that lead to the use of too short timesteps. This has been corrected and now the algorithm proposed in [1] is used.
• Furthermore the current implementation uses linear interpolation to evaluate the solution at the user specified times. This leads to a considerable loss in accuracy and is not consistent with MATLAB (which guarantees the same order of accuracy of the scheme used). In [1] different methods are proposed for the dense output: these will be used as a reference for the implementation of a better solution.
• In the released version of odepkg some of the solvers perform local extrapolation, that is the higher-order estimate is chosen as the solution. With the new stepper structure, as it is now, this choice effects all the solvers. It have to be decided whether to perform it or not (MATLAB doesn't seem to do it, thus I suggest to avoid it).
• MATLAB implementation of ode45 uses the Dormand-Prince (DP) coefficients. In the released version of odepkg there exits two solvers: ode45 that uses the Fehlberg coefficients and ode54 that uses the DP coefficients. To be consistent with MATLAB, ode45 now uses the DP method. This makes the runge_kutta_45_fehlberg stepper and the ode54 solver, as it is now, redundant. Either their elimination or a change of the solver might be considered. However one of the advantages of DP coefficients is the possibility to reuse the last function evaluation at a given step as the first evaluation of the subsequent one. This is not easily done with the stepper structure introduced by Roberto.

CHANGES TO THE OPTIONS

• InitialStep option has been modified to be MATLAB compatible (it must be a positive scalar).
• RelTol defalut value has been changed to 1e-3 instead of 1e-6 to be MATLAB compatible.
• MaxStep option has been implemented.
• NormControl option has been implemented.

TODO

In addition to the general plan, a couple of issues need to be addressed:

• Clarify the relation between ode_struct_value_check and odepkg_structue_check.
• Decide if local extrapolation has to be used or not. My opinion (and the current implementation) is to avoid it to be compliant to what MATLAB seems to be doing.
• Solve the dense output problem in a way that guarantees the consistency with MATLAB.
• Consider if it's possible to reduce the number of function evaluation for the Dormand-Prince stepper (ode45) and the Bogacki-Shampine stepper (ode23) exploiting the FSAL property (first same as last).
• Decide if in the future releases of odepkg ode54 has to be removed or maybe changed to become a Fehlberg solver.

[1] E. Hairer, S.P. N{\o}rsett, G. Wanner, Solving Ordinary Differential Equations, 1993, Springer.

SoCiS 2014 - New timeline

On the occasion of SoCiS 2014 I will take over the work that has been done last year on odepkg and continue it. The final goal is to release a new stable version of odepkg and to insert the most common solver in core-Octave in such a way that everything is MATLAB-compatible.

The following list explains the main points of the timeline for my SOCIS project:

1. Check of the current status of the code, in particular with respect to the current release on SourceForge. The two repository will be merged so that every test present in the old version will be added to the new code. Verify if there are missing features and add them if necessary.
2. Comparison of the performance between the old and the new structure. In particular we expect that the introduction of the Levenshtein algorithm for the string comparisons will be a critical issue. If necessary implement levenshtein.m and fuzzy_compare.m in C++.
3. Verify that the functions odeset and odeget are MATLAB-compatible and compliant to Octave core. Add the two functions to the core.
4. Move ode45, ode23 and ode23s to Octave core.
5. Implement ode15s solver. This solver is still missing in odepkg but is highly suitable for stiff problems.
6. Move ode15s to Octave core.
7. New release of odepkg.

Documentation

This is the link for the documentation file :
SOCIS 2013 - Documentation for new functions added to odepkg

This is the link for the implemented code :
OdePkg with my added functions

RATTLE algorithm

RATTLE

Constrained mechanical systems form an important class of differential equations on manifolds. For all the theory that I present here and I've used for the implementation I refer to "Geometric Numerical Integration" by Hairer, Lubich and Wanner.

Consider a mechanical system described by position coordinates \(q_1,\dots q_d\) and suppose that the motion is constrained to satisfy \(g(\mathbf q) = \mathbf 0\) where \(g:\mathbb R^d \rightarrow \mathbb R^m \), with \(m \lt d \). So that, the equations of motion governing the system become: $$ \left\{ \begin{array}{l} \dot{\mathbf q} = \frac{\partial H}{\partial \mathbf p} \\ \dot{\mathbf p} = -\frac{\partial H}{\partial \mathbf q} - G(q)^T\lambda \\ g(\mathbf q) = \mathbf 0 \end{array}\right. $$ where \(G(\mathbf q) = \dfrac{\partial g(\mathbf q)}{\partial \mathbf q} \).

Symplectic Euler method can be extended to constrained systems but we focus on SHAKE and RATTLE algorithms. SHAKE is a 2-steps algorithm so that, since I'm implementing only 1-step algorithms and the overall structure of solvers and integrators is made for 1-step solvers, I implemented just RATTLE algorithm.

odeRATTLE

The RATTLE algorithm implemented works with any general Hamiltonian \(H(\mathbf q,\mathbf p \) and is defined as follows: $$ \left\{\begin{array}{l} \mathbf p_{n+\frac{1}{2}} = \mathbf p_n -\frac{h}{2}\big(H_{\mathbf q}(\mathbf q_n,\mathbf p_{n+\frac{1}{2}}) + G(\mathbf q_n)^T \mathbf {\lambda}_n \big) \\ \mathbf q_{n+1} = \mathbf q_n +\frac{h}{2} \big( H_{\mathbf p}(\mathbf q_n,\mathbf p_{n+\frac{1}{2}}) + H_{\mathbf p}(\mathbf q_{n+1},\mathbf p_{n+\frac{1}{2}}) \big) \\ g(\mathbf q_{n+1}) = \mathbf 0 \\ \mathbf p_{n+1} = \mathbf p_{n+\frac{1}{2}} -\frac{h}{2}\big(H_{\mathbf q}(\mathbf q_{n+1},\mathbf p_{n+\frac{1}{2}}) + G(\mathbf q_{n+1})^T \mathbf{\mu}_n \big) \\ G(\mathbf q_{n+1}) H_{\mathbf p}(\mathbf q_{n+1},\mathbf p_{n+1}) = \mathbf 0 \end{array}\right. $$ where \( h=\Delta t=t_{k+1}-t_k\) and \(\mathbf{\mu}_n \),\( \mathbf{\mu}_n \) are Lagrangian multipliers nedded to impose the constraints.

It can be demonstrated that this numerical method is symmetric, symplectic and convergent of order 2.

The following code represent it's implementation:

function [t_next,x_next,err]=rattle(f,t,x,dt,options)
  H = odeget(options,'HamiltonianHessFcn',[],'fast');
  GG = odeget(options,'ConstraintHessFcn',[],'fast');
  if( ~isempty(H) && ~isempty(GG) )
    fsolve_opts = optimset('Jacobian','on');
  else
    fsolve_opts = optimset('Jacobian','off');
  end
  g = odeget(options,'ConstraintFcn',[],'fast');
  G = odeget(options,'ConstraintGradFcn',[],'fast');
  c_nb = odeget(options,'ConstraintsNb',[],'fast');

  dim = length(x)/2;
  q0 = x(1:dim);
  p0 = x(dim+1:end);

  RATTLE = @(y)constr_sys(y,dim,c_nb,f,H,g,G,GG,t,q0,p0,dt);
  y0 = [q0;p0;p0;zeros(2*c_nb,1)];
  y0 = fsolve(RATTLE,y0,fsolve_opts);
  t_next = t+dt;
  x_next = [y0(1:dim);y0(2*dim+1:3*dim)];

  if(nargout==3)
    dt = dt/2;
    RATTLE = @(y)constr_sys(y,dim,c_nb,f,H,g,G,GG,t,q0,p0,dt);
    y0 = [q0;p0;p0;zeros(2*c_nb,1)];
    y0 = fsolve(RATTLE,y0,fsolve_opts);

    q0 = y0(1:dim);
    p0 = y0(2*dim+1:3*dim);
    t = t+dt;

    RATTLE = @(y)constr_sys(y,dim,c_nb,f,H,g,G,GG,t,q0,p0,dt);
    y0 = [q0;p0;p0;zeros(2*c_nb,1)];
    y0 = fsolve(RATTLE,y0,fsolve_opts);

    x_est = [y0(1:dim);y0(2*dim+1:3*dim)];
    err = norm(x_est-x_next,2);
  end
end

function [F,J] = constr_sys(y,dim,c_nb,f,H,g,G,GG,t,q0,p0,dt)
  F = zeros(3*dim+2*c_nb,1);
  F(1:dim) = y(1:dim) - q0 - (dt/2).*(f(t,[q0; ...
        y(dim+1:2*dim)])(1:dim) + f(t+dt,y(1:2*dim))(1:dim));
  F(dim+1:2*dim) = y(dim+1:2*dim) - p0 - (dt/2).*(f(t,[q0; ...
   y(dim+1:2*dim)])(dim+1:end) - G(q0)'*y(3*dim+1:3*dim+c_nb));
  F(2*dim+1:3*dim) = y(2*dim+1:3*dim) - y(dim+1:2*dim) - ...
             (dt/2)*(f(t+dt,y(1:2*dim))(dim+1:end) - ...
             G(y(1:dim))'*y(3*dim+c_nb+1:end));
  F(3*dim+1:3*dim+c_nb) = g(y(1:dim));
  F(3*dim+c_nb+1:end) = G(y(1:dim))*(f(t+dt,[y(1:dim); ...
                                 y(2*dim+1:3*dim)])(1:dim));

  if( nargout==2 )
    J = zeros(3*dim+2*c_nb,3*dim+2*c_nb);
    J(1:dim,1:dim) = eye(dim) - ...
              (dt/2)*(H(t+dt,y(1:2*dim))(dim+1:end,1:dim));
    J(1:dim,dim+1:2*dim) = -(dt/2)*(H(t,[q0; ...
          y(dim+1:2*dim)])(dim+1:end,dim+1:end) + ...
            H(t+dt,y(1:2*dim))(dim+1:end,dim+1:end));
    J(dim+1:2*dim,dim+1:2*dim) = eye(dim) + ...
        (dt/2)*(H(t,[q0;y(dim+1:2*dim)])(1:dim,dim+1:end));
    J(dim+1:2*dim,3*dim+1:3*dim+c_nb) = (dt/2)*G(q0)';
    J(2*dim+1:3*dim,1:dim) = (dt/2)*(H(t+dt, ...
                            y(1:2*dim))(1:dim,1:dim));
    for k = 1:1:c_nb
      J(2*dim+1:3*dim,1:dim) = J(2*dim+1:3*dim,1:dim) - ...
            (dt/2)*(y(3*dim+c_nb+k)*(GG(y(1:dim))(:,:,k)));
    end
    J(2*dim+1:3*dim,dim+1:2*dim) = -eye(dim) + ...
               (dt/2)*(H(t+dt,y(1:2*dim))(1:dim,dim+1:end));
    J(2*dim+1:3*dim,2*dim+1:3*dim) = eye(dim) + ...
               (dt/2)*(H(t+dt,y(1:2*dim))(1:dim,dim+1:end));
    J(2*dim+1:3*dim,3*dim+c_nb+1:end) = (dt/2)*G(y(1:dim))';
    J(3*dim+1:3*dim+c_nb,1:dim) = G(y(1:dim));
    J(3*dim+c_nb+1:end,1:dim) = G(y(1:dim))* ...
      (H(t+dt,[y(1:dim);y(2*dim+1:3*dim)])(dim+1:end,1:dim));
    for k = 1:1:c_nb
      J(3*dim+c_nb+k,1:dim) = J(3*dim+c_nb+k,1:dim) + ...
          ((GG(y(1:dim))(:,:,k))*(f(t+dt,[y(1:dim); ...
                     y(2*dim+1:3*dim)])(1:dim)))';
    end
    J(3*dim+c_nb+1:end,2*dim+1:3*dim) = G(y(1:dim))* ...
        (H(t+dt,[y(1:dim);y(2*dim+1:3*dim)]) ...
                  (dim+1:end,dim+1:end));
  end
end

It works with any number of constraint, unless this is equal or greater to system dimension. As usual, all the source code is available at my public repository octave-odepkg.

Spectral Variational Integrators

Variational integrators

The variational integrators are a class of numerical methods for mechanical systems which comes from the discrete formulation of Hamilton's principle of stationary action. They can be applied to ODEs, PDEs and to both conservative and forced systems. In absence of forcing terms these methods preserve momenta related to symmetries of the problem and don't dissipate energy. So that they exhibit long-term stability and good long-term behaviour.

Considering a configuration manifold V, the discrete Lagrangian is a function from V to the real numbers space, which represents an approximation of the action between two fixed configurations: $$ L_d(\mathbf q_k,\mathbf q_{k+1}) \approx \int_{t_k}^{t_{k+1}} L(\mathbf q,\dot{\mathbf q};t) dt \hspace{0.5cm}\text{with}\hspace{0.4cm}\mathbf q_{k},\mathbf q_{k+1}\hspace{0.3cm}\text{fixed.}$$ From here, applying the Hamilton's principle, we can get the Euler-Lagrange discrete equations: $$D_2L_d(\mathbf q_{k-1},\mathbf q_k) + D_1 L_d(\mathbf q_k,\mathbf q_{k+1}) = 0\ ,$$ and thanks to the discrete Legendre transforms we get the discrete Hamilton's equation of motion: $$\left\{\begin{array}{l} \mathbf p_k = -D_1 L_d(\mathbf q_k,\mathbf q_{k+1}) \\ \mathbf p_{k+1} = D_2 L_d(\mathbf q_k,\mathbf q_{k+1}) \end{array}\right.\ ,$$ so that we pass from a 2-step system of order N to a 1-step system of order 2N. This system gives the updating map: $$ (\mathbf q_k,\mathbf p_k)\rightarrow(\mathbf q_{k+1},\mathbf p_{k+1})\ .$$ For all the theory behind this I refer to "Discrete mechanics and variational integrators" by Marsden and West.

Spectral variational integrators

To create a spectral variational integrator I considered a discretization of the configuration manifold on a n-dimensional functional space generated by the orthogonal basis of Legendre polynomials. So that, after rescaling the problem from \([t_k,t_{k+1}]\) (with \(t_{k+1}-t_k=h\)) to \([-1,1]\), we get: $$ \begin{array}{l} \mathbf q_k(z) = \sum_{i=0}^{n-1}\mathbf q_k^i l_i(z)\\ \dot{\mathbf q}_k(z) = \frac{2}{h} \sum_{i=0}^{n-1}\mathbf q_k^i \dot l_i(z) \end{array} \ .$$ Then I approximate the action using the Gaussian quadrature rule using \(m\) internal nodes, so putting all together we have: $$ \int_{t_k}^{t_{k+1}} L(\mathbf q,\dot{\mathbf q};t)dt\hspace{0.5cm} \approx\hspace{0.5cm} \frac{h}{2}\sum_{j=0}^{m-1} \omega_j L\big( \sum_{i=0}^{n-1}\mathbf q_k^i l_i , \frac{2}{h} \sum_{i=0}^{n-1}\mathbf q_k^i \dot l_i \big) $$ Now imposing Hamilton's principle of stationary action under the constraints: $$ \mathbf q_k = \sum_{i=0}^{n-1}\mathbf q_k^i l_i(-1) \hspace{1.5cm} \mathbf q_{k+1} = \sum_{i=0}^{n-1}\mathbf q_k^i l_i(1)\ ,$$ we obtain the system: $$ \left\{ \begin{array}{l} \sum_{j=0}^{m-1}\omega_j \bigg[ \mathbf p_j \dot{l}_s(z_j) - \frac{h}{2} l_s(z_j) \frac{\partial H}{\partial \mathbf q} \bigg ( \sum_{i=0}^{n-1}\mathbf q_k^i l_i(z_j),\mathbf p_j \bigg ) \bigg ] + l_s(-1)\mathbf p_k - l_s(1)\mathbf p_{k+1} = 0 \hspace{1cm} \forall s=0,\dots,n-1\\ \frac{\partial H}{\partial \mathbf p}\bigg (\sum_{i=0}^{n-1}\mathbf q_k^i l_i(z_r),\mathbf p_r\bigg ) -\frac{2}{h}\ \sum_{i=0}^{n-1} \mathbf q_k^i \dot{l}_i(z_r)=0 \hspace{1cm} \forall r=0,\dots,m-1 \\ \sum_{i=0}^{n-1} \mathbf q_k^i l_i(-1) - \mathbf q_k = 0\\ \sum_{i=0}^{n-1} \mathbf q_k^i l_i(1) - \mathbf q_{k+1} = 0 \end{array}\right. $$

odeSPVI

Within odeSPVI I implemented a geometric integrator for Hamiltonian systems, like odeSE and odeSV, which uses spectral variational integrators with any order for polynomials of the basis and with any number of internal quadrature nodes, for both unidimensional and multidimensional problems.

[T,Y]=odeSPVI(@hamilton_equations,time,initial_conditions,options)

This solver just uses the stepper spectral_var_int.m which is the function where I implemented the resolution of the system displayed above; hamilton_equations must be a function_handle or an inline function where the user has to implement Hamilton's equations of motion: $$ \left\{\begin{array}{l} \dot{\mathbf q} = \frac{\partial H}{\partial \mathbf p}(\mathbf q,\mathbf p)\\ \dot{\mathbf p} = - \frac{\partial H}{\partial \mathbf q}(\mathbf q,\mathbf p) + \mathbf f(\mathbf q,\mathbf p)\end{array}\right.$$ where \(H(\mathbf q,\mathbf p)\) is the Hamiltonian of the system and \(\mathbf f(\mathbf q,\mathbf p)\) is an optional forcing term.

This below is the implementation of the stepper:

function [t_next,x_next,err]=spectral_var_int(f,t,x,dt,options)
  fsolve_opts = optimset('Jacobian','on');

  q_dofs = odeget(options,'Q_DoFs',[],'fast');
  p_dofs = odeget(options,'P_DoFs',[],'fast');
  nodes = odeget(options,'Nodes',[],'fast');
  weights = odeget(options,'Weights',[],'fast');
  leg = odeget(options,'Legendre',[],'fast');
  deriv = odeget(options,'Derivatives',[],'fast');
  extremes = odeget(options,'Extremes',[],'fast');

  dim = length(x)/2;
  N = q_dofs+p_dofs+2;
  q0 = x(1:dim)(:)';
  p0 = x(dim+1:end)(:)';

  SPVI = @(y)svi_system(y,q_dofs,p_dofs,weights,leg, ...
                deriv,extreme,dim,N,f,t,q0,p0,dt);
  y0 = [q0;zeros(q_dofs-1,dim);ones(p_dofs,1)*p0;q0;p0];
  y0 = reshape(y0,dim*N,1);
  z = fsolve(SPVI,y0,fsolve_opts);
  %z = inexact_newton(SVI,y0,new_opts) %new_opts must be set 
  z = reshape(z,N,dim);

  y = [leg*z(1:q_dofs,1:dim),z((q_dofs+1):(end-2),1:dim)];
  x_next = [y;z(end-1,:),z(end,:)]';
  t_next = [t+(dt/2)*(nodes+1), t+dt]';

  if(nargout==3)
    q_dofs_err = odeget(options,'Q_DoFs_err',[],'fast');
    p_dofs_err = odeget(options,'P_DoFs_err',[],'fast');
    nodes_err = odeget(options,'Nodes_err',[],'fast');
    weights_err = odeget(options,'Weights_err',[],'fast');
    leg_err = odeget(options,'Legendre_err',[],'fast');
    deriv_err = odeget(options,'Derivatives_err',[],'fast');
    extreme_err = odeget(options,'Extremes_err',[],'fast');

    N_err = q_dofs_err+p_dofs_err+2;
    q1 = x_next(1:dim,end)(:)';
    p1 = x_next(dim+1:end,end)(:)';

    SPVI = @(y)svi_system(y,q_dofs_err,p_dofs_err, ...
         weights_err,leg_err,deriv_err,extreme_err, ...
         dim,N_err,f,t,q0,p0,dt);
    p_interp = zeros(p_dofs_err,dim);

    p_lower_order = [p0;z(q_dofs+1:q_dofs+p_dofs,:);p1];
    for i=1:1:p_dofs_err
      p_interp(i,:) = .5*(p_lower_order(i,:) + ...
                            p_lower_order(i+1,:));
    end

    y0 = [z(1:q_dofs,:);zeros(1,dim);p_interp;q1;p1];
    y0 = reshape(y0,dim*N_err,1);
    z = fsolve(SPVI,y0,fsolve_opts);
    %z = inexact_newton(SVI,y0,new_opts) %new_opts must be set 
    z = reshape(z,N_err,dim);

    x_est = [z(end-1,:),z(end,:)]';
    err = norm(x_est-x_next(:,end),2);
  end
end

At lines 22 and 56 I put a comment line to show that it is possible to solve the implicit system also with inexact_newton (which is the one explained in the previous post on backward Euler), so it's not mandatory to use only fsolve. I will make it an option in order to let the user choose which solver to use.

Another important aspect to point out is that in case of an adaptive timestep the error is estimated using a new solution on the same timestep but with polynomials one order higher and one more internal node for the quadrature rule. Furthermore, for this last solution, the starting guess for fsolve is chosen in a non-trivial way: at line 53 we see that y0 has in the first q_dofs_err-1 rows the same modal values calculated before for the new solution x_next and a row of zeros just below. Then the starting nodal values for the momenta are set (lines 47-51) as a trivial average of new solution nodal values. This can seem wrong but I empirically stated that the number of iterations done by fsolve are the same of cases in which the reinitialization of nodal values is more sophisticated, so that is less computationally expensive to do a trivial average.

In the following code is shown the implementation of the system to be solved:

function [F,Jac] = svi_system(y,q_dofs,p_dofs,w,L,D,C, ...
          dim,N,f,t,q0,p0,dt)
  F = zeros(N*dim,1);
  V = zeros(p_dofs,dim*2);
  X = zeros(dim*2,1);
  W = reshape(y,N,dim);
  for i = 1:1:p_dofs
    X = [L(i,:)*W(1:q_dofs,:),W(i+q_dofs,:)]';
    V(i,:) = f(t,X);
  end
  for i = 1:1:dim
    F((1+N*(i-1)):(q_dofs+N*(i-1)),1) = (ones(1,p_dofs)* ...
     (((w.*y((q_dofs+1+N*(i-1)):(q_dofs+p_dofs+N*(i-1))))* ...
     ones(1,q_dofs)).*D + (((dt/2).*w.*V(:,i+dim))* ...
     ones(1,q_dofs)).*L) + (p0(i)*ones(1,q_dofs)).*C(1,:) - ...
      (y(N*i)*ones(1,q_dofs)).*C(2,:))';    
    F(1+N*(i-1)+q_dofs:N*(i-1)+q_dofs+p_dofs,1) = V(:,i) - ...
            (2.0/dt)*(D*y((1+N*(i-1)):(q_dofs+N*(i-1))));
    F(N*i-1) = C(2,:)*y((1+N*(i-1)):(q_dofs+N*(i-1)))-y(N*i-1);
    F(N*i) = C(1,:)*y((1+N*(i-1)):(q_dofs+N*(i-1))) - q0(i);
  end
  if(nargout>1)
    warning('off','Octave:broadcast');
    flag = 0;
    try
      [~,H]=f(0,zeros(2*dim,1));
    catch
      flag = 1;
      warning();
    end
    DV = zeros((dim*2)^2,p_dofs);
    if( flag==1 )
      for i = 1:1:p_dofs
        X = [L(i,:)*W(1:q_dofs,:),W(i+q_dofs,:)]';
        DV(:,i) = differential(f,t,X);
      end
    else
      for i = 1:1:p_dofs
        X = [L(i,:)*W(1:q_dofs,:),W(i+q_dofs,:)]';
        [~,DV(:,i)] = f(t,X);
      end
    end
    DV = DV';
    Jac=zeros(N*dim,N*dim);
    for u=1:1:dim
      [...]
    end
  end
end 

Here, at line 46, I don't show the implementation of the Jacobian of the implicit system because, with this visualization style, it may appear very chaotic. The important aspect is that the user can use the implemented Jacobian to speedup the computation. I stated that it really allows to speedup the computation when using fsolve as solver. Furthermore the code tries to see if the user has implemented, in the function defining the Hamilton's equations, also the Hessian of the Hamiltonian (which is needed int the computation of the Jacobian of the system). If the function passes as second parameter the Hessian, then the program uses that one and the computation is faster, otherwise it approximates the Hessian with the following function and the computation will be slower:

function Hamilt = differential(f,t,x)
  f_x = f(t,x);
  dim_f = length(f_x);
  dim_x = length(x);
  if( dim_f ~= dim_x )
    error('not implemented yet');
  end
  Hamilt = zeros(dim_f*dim_x,1);
  delta = sqrt(eps);
  for i = 1:1:dim_f
    for j = i:1:dim_x
      Hamilt(i+(j-1)*dim_f) = (1/delta)*(f(t,x+delta* ...
         [zeros(j-1,1);1;zeros(dim_x-j,1)])(i) - f_x(i));
      Hamilt(j+(i-1)*dim_x) = Hamilt(i+(j-1)*dim_f);
    end
  end
end

The Hessian of the Hamiltonian must be stored in a particular way (that I have to optimize yet, but the actual one works fine too) which is showed in the following example which is the definition of the Hamilton's equations for the armonic oscillator:

function [dy,ddy] = armonic_oscillator(t,y)
  dy = zeros(size(y));
  d = length(y)/2;
  q = y(1:d);
  p= y(d+1:end);

  dy(1:d) = p;
  dy(d+1:end) = -q;

  H = eye(2*d);
  ddy = transf(H);
end

function v = transf(H)
  [r,c] = size(H);
  v = zeros(r*c,1);
  for i=1:1:r
    for j=1:1:c
      v(i+(j-1)*r) = H(i,j);
    end
  end
end

Symplectic Euler (semi-implicit Euler), Velocity-Verlet and Stormer-Verlet methods

Symplectic Euler method

The symplectic Euler method is a modification of the Euler method and is useful to solve Hamilton's equation of motion, that is a system of ODE where the unknowns are the generalized coordinates q and the generalized momenta p. It is of first order but is better than the classical Euler method because it is a symplectic integrator, so that it yelds better results.

Given a Hamiltonian system with Hamiltonian H=H(t;q,p) then the system of ODE to solve writes: $$\left\{\begin{array}{l} \dot{q} = \frac{dH}{dp}(t;q,p) = f(t;q,p)\\ \dot{p}=-\frac{dH}{dq}(t;q,p) = g(t;q,p) \end{array}\right. $$

From E. Hairer, C. Lubich and G. Wanner, "Geometric Numerical Integration" we can state that the semi-implicit Euler scheme for previous ODEs writes: $$ \left\{\begin{array}{l} p_{k+1} = p_k + h g(t_k;q_k,p_{k+1}) \\ q_{k+1}=q_k + h f(t_k;q_k,p_{k+1}) \end{array}\right. $$ If g(t;q,p) does not depend on p, then this scheme will be totally explicit, otherwise the first equation will be implicit and the second will be explicit.

The function that I created to make use of this integrator is odeSE.m and it passes to the integrate functions the stepper defined in the function symplectic_euler.m. The signature of odeSE is similar to those of the other ode solvers:

[t,y]=odeSE(@ode_fcn,tspan,y0)
[t,y]=odeSE(@ode_fcn,tspan,y0,options)

It's important to know that y0 must be a vector containing initial values for coordinates in its first half and initial values for momenta in its second half; the function handle (or inline function) ode_fcn must take exactly two input arguments (the first is the time and the second is the vector containing coordinates and momenta such as in y0) and returns a vector containing dq/dt in its first half and dp/dt in the second half.

options variable can be set with odeset and, if the system of ODE is explicit, the field options.Explicit can be set to 'yes' in order to speedup the computation. If tspan has only one element, then options.TimeStepNumber and options.TimeStepSize must not be empty, so that it will be used the integrate function integrate_n_steps.

This is the code of the stepper symplectic_euler.m:

function [x_next,err] = symplectic_euler(f,t,x,dt,options)
  dim = length(x)/2;
  q = x(1:dim);
  p = x(dim+1:end);
  if( strcmp(options.Explicit,'yes') )
    p_next = p + dt*(f(t,[q; p])(dim+1:end));
    q_next = q + dt*f(t,[q; p_next])(1:dim);
    x_next = [q_next; p_next];
    if(nargout == 2)
      dt_new = dt/2;
      p_next = p + dt_new*(f(t,[q; p])(dim+1:end));
      q_next = q + dt_new*f(t,[q; p_next])(1:dim);
      q = q_next;
      p = p_next;
      t = t+dt_new;
      p_next = p + dt_new*(f(t,[q; p])(dim+1:end));
      q_next = q + dt_new*f(t,[q; p_next])(1:dim);
      x_est = [q_next;p_next];
      err = norm(x_next-x_est,2);
    end
  else
    SE1 = @(y) (y-p-dt*(f(t,[q; y])(dim+1:end)));
    p_next = fsolve(SE1,zeros(size(p)));
    q_next = q + dt*f(t,[q; p_next])(1:dim);
    x_next = [q_next; p_next];
    if(nargout == 2)
      dt_new = dt/2;
      SE1 = @(y) (y-p-dt_new*(f(t,[q; y])(dim+1:end)));
      p_next = fsolve(SE1,zeros(size(p)));
      q_next = q + dt_new*f(t,[q; p_next])(1:dim);
      q = q_next;
      p = p_next;
      t = t+dt_new;
      SE1 = @(y) (y-p-dt_new*(f(t,[q; y])(dim+1:end)));
      p_next = fsolve(SE1,zeros(size(p)));
      q_next = q + dt_new*f(t,[q; p_next])(1:dim);
      x_est = [q_next;p_next];
      err = norm(x_next-x_est,2);
    end
  end
end

Velocity-Verlet method

The velocity-Verlet method is a numerical method used to integrate Newton's equations of motion. The Verlet integrator offers greater stability, as well as other properties that are important in physical systems such as preservation of the symplectic form on phase space, at no significant additional cost over the simple Euler method.

If we call x the coordinate, v the velocity and a the acceleration then the equations of motion write: $$\left\{ \begin{array}{l} \frac{dx}{dt} = v(t,x)\\ \frac{dv}{dt} = a(t,x) \end{array} \right. $$

So that, given the initial conditions (coordinates and velocities), the velocity-verlet scheme writes: $$ \left\{ \begin{array}{l} x_{k+1} = x_k + h v_k + 0.5 h^2 a_k\\ v_{k+1} = v_k + 0.5 h (a_k + a_{k+1}) \end{array}\right. $$ where $$ a_{k+1} = a(t_{k+1},x_{k+1})\ .$$

This method is one order better than the symplectic Euler method. The global error of this method is of order two. Furthermore, if the acceleration indeed results from the forces in a conservative mechanical or Hamiltonian system, the energy of the approximation essentially oscillates around the constant energy of the exactly solved system, with a global error bound again of order two.

The function that uses velocity-Verlet scheme is odeVV.m and the stepper called at each iteration is velocity_verlet.m. The signature of odeVV is the same of those of others ODE solvers:

[t,y]=odeVV(@ode_fcn,tspan,y0)
[t,y]=odeVV(@ode_fcn,tspan,y0,options)

The documentation of input arguments is the same descripted in the previous symplectic Euler section, but there is the difference that now the function ode_fcn must return a vector containing the velocities in its first half and the expression of the acceleration in the second half.

This is the code of the stepper symplectic_euler.m:

function [x_next,err] = velocity_verlet(f,t,x,dt,options)
  dim = length(x)/2;
  q = x(1:dim);
  v = x(dim+1:end);
  a = f(t,x);
  q_next = q + dt*v + .5*dt*dt*a(dim+1:end);
  v_next = v + .5*dt*((a+f(t+dt,[q_next; v]))(dim+1:end));
  x_next = [q_next; v_next];
  if(nargout == 2)
    dt_new = dt/2;
    q_next = q + dt_new*v + .5*dt_new*dt_new*a(dim+1:end);
    v_next = v + .5*dt_new*((a + ...
                      f(t+dt_new,[q_next;v]))(dim+1:end));
    t = t+dt_new;
    q = q_next;
    v = v_next;
    a = f(t,[q; v]);
    q_next = q + dt_new*v + .5*dt_new*dt_new*a(dim+1:end);
    v_next = v + .5*dt_new*((a + ...
                      f(t+dt_new,[q_next;v]))(dim+1:end));
    x_est = [q_next; v_next];
    err = norm(x_next - x_est,2);
  end
end

Stormer-Verlet method

The Stormer-Verlet scheme is a symplectic integrator of order two, useful to integrate Hamiltonian systems in the form described in the previous symplectic Euler section. It is characterized by the approximation of the integral defining the discrete Lagrangian $$ L_d(q_k,q_{k+1})\approx\int_{t_k}^{t_{k+1}}L(t;q,\dot{q})dt $$ with the trapezoidal rule. So that $$ L_d(q_k,q_{k+1}) = \frac{h}{2}\bigg( L\Big(q_k,\frac{q_{k+1}-q_k}{h}\Big) + L\Big(q_{k+1},\frac{q_{k+1}-q_k}{h}\Big) \bigg)\ .$$

Considering again the system: $$ \left\{\begin{array}{l} \dot{q} = \frac{dH}{dp}(t;q,p) = f(t;q,p)\\ \dot{p}=-\frac{dH}{dq}(t;q,p) = g(t;q,p) \end{array}\right. \ ,$$ the scheme implemented for this integrator (described in E. Hairer, C. Lubich and G. Wanner, "Geometric Numerical Integration") writes as follows: $$ \left\{\begin{array}{l} p_{k+\frac{1}{2}} = p_k + \frac{h}{2}g(t_k;q_k,p_{k+\frac{1}{2}}) \\ q_{k+1}=q_k+\frac{h}{2}\Big( f(t_k;q_k,p_{k+\frac{1}{2}}) + f(t_{k+1};q_{k+1},p_{k+\frac{1}{2}})\Big) \\ p_{k+1} = p_{k+\frac{1}{2}} + \frac{h}{2} g(t_{k+\frac{1}{2}};q_{k+1},p_{k+\frac{1}{2}}) \end{array}\right. $$

The function in which this method is implemented is odeSV.m and it calls the stepper stormer_verlet.m. The documentation is the same of odeSE.m. Its implementation is the following:

function [x_next,err] = stormer_verlet(f,t,x,dt,options)
  dim = length(x)/2;
  q = x(1:dim);
  p = x(dim+1:end);
  if( strcmp(options.Explicit,'yes') )
    p_mid = p + .5*dt*(f(t,[q; p])(dim+1:end));
    q_next = q + .5*dt*((f(t,[q; p_mid])(1:dim))+ ...
                         (f(t+dt,[q;p_mid])(1:dim)));
    p_next = p_mid +.5*dt*(f(t+dt,[q_next;p_mid])(dim+1:end));
    x_next = [q_next; p_next];
    if(nargout == 2)
      dt_new = dt/2;
      p_mid = p + .5*dt_new*(f(t,[q; p])(dim+1:end));
      q_next = q + .5*dt_new*((f(t,[q; p_mid])(1:dim))+ ...
                           (f(t+dt_new,[q;p_mid])(1:dim)));
      p_next = p_mid + .5*dt_new* ...
                   (f(t+dt_new,[q_next;p_mid])(dim+1:end));
      q = q_next;
      p = p_next;
      t = t+dt_new;
      p_mid = p + .5*dt_new*(f(t,[q; p])(dim+1:end));
      q_next = q + .5*dt_new*((f(t,[q; p_mid])(1:dim))+ ...
                           (f(t+dt_new,[q;p_mid])(1:dim)));
      p_next = p_mid + .5*dt_new* ...
                   (f(t+dt_new,[q_next;p_mid])(dim+1:end));
      x_est = [q_next; p_next];
      err = norm(x_est - x_next,2);
    end
  else
    SV1 = @(y) (y - p - .5*dt*(f(t,[q;y])(dim+1:end)));
    p_mid = fsolve(SV1,zeros(size(p)));
    SV2 = @(y) (y - q - .5*dt*((f(t,[q;p_mid])(1:dim))+ ...
                              (f(t+dt,[y;p_mid])(1:dim))));
    q_next = fsolve(SV2,q);
    p_next = p_mid + .5*dt* ...
                         f(t+dt,[q_next;p_mid])(dim+1:end);
    x_next = [q_next;p_next];
    if(nargout == 2)
      dt_new = dt/2;
      SV1 = @(y) (y - p - .5*dt_new*(f(t,[q;y])(dim+1:end)));
      p_mid = fsolve(SV1,zeros(size(p)));
      SV2 = @(y) (y - q - .5*dt_new* ...
   ((f(t,[q;p_mid])(1:dim))+(f(t+dt_new,[y;p_mid])(1:dim))));
      q_next = fsolve(SV2,q);
      p_next = p_mid + .5*dt_new* ...
                       f(t+dt_new,[q_next;p_mid])(dim+1:end);
      q = q_next;
      p = p_next;
      t = t+dt_new;
      SV1 = @(y) (y - p - .5*dt_new*(f(t,[q;y])(dim+1:end)));
      p_mid = fsolve(SV1,zeros(size(p)));
      SV2 = @(y) (y - q - .5*dt_new* ...
   ((f(t,[q;p_mid])(1:dim))+(f(t+dt_new,[y;p_mid])(1:dim))));
      q_next = fsolve(SV2,q);
      p_next = p_mid + .5*dt_new* ...
                       f(t+dt_new,[q_next;p_mid])(dim+1:end);
      x_est = [q_next; p_next];
      err = norm(x_next-x_est,2);  
    end
  end
end

Backward Euler with Inexact Newton solver

Backward Euler

The backward Euler method is one of the most common methods used to solve ODEs. It is an implicit method of order one but its strenght lies in its A-stability property.

Given a system of ODEs of the form: $$\frac{d\mathbf{y}}{dt} = \mathbf{f}(t,\mathbf{y}) \\ \mathbf{y}(t_0) = \mathbf{y}_0\ ,$$ the backward Euler method is written as: $$ \mathbf{y}_{k+1} = \mathbf{y}_k + \Delta t\ \mathbf{f}(t_{k+1},\mathbf{y}_{k+1})\ .$$

If we define u to be the exact solution then we have that: $$ \frac{d\mathbf{u}}{dt} = \mathbf{f}(t,\mathbf{u}) \\ \mathbf{u}(t_k) = \mathbf{y}_k $$ So that $$\mathbf u(t_{k+1}) = \mathbf u(t_k) + \int_{t_k}^{t_{k+1}} \mathbf{f}(t,\mathbf u) dt\ ,$$ which, under some conditions of regularity, can be rewritten as $$\mathbf u(t_{k+1}) = \Delta t\ \mathbf f(t_{k+1},\mathbf u(t_{k+1})) - \frac{{\Delta t}^2}{2} \mathbf f'(\xi)\ .$$

From these results it can be derived that the local truncation error (LTE) is $$ \mathbf e_{k+1} = \mathbf y_{k+1} - \mathbf u(t_{k+1}) = \mathbf o({\Delta t}^2)\ .$$

Inexact Newton

Now we have to solve the non-linear system built up with the backward Euler method. We want to use an inexact Newton solver. The classical Newton method is written as: given a non-linear functions system F(x) = 0 with a given initial condition, solve the linearized system: $$ \mathbf F'(\mathbf x_k)\mathbf s_k + \mathbf F(\mathbf x_k) = \mathbf 0 \\ \mathbf x_{k+1} = \mathbf x_k + \mathbf s_k\ .$$ Go on iterating this method until a given tolerance is reached.

In many cases, especially when we don't have an explicit expression for the Jacobian or it can't be inverted, we can use an inexact Newton method: $$ || \mathbf F'(\mathbf x_k)\mathbf s_k + \mathbf F(\mathbf x_k) || \leq \eta_k || \mathbf F(\mathbf x_k)||\ ,$$ where \(\eta_k\) is said to be the forcing term. The linearized system can be solved with an iterative method such as GMRES or preconditioned conjugate gradient.

There are two possible choices for the forcing term:

1. first choice $$ \eta_k = \frac{|| \mathbf F(\mathbf x_k)-\mathbf F(\mathbf x_{k-1}) -\mathbf F'(\mathbf x_{k-1})\mathbf s_{k-1} ||}{||\mathbf F(\mathbf x_{k-1}) ||}\ ;$$
2. second choice $$ \eta_k = \gamma \Big(\frac{|| \mathbf F(\mathbf x_k) ||}{|| \mathbf F(\mathbf x_{k-1}) ||}\Big)^{\alpha}\ .$$

Choosing the forcing term with one of the two previous possibilities, rather than imposing a fixed tolerance, it is possible to avoid the phenomenon of oversolving. Another advantage is that if we use the GMRES as the linear solver it's not necessary to know the Jacobian nor to datermine it, because we can approximate the term $$ \mathbf F'(\mathbf x_{k})\mathbf s_k $$ with a first order finite differences formulae: $$ \mathbf F'(\mathbf x_k)\mathbf s_k \approx \frac{\mathbf F(\mathbf x_k + \delta\mathbf s_k) - \mathbf F(\mathbf x_k)}{\delta} \ .$$

The signature of inexact_newton is the same of fsolve:

[x,fval,exitflag,output,jacobian]=inexact_newton(fun,x0,options)

Ode solvers options

Ode solvers options

What I want to do at this moment of the project is to add to all the solvers I've written since now, the ability to manage all the possible options of the ODE structure (defined previously in odeset structure) that are already managed in corresponding MATLAB ode solvers; so that the solvers I've written will be totally MATLAB-compatible and this will be the same also for the new solvers I'm going to implement.

Timeline

Timeline

The following list explains the main points of the timeline for my SOCIS project:

1. adapt odeset and odeget from Octave odepkg to be completely MATLAB-compatible: all option names supported by MATLAB should be available; furthermore we want the possibility to add new option names that will be useful for example for symplectic integrators;
2. implement a new code structure for one-step integrators ode12, ode23, ode45, odefwe (forward Euler), odebwe (backward Euler) and inexact solvers, subdividing the execution into 'blocks' in order to be able to do an easier optimization and an easier extension to new functionalities; in particular implement the steppers which will be the functions called to do one integration step and then implement the integrate functions which will integrate, in different ways, over the whole time interval, calling the respective steppers. Afterwards adapt the interface of ode45 to be completely MATLAB-compatible;
3. implement MATLAB-compatible version of deval;
4. implement the following geometric integrators as new steppers:
   • symplectic Euler;
   • Stormer-Verlet;
   • velocity-Verlet;
   • spectral variational integrators;
   • SHAKE;
   • RATTLE.
5. optimization, error checking and documentation of all the code written.

At the moment point 1 is completed; in point 2 I still have to adapt the interface of ode45 and to refine the implemented code. Point 3 has to be done. In point 4 I have implemented sympelctic Euler, Stormer-Verlet and spectral variational integrators but these need to be adapted to the code structure 'steppers', 'integrate functions'; the other integrators must be implemented. Finally, during the implementation I will start to write a bit of documentation and at the end I will do optimizations and debugging.

Odeset & odeget

Odeset & odeget

Odeset and odeget functions allow to build a valid ODE options structure. They are already available in Octave odepkg but they are not perfectly compatible with MATLAB odeset and odeget functions. Furthermore, for geometric integrators like Spectral Variational Integrators,
I will need new options into the ODE options structure, which now are not admitted in Octave odepkg.

So that I have written my own versions of odeget.m, odeset.m and the new function ode_struct_value_check.m in order to have full compatibility with MATLAB odeset and odeget (their same behaviour on the basis of their official documentation) and also the possibility to add new field names which i will need for future solvers of this project.

The new fields are the union of MATLAB ODE options and Octave ODE options, plus my new options (default values are in square brackets):

'             AbsTol:  scalar or vector, >0, [1.e-6]'
'                BDF:  binary, {on, [off]}'
'             Events:  '
'           Explicit:  binary, {yes, no, []}'
'      InexactSolver:  switch, {inexact_newton, []}'
'       InitialSlope:  vector, []'
'        InitialStep:  scalar, >0, []'
'           Jacobian:  matrix or function_handle, []'
'          JConstant:  binary, {on, [off]}'
'           JPattern:  sparse matrix, []'
'               Mass:  matrix or function_handle, []'
'       MassConstant:  binary, {on, [off]}'
'       MassSingular:  switch, {yes, [maybe], no}'
'MaxNewtonIterations:  scalar, integer, >0, [1.e3]'
'           MaxOrder:  switch, {1, 2, 3, 4, [5]}'
'            MaxStep:  scalar, >0, []'
'   MStateDependence:  switch, {none, [weak], strong}'
'          MvPattern:  sparse matrix, []'
'          NewtonTol:  scalar or vector, >0, []'
'        NonNegative:  vector of integers, []'
'        NormControl:  binary, {on, [off]}'
'          OutputFcn:  function_handle, []'
'         OutputSave:  scalar, integer, >0, []'
'          OutputSel:  scalar or vector, []'
'   PolynomialDegree:  scalar or vector, integer, >0, []'
'    QuadratureOrder:  scalar or vector, integer, >0, []'
'             Refine:  scalar, integer, >0, []'
'             RelTol:  scalar, >0, [1.e-3]'
'              Stats:  binary, {on, [off]}'
'     TimeStepNumber:  scalar, integer, >0, []'
'       TimeStepSize:  scalar, >0, []'
'         Vectorized:  binary, {on, [off]}'

The usage of odeset will be one of the following:

odeset
opt = odeset()
opt = odeset(old_ODE_options,new_ODE_options)
opt = odeset(old_ODE_options,'field1',value1,'field2',value2,...)
opt = odeset('field1',value1,'field2',value2,'field3',value3,...)

The usage of odeget will be one of the following:

option_value = odeget(ODE_structure,'field')
option_value = odeget(ODE_structure,'field',default_value)
option_value = odeget(any_struct,'field',default_value,'fast')

The last usage is needed for MATLAB compatibility and represents a fast access to the given structure with no error checking on options names.

Fuzzy string comparison

Users may do mistakes while typing ODE option names or they may want to write everything in low cases for a faster coding. So that my odeset and odeget functions behave in a user-friendly way, calling the fuzzy_compare function which determines the distance of each structure option name from the given word and returns the indices of those options whose distance is under a given tolerance.

The idea is that word cases are not relevant, if we want an exact matching then the tolerance will be 0, if no index is returned then it will definitely be a typing error, if one index is returned but words don't match exactly a warning will be displayed saying that the program is going on assuming that the user was intending the closest option name, otherwise all the fields whose distance is under the tolerance will be displayed and the user will be asked to insert the name again.

Signature of this function follows.

res = fuzzy_compare(string1,string_set,correctness)

• string1 must be a string and is the option name we're looking for;
• string_set must be a cell array of strings or a column vector of strings; it represents the set of option names;
• correctness is the tolerance we want. It is an optional input argument;

We may set correctness equal to 'exact' or 0 and in this case fuzzy_compare will look only for exact matching; if we set correctness to any positive integer it will be the tolerance of the method (that is the maximum distance, from the given string, accepted). If we don't specify anything for correctness then the tolerance will be calculated as a function of the minimum distance and of the length of the given string: the less the distance, the less the tolerance; the greater the length, the greater the tolerance.

tolerance = minimus + floor(length(string1)/4)*floor(minimus/3)

There exist many definitions of distance between strings. I've chosen the Levenshtein distance.

Levenshtein distance

The Levenshtein distance is a string metric and is equal to the minimum number of single-characters edit (insertion, deletion and substitution) required to change one word into the other. This is the main algorithm of levenshtein function:

function [dist,d] = levenshtein(string1,string2)
  m = length(string1);
  n = length(string2);

  d = zeros(m+1,n+1);
  d(2:m+1,1) = [1:1:m];
  d(1,2:n+1) = [1:1:n];

  for j=2:1:n+1
    for k=2:1:m+1
      if(string1(k-1)==string2(j-1))
        d(k,j)=d(k-1,j-1);
      else
        d(k,j)=min(d(k-1,j)+1,min(d(k,j-1)+1,d(k-1,j-1)+1));
      end
    end
  end

  dist = d(m+1,n+1);
end

All the code is available under the terms of GNU General Public License at my_octave-odepkg.

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)

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)