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

簡単にわかるように、同様に爆発します。私のコードには (私には見えない) 問題があるはずです。

前もって感謝します！