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