Mechanical vibrations without forcing

We have a mass , a spring with spring constant , and a dashpot with damping constant (with there is no damping).

or

We define the natural frequency and the damping rate . So we can write the ODE as

syms p(r) % define symbolic variable r and function p(r)

p(r) = r^2 + 2*mu*r + omega0^2

p(r) = 

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

rs = 

For plotting let us consider the case :

rs0 = subs(rs,omega0,1)

rs0 = 

Let us plot the real and imaginary parts of :

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 we get the same graph, but the unit on the axes is replaced with :

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 . We want to solve the initial value problem

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) )

sol = 

This is correct for , but Matlab forgot about the case .

Case :

ode = D2y + omega0^2*y == cos(omega0*t);

sol = simplify( dsolve(ode,y(0)==0,Dy(0)==0) )

sol = 

We now consider the case . What happens with the solution when the forcing frequency ω gets closer and closer to the natural frequency ? In the graphs below we also plot (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

sol = 

sol = 

sol = 

• for ω close to we obtain so-called "beats": we get oscillations with frequency multiplied by an amplitude changing with the slow frequency
• for we have "resonance" and obtain displacements . Now the amplitude keeps growing forever.

Forced vibrations with damping

with .

For the homogeneous problem we obtain exponentially decaying solutions .

For the particular solution we use the method of undetermined coefficients and obtain , yielding the general solution

Finding the particular solution

We consider a complex valued function satisfying

Taking the real part of this equation gives satisfying the ODE .

The method of undetermined coefficients gives a particular solution . Plugging this in gives

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

A = 

Let . The amplitude A has a maximum where q has a minimum and

q = simplify( abs(p(i*omega))^2 ,100)

q = 

dq = diff(q,omega)

dq = 

solution = solve(dq==0,omega,'ReturnConditions',true)

solution = struct with fields:

omega: (omega0^2 - 2*mu^2)^(1/2) parameters: [1×0 sym] conditions: 2*mu^2 < omega0^2

omegamax = solution.omega

omegamax = 

solution.conditions

ans = 

Result: For small damping the amplitude A has a maximum for external frequency , and the value of A at the maximum is

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

Amax = 

Let . We consider small damping values with and plot the amplitude A vs. the external frequency ω.

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 axes through the origin.