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