MATLAB: Verlet Algorithm -

Question

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.stack.imgur.com/1rAVc.png

Here is 200 runs of 0.0001 step size: https://i.stack.imgur.com/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!

I'd expect the output to be a sinusoidal graph, with a relatively constant amplitude - about 2 (metres), similar to the initial amplitude. — james.sw.clark, Apr 17 '15 at 14:16
My analytical solution looks like this, over 5 seconds: http://i.imgur.com/8gwC2Zd.png — james.sw.clark, Apr 17 '15 at 14:19

Lutz Lehmann · Accepted Answer · 2015-04-17T16:41:00.100

3

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

edited Apr 17 '15 at 16:41

answered Apr 17 '15 at 15:22

Lutz Lehmann

25,219
2
22
51

Thank you very much. My Verlet algorithm is now working as expected! May I ask how you spotted my error? – james.sw.clark Apr 17 '15 at 18:31
First comparing that the actual computations were correct, as they are, then looking for a line that did not conform to expectations, esp. the argument sequences in the function calls. `n` followed by `n-1` was a pattern that did not fit. – Lutz Lehmann Apr 17 '15 at 18:35
Thank you, again. Do you have any advice on how to have a clarity of mind and incisive mindset when debugging? – james.sw.clark Apr 17 '15 at 18:37

MATLAB: Verlet Algorithm -

1 Answers1

Error

Linked