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