MATLAB: Verlet Algorithm -

Below is my code for the Verlet function, to be called from my main script.

% verlet.m
% uses the verlet step algorithm to integrate the simple harmonic
% oscillator.

% stepsize h, for a second-order ODE

function vout = verlet(vinverletx,h,params)

% vin is the particle vector (xn,yn)
x0 = vinverletx(1);
x1 = vinverletx(2);


% find the verlet coefficients (F=0)
D = (2*params(1))+(params(3)*h);
A = (2/D)*((2*params(1))-(params(2)*h^2));
B=(1/D)*((params(3)*h)-(2*params(1)));

x2 = (A*x1)+(B*x0);

vout = x2;

% vout is the particle vector (xn+1,yn+1)
end 

I made a script to test this function. The context is simple harmonic motion, and the Verlet algorithm is going to be tested for relative accuracy compared to other algorithms.

Here is my testing script:

% verlet test
clear all
close all

% don't define fixed paramaters every loop
h = 0.001;
m = 7.4; % Mass
k = 7.7; % Return force
b = 0; % Drag
params = [m,k,b];


% verlet
x=2; % define initial values and put into vector form
v=0;
vin = [x,v];
vstorex = vin(1);
vstorev = vin(2);

for n=1:200
   if n == 1 
      vnext = eulerintegrate(vin,n,h,params); % the next position and velocity
      vstorex = [vstorex;vnext(1)]; %#ok<*AGROW> % store xn and vn+1
      vinverletx = [vin(1),vnext(1)]; % now we have two x values for the verlet algorithm!
   else if n ==2
      xnext=verlet(vinverletx,h,params);
      vstorex = [vstorex;xnext];
       else
            vinverletx = [vstorex(n),vstorex(n-1)];
            xnext=verlet(vinverletx,h,params); 
            vstorex = [vstorex;xnext];
       end
   end
end

plot(vstorex);

The plot produced blows up hugely for 200 runs of 0.001 step size - https://i.sstatic.net/1rAVc.png

Here is 200 runs of 0.0001 step size: https://i.sstatic.net/yqJkt.png

It blows up similarly, as you can easily tell. There must be a problem (which I can't see) in my code.

Thanks in advance!

Your differential equation is x''=a(x)=-k/m*x, with the midpoint formula of the basic Verlet method

x0-2*x1+x2= h*h*a(x1)

you get

x2 = -x0+(2-h*h*k/m)*x1

To get the correct error order, you need the best initialization possible, which is

x1 = x0 + v0*h + 0.5*a(x0)*h*h

---

You can not use the Verlet method in the presence of drag. Or at least you can not expect it to have the advertized properties. Those only hold for conservative systems, where the force results from a potential field, and only from that.

---

Error

In the procedure you expect the two values in increasing index order. In the function call, you construct the input in decreasing index order. Apart from correcting that error, I would change the whole loop to simplify it to

vin = [x,v];
vnext = eulerintegrate(vin,n,h,params); % the next position and velocity
vstorex = [vin(1), vnext(1)]; % or to the same effect: [x, x+h*v]


for n=2:200
   vinverletx = [vstorex(n-1),vstorex(n)];
   xnext=verlet(vinverletx,h,params); 
   vstorex = [vstorex;xnext];
end