explanation of ode45 for constant force

2 次查看(过去 30 天)
I am implementing the shm equation ( m x'' + c x' + k x = F0) in matlab using ode45. The RHS is a constant force (F0). I am still getting an oscillatory response. The steady state response is equal to the static deformation. What is confusing me are the oscillations. What a force with constant magnitude result in an oscillatory response. Is my implementation wrong or my understanding of physics?
The sample code is shown below.
function dydt = sdof_fun(t,y,m,k,c,F)
dydt(1) = y(2);
dydt(2) = (F/m) - (k/m)*y(1) - (c/m)*y(2);
dydt = dydt(:);
end
[t,w] = ode45(@(t,y) sdof_fun(t,y,m,k,c,F),tspan,ic);
  2 个评论
Ajinkya
Ajinkya 2023-1-23
Oh, Is that the free vibration response (response of the Homogeneous equation)?
Torsten
Torsten 2023-1-23
编辑:Torsten 2023-1-23
It depends on the zeros of the characteristic equation
m*x^2 + c*x + k = 0
what kind of response you get (especially the sign of the discriminant disc = c^2 - 4*k*m).
So check your constants.

请先登录,再进行评论。

回答(1 个)

Sam Chak
Sam Chak 2023-1-24
编辑:Sam Chak 2023-1-24
Analysis shows that the steady-state value of of the linear system is given by
.
In your case, I guess that the oscillatory response is observed because . If so, then try extending the simulation time long enough until you see the response converges to .
Another way is to check the value of the system's damping ratio, given by
.
If the damping ratio value is very small, then the oscillations die out very slowly. If , then a critically-damped response will be observed.
The following example shows the convergence, with so that the oscillations die out faster.
m = 11; % change to 49/20 to observe a critically-damped response
c = 7;
k = 5;
F = 3;
params = [m c k F];
tspan = 0:0.01:30;
ic = [-1 0];
[t, y] = ode45(@(t, y) sdof_fun(t, y, params), tspan, ic);
plot(t, y(:,1), 'linewidth', 1.5), grid on
xlabel('t'), ylabel('y')
title('System Response')
function dydt = sdof_fun(t, y, params)
m = params(1);
c = params(2);
k = params(3);
F = params(4);
dydt(1) = y(2);
dydt(2) = (F/m) - (k/m)*y(1) - (c/m)*y(2);
dydt = dydt(:);
end

类别

Help CenterFile Exchange 中查找有关 Programming 的更多信息

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by