1

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

正如你很容易看出的那样，它同样会爆炸。我的代码中一定有问题（我看不到）。

提前致谢！

4

1 回答 1

3

你的微分方程是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
于 2015-04-17T15:22:09.817 回答