MATLAB Implementation of AVI

I work with AVI in my research so I have included a MATLAB implementation of AVI for the 1-D harmonic oscillator. The code will solve the equation a + gamma * v + (k1 + k2 + k3) x = 0 with any initial conditions x(0) and v(0). Here the spring constant has been artifically split into three spring constants to simulate multiple potentials. If there is only one potential AVI simplifies to the usual Velocity Verlet integrator. The friction term is absorbed into the k1 term in the implementation. The main idea in the implementation is to construct the propagation matrix for the system (x,v) for the different potentials.

To use the code to solve the simple harmonic oscillator a + 3 x = 0 over two periods with the Velocity Verlet algorithm let

ts1 = ts2 = ts3 = 0.01; k1 = k2 = k3 = 1; x0 = 1, v0 = 0, gamma = 0, and T = 4*pi/sqrt(3).

To solve the same system with AVI you can choose

ts1 = 0.005, t2 = 0.01, t3 = 0.02; k1 = 1.5, k2 = 1, k3 = 0.5; and the same values for the other parameters.

For more information regarding AVI please look at Variational Time Integrators by Lew et al.