%% Mechanical vibrations without forcing
% We have a mass m>0, a spring with spring constant k>0, and a dashpot with damping constant c\ge 0 (with c=0 there is no damping).
% We have the homogeneous ODE
%         my''+cy'+ky=0    or    y''+\frac{c}my'+\frac{k}my=0
% We define the natural frequency \omega_0 = \sqrt\frac{k}m>0 and the damping rate \mu = \frac{c}{2m}\ge 0. So we can write the ODE as
%         y''+2\mu y'+\omega_0^2 y = 0

clearvars; clf
syms omega0 positive
syms mu real
assume(mu>=0)

% We get the characteristic polynomial

syms p(r)              % define symbolic variable r and function p(r)
p(r) = r^2 + 2*mu*r + omega0^2

% which has the two roots

rs = simplify( solve(p(r)==0) , 10 )

% For plotting let us consider the case \omega_0=1:

rs0 = subs(rs,omega0,1)

% Let us plot the real and imaginary parts of r_1,r_2:

f = figure;
f.Position(4) = 2*f.Position(4);      % make figure twice as high
fplot(real(rs0),[0,3]);       ax_origin; hold on
fplot(imag(rs0),[0 3],':');   hold off
xlabel('\mu','FontSize',18)
legend('Re r_1','Re r_2','Im r_1','Im r_2','Location','E','FontSize',12)
title('$r_1,r_2$ (real parts solid, imaginary parts dotted)','Interpreter','latex','FontSize',18)
axis equal

%
% For general \omega_0 we get the same graph, but  the unit on the axes is replaced with \omega_0:

xticks(0:3); yticks(-5:5);
xticklabels({'0','\omega_0','2\omega_0'})
yticklabels({'','-4\omega_0','-3\omega_0','-2\omega_0','-\omega_0','0','\omega_0'});

%
%% Forced vibrations
%% Forced vibrations without damping
% We now have \mu=0. We want to solve the initial value problem
%        y''+\omega_0^2 y=\cos(\omega t),
% \qquad y(0)=0,\quad y'(0)=0
% where \omega is the frequency of the external forcing.

syms y(t)
syms omega positive
Dy = diff(y); D2y = diff(Dy);
ode =  D2y + omega0^2*y == cos(omega*t);
sol =  simplify( dsolve(ode,y(0)==0,Dy(0)==0) )

% This is correct for \omega\ne\omega_0, but Matlab forgot about the case \omega=\omega_0.
% Case \omega=\omega_0:

ode =  D2y + omega0^2*y == cos(omega0*t);
sol = simplify( dsolve(ode,y(0)==0,Dy(0)==0) )

% We now consider the case \omega_0=10. What happens with the solution when the forcing frequency \omega gets closer and closer to the natural frequency \omega_0? In the graphs below we also plot \pm\frac{t}{2\omega_0} (orange lines).

for om = [9.8 9.9 10]
  f = figure;                                      % generate new figure
  f.Position(3) = 2*f.Position(3);                 % make figure twice as wide
  ode1 =  D2y+100*y==cos(om*t);
  sol = simplify( dsolve(ode1,y(0)==0,Dy(0)==0)) 
  fplot(sol,[0 40]); hold on                       % plot solution
  fplot(@(t)[t/20,-t/20],[0,40],'Color',[1 .6 0]); hold off  % plot functions t/(2*omega0),-t/(2*omega0)
  ax_origin; title(sprintf('\\omega = %g',om),'FontSize',20)
  ylim([-2 2]);                                     % use the same y-range for all three plots
end

% Result:
% for \omega close to \omega_0 we obtain so-called "beats": we get oscillations with frequency \frac{\omega+\omega_0}2 multiplied by an amplitude changing with the slow frequency \frac{\omega-\omega_0}2
% for \omega=\omega_0 we have "resonance" and obtain displacements c\cdot t\cdot \sin(\omega_0 t). Now the amplitude keeps growing forever.
%
%% Forced vibrations with damping
% We consider ODE
%         y''+2\mu y'+\omega_0^2 y = \cos(\omega t)      (*)
% with \mu>0.
% For the homogeneous problem we obtain exponentially decaying solutions Y_1(t),\;Y_2(t).
% For the particular solution we use the method of undetermined coefficients and obtain Y(t)=a\cos(\omega t)+b\sin(\omega t), yielding the general solution
%
% y(t) = \underbrace{a \cos(\omega t) + b\sin(\omega t)}+ \underbrace{c_1 Y_1(t) + c_2 Y_2(t) } \\
% \qquad\quad\text{steady state solution}\qquad\text{ transient solution}
%% Finding the particular solution
% We want to find the steady state solution
%         Y(t) = a\cos(\omega t) + b\sin(\omega t) = A\cos(\omega t-\phi)
% We want to understand how the amplitude A and the phase \phi depend on the external frequency \omega.
% We consider a complex valued function u(t) satisfying
%         u''+2\mu u'+\omega_0^2 u = e^{i\omega t}
% Taking the real part of this equation gives y(t)=\text{Re }u(t) satisfying the ODE (*).
% The method of undetermined coefficients gives a particular solution U(t)=Ce^{i\omega t}. Plugging this in gives
%          p(i\omega)C=1 \\
% (-\omega^2+2\mu\cdot i\omega+\omega_0^2)C=1
% yielding
%         C =  \frac1{p(i\omega)}=\frac1{(\omega_0^2-\omega^2)+i\cdot 2\mu\omega}
% For the amplitude we obtain
%        A = |C| = \frac1{|p(i\omega)|}

A = rewrite( 1/abs(p(i*omega)) , 'sqrt')

% Let q = |p(i\omega)|^2. The amplitude A has a maximum where q has a minimum and \partial_\omega q = 0

q = simplify( abs(p(i*omega))^2 ,100)
dq = diff(q,omega)
solution = solve(dq==0,omega,'ReturnConditions',true)
omegamax = solution.omega
solution.conditions

% Result: For small damping \mu<\frac{\omega_0}\sqrt2 the amplitude A has a maximum for external frequency \omega_\max = \sqrt{\omega_0^2-2\mu^2}, and the value of A at the maximum is

Amax = simplify( subs(A,omega,omegamax) )

% Let \omega_0=1. We consider small damping values \mu = \frac{j}8\cdot \frac{\omega_0}\sqrt2 with j=0,\ldots,8 and plot the amplitude A vs. the external frequency \omega.
% The black dots show the location of the maximum.

figure
for d = (0:8)/8*1/sqrt(2)
  fplot(subs(A,{omega0,mu},{1,d}),[0 2]); hold on
  if d~=0
    plot(subs(omegamax,{omega0,mu},{1,d}),subs(Amax,{omega0,mu},{1,d}),'k.','MarkerSize',12)
  end
end
hold off; ylim([0 6]); ax_origin
xlabel('\omega'); title('amplitude A of steady state response for \omega_0=1') 

%
%
% Use ax_origin to have x,y axes through the origin.

function ax_origin
  set(gca,'XAxisLoc','origin','YAxisLoc','origin'); box off
end