MATLAB：verlet的算法 - [英] MATLAB: Verlet Algorithm -

问题描述

下面是我的code为verlet的功能，被称为从我的主脚本。

Below is my code for the Verlet function, to be called from my main script.

% 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的算法将被测试的相对精度相比其他算法

I made a script to test this function. The context is simple harmonic motion, and the Verlet algorithm is going to be tested for relative accuracy compared to other algorithms.

下面是我的测试脚本：

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

产生的阴谋炸毁了巨大的200奔跑0.001步长 - HTTP：//i.imgur。 COM / GF2Zdvu.png

下面是200奔跑0.0001步长： http://i.imgur.com/u0zCUWS.png

Here is 200 runs of 0.0001 step size: http://i.imgur.com/u0zCUWS.png

这打击了类似的，你可以很容易地分辨。必须有一个问题（我看不到），在我的code。

It blows up similarly, as you can easily tell. There must be a problem (which I can't see) in my code.

在此先感谢！

推荐答案

您微分方程是 X''= A（X）= - K / M * X ，与基本verlet的方法的中点式

Your differential equation is x''=a(x)=-k/m*x, with the midpoint formula of the basic Verlet method

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

你

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

要得到正确的错误命令，你需要最好的初始化可能的，这是

To get the correct error order, you need the best initialization possible, which is

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

---

您不能使用拖动的presence的verlet的方法。或者至少，你不能指望它有做广告的属性。那些只为保守系统，其中力的结果从一个潜在的领域，仅从这一点。

---

You can not use the Verlet method in the presence of drag. Or at least you can not expect it to have the advertized properties. Those only hold for conservative systems, where the force results from a potential field, and only from that.

在过程中，您期望在增加索引顺序的两个值。在函数调用，构造输入递减索引顺序。除了纠正这个错误，我会改变整个循环将其简化为

In the procedure you expect the two values in increasing index order. In the function call, you construct the input in decreasing index order. Apart from correcting that error, I would change the whole loop to simplify it to

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

这篇关于MATLAB：verlet的算法 - 的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！