简单的谐波运动 - Verlet - 外力 - Matlab

Question

我浏览了我之前在没有力的情况下为 Verlet 方法所做的代数 - 这导致了与您在下面看到的相同的代码，但是当我忽略外力时，"+（2*F/D）"项丢失了。该算法如预期的那样准确工作，但是对于以下参数：

m = 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的因素。