Matlab的：ODE45輸出不正確的用力彈簧質量阻尼器

Question

這裏是我的Matlab代碼來解決二階常微分方程的質量彈簧阻尼系統：

function Spring 

clear all; 
close all; 
options=odeset('RelTol',1e-6); 
p0 = [1 0]; %initial position and velocity 
[t,p] = ode45(@SpringFunction, [0 20], p0, options); 

plot(t,p(:,1)) %plot the first column (x) vs. time 

return 

function pdot = SpringFunction(t,p) 

c = 5; w = 2; 
g = sin(t); % forcing function 
pdot = zeros(size(p)); 
pdot(1) = p(2); 
pdot(2) = g-c*p(2)-(w^2)*p(1); 

return 

我相信我得到的輸出是錯誤的，因爲對於這種情況，我認爲位移與時間的關係應該看起來像是一個幅度減小的正弦曲線。相反，它看起來像這樣：

這對我來說似乎不正確，但請糾正我，如果我錯了。我看不到代碼中的錯誤。

Answer 1

你正在強制阻尼系統，所以你應該期望穩定狀態是一個正弦曲線。這裏的a nice vibrations tutorial（PDF） - 見第448-450頁關於阻尼正弦強迫。嘗試消除強迫或大幅度降低幅度。此外，它看起來像你有很多阻尼。您的damping ratio，ζ（zeta），似乎是c/(2*w) = 5/4。這意味着你的系統的非受迫形式被過度阻尼 - 沒有強迫你不會看到振盪。

另外，使用振盪系統時要小心使用ode45。如果您的系統剛好太硬，您可能需要調整容差或使用硬解算器（如ode15s）來獲得準確的結果。你總是可以嘗試使用更嚴格的公差來檢查它們是否產生定性相似的結果（相同週期，幅度，增長率/衰減率等）。