### 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