減衰強制振動の過渡応答

clear; close all;
m  = 1;
k = 1;
wn = sqrt(k/m);
wd = wn;
c = 0.1;
F = 1;
zeta = c/(2*sqrt(m*k));
beta = zeta*wn;
t = linspace(0,100,500);

ws = [1,0.2,5];
for idx=1:length(ws)
  w = ws(idx);
A = 1/wd*(m*w^2-k+beta*c)*w/((-m*w^2+k^2)^2+c^2*w^2)*F;
B = c*w/((-m*w^2+k)^2+c^2*w^2)*F;
C = (-m*w^2+k)/((-m*w^2+k)^2+c^2*w^2)*F;
D = -c*w/((-m*w^2+k)^2+c^2*w^2)*F;

x = exp(-beta*t).*(A*sin(wd*t)+B*cos(wd*t))+...
    C*sin(w*t)+D*cos(w*t);

subplot(3,1,idx);
plot(t,x,'DisplayName',sprintf('omega=%.1g',w));
legend();

end

big;