0

I ran through the algebra which I had previously done for the Verlet method without the force - this lead to the same code as you see below, but the "+(2*F/D)" term was missing when I ignored the external force. The algorithm worked accurately, as expected, however for the following parameters:

m = 7 ; k = 8 ; b = 0.1 ; params = [m,k,b];

(and step size h = 0.001)

a force far above something like 0.00001 is much too big. I suspect I've missed a trick with the algebra.

My question is whether someone can spot the flaw in my addition of a force term in my Verlet method

% 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,F)

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


% find the verlet coefficients
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)+(2*F/D);

vout = x2;

% vout is the particle vector (xn+1,yn+1)
end
james.sw.clark
  • 357
  • 1
  • 8
  • 20
  • Your question is not clear, please rephrase it with more detail - see http://stackoverflow.com/help/mcve – Buck Thorn Apr 23 '15 at 21:57
  • You could have provided a link to the previous question https://stackoverflow.com/q/29701350/3088138 where there is more of the code. – Lutz Lehmann Apr 24 '15 at 13:49
  • It would make the code more readable if you expanded the parameters into `m=param(1)`, `k=param(2)`, `b=param(3)`. – Lutz Lehmann Apr 24 '15 at 13:50

1 Answers1

1

As written in the answer to the previous question, the moment friction enters the equation, the system is no longer conservative and the name "Verlet" does no longer apply. It is still a valid discretization of

m*x''+b*x'+k*x = F

(with some slight error with large consequences).

The discretization employs the central difference quotients of first and second order

x'[k]  = (x[k+1]-x[k-1])/(2*h) + O(h^2)
x''[k] = (x[k+1]-2*x[k]+x[k-1])/(h^2) + O(h^2)

resulting in

(2*m+b*h)*x[k+1] - 2*(2*m+h^2*k) * x[k] + (2*m-b*h)*x[k-1] = 2*h^2 * F[k] + O(h^4)

Error: As you can see, you are missing a factor h^2 in the term with F.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51