Простое гармоническое движение - Verlet - Внешняя сила - Matlab

Question

Я пробежался по алгебре, которую я ранее делал для метода Верлета, без силы - это привело к тому же коду, как вы видите ниже, но термин "+(2*F/D)" отсутствовал, когда я игнорировал внешнюю силу, Алгоритм работал точно, как и ожидалось, однако для следующих параметров:

м = 7; к = 8; б = 0,1; params = [m, k, b];

(и размер шага h = 0,001)

сила, намного превышающая что-то вроде 0,00001, слишком велика. Я подозреваю, что пропустил трюк с алгеброй.

Мой вопрос заключается в том, может ли кто-то заметить изъян в моем добавлении силового термина в моем методе Верле

% 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

Answer 1

Как написано в ответе на предыдущий вопрос, в момент, когда трение входит в уравнение, система перестает быть консервативной и имя "Verlet" больше не применяется. Это все еще действительная дискретизация

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

(с небольшой ошибкой с большими последствиями).

Дискретизация использует центральные разности первого и второго порядка

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)

в результате чего

(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)

---

Ошибка: как видите, у вас отсутствует фактор h^2 в срок с F,