MATLAB: алгоритм Верле -

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 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

Он взрывается так же, как вы можете легко сказать. В моем коде должна быть проблема (которую я не вижу).

Заранее спасибо!

Answer 1

Ваше дифференциальное уравнение x''=a(x)=-k/m*x, с формулой средней точки основного метода Верле

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