0

我在没有力的情况下运行了我之前为 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
4

1 回答 1

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.

于 2015-04-24T14:06:47.790 回答