algorithm - MATLAB:Verlet 算法 -

Question

下面是我的 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)

% vin is the particle vector (xn,yn)
x0 = vinverletx(1);
x1 = vinverletx(2);


% find the verlet coefficients (F=0)
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);

vout = x2;

% vout is the particle vector (xn+1,yn+1)
end 

我写了一个脚本来测试这个功能。上下文是简谐运动，将测试 Verlet 算法与其他算法相比的相对精度。

这是我的测试脚本:

% verlet test
clear all
close all

% don't define fixed paramaters every loop
h = 0.001;
m = 7.4; % Mass
k = 7.7; % Return force
b = 0; % Drag
params = [m,k,b];


% verlet
x=2; % define initial values and put into vector form
v=0;
vin = [x,v];
vstorex = vin(1);
vstorev = vin(2);

for n=1:200
   if n == 1 
      vnext = eulerintegrate(vin,n,h,params); % the next position and velocity
      vstorex = [vstorex;vnext(1)]; %#ok<*AGROW> % store xn and vn+1
      vinverletx = [vin(1),vnext(1)]; % now we have two x values for the verlet algorithm!
   else if n ==2
      xnext=verlet(vinverletx,h,params);
      vstorex = [vstorex;xnext];
       else
            vinverletx = [vstorex(n),vstorex(n-1)];
            xnext=verlet(vinverletx,h,params); 
            vstorex = [vstorex;xnext];
       end
   end
end

plot(vstorex);

生成的图在步长为 0.001 的 200 次运行中大幅膨胀 - http://i.imgur.com/GF2Zdvu.png

这是步长为 0.0001 的 200 次运行:http://i.imgur.com/u0zCUWS.png

正如您很容易看出的那样，它也会以类似的方式爆炸。我的代码中一定有问题(我看不到)。

提前致谢!

Answer 1

你的微分方程是 x''=a(x)=-k/m*x，带有基本 Verlet 方法的中点公式

x0-2*x1+x2= h*h*a(x1)

你得到了

x2 = -x0+(2-h*h*k/m)*x1

要获得正确的错误顺序，您需要尽可能最好的初始化，即

x1 = x0 + v0*h + 0.5*a(x0)*h*h

---

您不能在存在拖动的情况下使用 Verlet 方法。或者至少您不能指望它具有广告宣传的特性。这些仅适用于保守系统，其中力来自势场，并且仅来自势场。

---

错误

在该过程中，您希望这两个值按递增的索引顺序排列。在函数调用中，您按索引降序构造输入。除了更正该错误之外，我还会更改整个循环以将其简化为

vin = [x,v];
vnext = eulerintegrate(vin,n,h,params); % the next position and velocity
vstorex = [vin(1), vnext(1)]; % or to the same effect: [x, x+h*v]


for n=2:200
   vinverletx = [vstorex(n-1),vstorex(n)];
   xnext=verlet(vinverletx,h,params); 
   vstorex = [vstorex;xnext];
end