简谐运动 - Verlet - 外力 - Matlab

Question

我运行通过代数，我之前在没有力的情况下为 Verlet 方法完成了 - 这导致了与您在下面看到的相同的代码，但是 "+(2*F/D )" 当我忽略外力时，术语丢失了。该算法如预期的那样准确地工作，但是对于以下参数：

米=7; k = 8 ; b = 0.1 ; 参数 = [m,k,b];

（步长 h = 0.001）

远高于 0.00001 之类的力太大了。我怀疑我错过了代数技巧。

我的问题是是否有人能发现我在 Verlet 方法中添加力项的缺陷

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

错误： 如您所见，您在 F.

项中缺少一个因子 h^2

简谐运动 - Verlet - 外力 - Matlab

Simple Harmonic Motion - Verlet - External force - Matlab

matlab

verlet-integration