簡單的諧波運動 - Verlet - 外力 - Matlab

我跑過代碼，我以前在沒有力的情況下爲Verlet方法完成 - 這導致了與您在下面看到的相同的代碼，但「+（2 * F/D）「當我忽視外力時，失去了這個術語。如預期的那樣，該算法準確地工作，但是對於以下參數：簡單的諧波運動 - Verlet - 外力 - Matlab

m = 7; k = 8; b = 0.1; params = [m，k，b];

（和步長H = 0.001）

遠高於類似0.00001的力是太大了。我懷疑我已經錯過了代數的一個技巧。

我的問題是，是否有人能發現這個漏洞在我另外一個力項在我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,F) 

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


% find the verlet coefficients 
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)+(2*F/D); 

vout = x2; 

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

來源

2015-04-23 Idios

您的問題尚不清楚，請更詳細地說明它 - 請參閱http://stackoverflow.com/help/mcve –

您可以提供指向上一個問題的鏈接https://stackoverflow.com/q/29701350/3088138哪裏有更多的代碼。 – LutzL

如果將參數擴展爲'm = param（1）'，'k = param（2）'，'b = param（3）'，它會使代碼更具可讀性。 – LutzL

寫在回答前一個問題，那一刻摩擦進入方程，該系統不再保守，名稱「Verlet」不再適用。它仍然是一個有效的離散化的

m*x''+b*x'+k*x = F

（有一些輕微的錯誤後果很大）。

離散化採用一階和二階

x'[k] = (x[k+1]-x[k-1])/(2*h) + O(h^2) 
x''[k] = (x[k+1]-2*x[k]+x[k-1])/(h^2) + O(h^2)

導致

(2*m+b*h)*x[k+1] - 2*(2*m+h^2*k) * x[k] + (2*m-b*h)*x[k-1] = 2*h^2 * F[k] + O(h^4)

錯誤的中心差商：正如你可以看到，你缺少的因素h^2在與F的期限。

來源

2015-04-24 14:06:47 LutzL

簡單的諧波運動 - Verlet - 外力 - Matlab

回答

相關問題