Application Center - Maplesoft
  Submit your work
Application Center Applications Visualization Free and Forced Harmonic Oscillations Preview

App Preview:

Visualization Free and Forced Harmonic Oscillations

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application



 

 

 

     Ali_Worksheet :Visualization  Free and Forced Harmonic Oscillations

                                                Using Maple 9

                                   By:  Ali Mohamed Abu Oam

        International University of Africa ( Sudan )                   Faculty of Pure and Applied Science

 

                         aaalll20002000@yahoo.com  ******************   asam295@hotmail.com

 

           Abstract


            Introduction 


         Animation
Maple code for a mass-spring system
Maple code for a rectangular pulse


           Free Harmonic Oscillations 
Undamped Motion
Overdamped Motion
Critically Damped Motion
Underdamped Motion 


          Forced Harmonic Oscillations
The response of a mass-spring system to an impulse function
The displacement resonance in forced oscillations


         

Abstract

     This paper focuses on how the symbolic, numerical and graphical power of the computer algebra system Maple 9 can be used to explore and visualize with animation free harmonic motion with and without damping and vibration of a mass-spring system forced by an external sharp blow. Dirac's delta function is used to model the external force that acts for a very short period of time. It is also demonstrated how displacement resonance in a forced oscillations readily can be obtained from the transfer function of the system represented by a second order differential equation. And the information contained in the system frequency response may also be conveniently displayed in animated graphical form.

 

Introduction

    I have incorporated the use of the computer algebra system Maple 9 in the engineering mathematics curriculum . With the powerful software program including graphical, symbolic and numerical techniques, Maple 9  is the ideal program to make mathematics more relevant and motivating for engineering students and an excellent demonstration tool for the classroom. The use of Maple has changed the way I understood mathematics . I now have the opportunity to emphasize the learning of concepts, the visualizatione of concepts and the solution of more realistic engineering applications.

  Many of the phenomena we observe in everyday life have a periodic motion. And in many cases the periodic motion is maintained by a periodic driving force. For these forced oscillations the amplitude depends on the frequency of the driving force. There is usually some frequency of the driving force at which the oscillations have their maximum amplitude. This occur when the natural frequency of the vibrating system and the frequency of the driving force are approximately equal. This phenomenon of resonance plays an important role in almost every branch of physics. And the avoidance of destructive resonance vibrations is an important factor in the design of mechanical systems of all types.

      Application of differential equations plays an important role in science and engineering. And the most important step in determining the natural frequency of vibration of a system is the formulation of its differential equation. The purpose of this article is to demonstrate how Maple can be used to investigate and visualize with animations:

- simple harmonic motion of a mass-spring system
- the response of a mechanical system (mass-spring system) to an impulse function
- the displacement resonance in forced oscillations

Animation

   Coding in Maple does not require expert programmimg skills, so writing a Maple program deals with putting a proc ( ) and an end around a sequence of commands. And it is pretty easy to put prewitten routines together from Maples powerful bulding blocks and plotting facilities. The plottools package provides many useful commands for producing plotting objects, which can be scaled, translated and rotated. The following Maple codes will be used in the coming sections. The first prosedure animate a spring-mass system with and without a dashpot, and the second animate the response to a rectangular pulse by a linear differential second order equation.


Because it is impossible to show animation on paper, the animation figures shows 4 consecutive frames. The animations in this paper can be seen live on the computer.


Maple code for a mass-spring system

m: mass, r: damping constant, k: spring constant, x1,x2 positions of the spring
x0 = x(0), xp0=D(x)(0): initial conditions, f: driving force, n: number of frames

> 

restart: with(plots): with(plottools): with(linalg): with(LinearAlgebra):

Warning, the name changecoords has been redefined

Warning, the name arrow has been redefined

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the assigned name GramSchmidt now has a global binding

> 

mass_spring:=proc(m,r,k,x1,x2,x0,xp0,tk,f,n )
local mass,deq,init,sol,xk,xu,plt1,plt2,plt,pltxk,rect,base,spring,
spring1,rod,cylinder, piston,fluid,dashpot;
with(plots):with(plottools):
spring:=proc(x1,x2,n) #procedure for the spring
local p1,p2,p3,p4,pn_1,pn,p;
p1:=[x1,0];
p2:=[x1+0.25,0];
p3:=[x2-0.25,0];
p4:=[x2,0];
p:=i->[x1+0.25+(x2-x1-0.5)/n*i,(-1)^(i+1)*0.5];
plot([p1,p2,seq(p(i),i=1..n-1),p3,p4],thickness=2,color=aquamarine);
end:
spring1:=x2->spring(x1,x2,12);
rect:=t->translate(rectangle([-0.2,-0.2],[0.2,0.2],color=red),t,xk(t)+x2+1);
#rect: display the position of the mass on the displacement curve x(t)
mass:=x2->rectangle([x2,-1],[x2+2,1],color=red):
if x0=0 and xk(0.1)=0 then
xu:=1.5
elif x0=0 and xk(0.1)<>0 then
xu:=5.0;
else
xu:=x0;
fi;
rod:=x2->plot([[x2+2,0],[x2+2.5+2*xu,0]],color=grey,thickness=2);
cylinder:=plot([[x2+2.5+xu,-1],[x2+2.5+xu,1],[x2+3.5+3*xu,1],[x2+3.5+3*xu,-1],
[x2+2.5+xu,-1]],color=tan,thickness=2);
piston:=x2->rectangle([x2+2.5+2*xu,-1],[x2+3+2*xu,1],color=grey);
fluid:=rectangle([x2+2.5+xu,-1],[x2+3.5+3*xu,1],color=blue);
dashpot:=x2->display(rod(x2),cylinder,piston(x2),fluid);
deq:=m*diff(x(t),t$2)+r*diff(x(t),t)+k*x(t)=f(t);
print(deq);
init:=x(0)=x0, D(x)(0)=xp0:
sol:=dsolve({deq,init},x(t));
print(combine(simplify(sol)));
xk:=unapply(rhs(sol),t);
#xk=x(t) : the displacement of mass
pltxk:=plot([xk(t)+x2+1,x2+1],t=0..max(tk,x2+4+3*xu)+0.5,
color=[blue,grey],numpoints=400);
base:=plot(-1.1,t=0..max(tk,x2+4+3*xu)+0.5,color=brown,thickness=4);
if r=0 then #undamped case
plt1:=x->translate(display(spring1(x2+x),mass(x2+x),base),0,-3);
else #with damping
plt1:=x->translate(display(spring1(x2+x),mass(x2+x),dashpot(x2+x),base),0,-3);
fi;
plt2:=i->display(pltxk,rect(tk/n*i));
plt:=i->display(plt1(xk(tk/n*i)),plt2(i));
display(seq(plt(i),i=0..n),insequence=true,args[11..nargs]);
end:

Maple code for a rectangular pulse

m: mass, r: damping constant, k: spring constant,
a: constant in the impulse function  , delta(t-a)  
eps: width  epsilon   of a rectangular pulse , n: number of frames

> 

impulse_func:=proc(m,r,k,a,eps,n)
local h,dlign1,dlign2,sol1,sol2,plt1,plt2,pltd,plt,txt,txtd,kloss,e;
with(plots):
alias(u=Heaviside):
h:=eps->(u(t-a)-u(t-a-eps))/eps;
dlign1:=eps->m*diff(y(t),t,t)+r*diff(y(t),t)+k*y(t)=5*h(eps);
dlign2:=m*diff(y(t),t,t)+r*diff(y(t),t)+k*y(t)=5*Dirac(t-a);
print(dlign1(eps), epsilon=eps);
print(dlign2);
sol1:=eps->simplify(dsolve({dlign1(eps),y(0)=0,D(y)(0)=0},y(t)));
sol2:=dsolve({dlign2,y(0)=0,D(y)(0)=0},y(t));
e:= textplot([6,1.3,convert([101],bytes)], font=[SYMBOL,12]):
txt:=eps->textplot([6.9,1.3,cat(` = `,convert(evalf(eps,2),string))]);
txtd:=textplot([4,1.4,`Response to Dirac's delta function`],align={RIGHT,ABOVE});
pltd:=eps->plot(h(eps),t=0..5*Pi,color=green,thickness=2):
plt1:=eps->plot(rhs(sol1(eps)),t=0..5*Pi,color=red,thickness=2):
plt2:=plot(rhs(sol2),t=0..5*Pi,color=blue,thickness=2):
plt:=eps->display(pltd(eps),txt(eps),e,txtd,plt1(eps),plt2);
display(seq(plt(eps-(eps-0.1)/n*i),i=0..n),insequence=true,args[7..nargs]);
end:

Free Harmonic Oscillations


Machines with rotating components commonly involve mass-spring systems or their equivalents in which the driving force is simple harmonic. The motion of a mass attached to a spring serves as a simple example of vibrations that occur in more complex mechanical systems.

It is suitable to consider a body of mass m attached to one end of a spring that resists compression as well as stretching. A rod attached to the mass carries a disk moving in an oil-filled cylinder (a dashpot). The other end of the spring could be attached to a fixed wall and vibrating horizontally. The resultant force on the body is the sum of the restoring force - k*x(t)   and the damping force - r*(diff(x(t), t))   , k    is the force constant,  x(t)    is the distance of the body of mass from its equilibrium position, t  is time and and  r  is the damping constant. We take  x(t)   > 0 when the spring is stretched. The differential equation of motion is therefore:


> 

deq1:=m*diff(x(t),t$2)+r*diff(x(t),t)+k*x(t)=0;

deq1 := m*(diff(x(t), `$`(t, 2)))+r*(diff(x(t), t))+k*x(t) = 0

(1)

If we set r = 0  in deq1    the motion is undamped. Otherwise the solution of  deq1   presents three distinct cases of damping according to whether d = r^2-4*m*k     is greater than, equal to , or less than zero.

> 

eq:=d=r^2-4*m*k;

eq := d = r^2-4*m*k

(2)
> 

solm:=solve(eq,m);

solm := (-d+r^2)/(4*k)

(3)

Let us replace  m   with a new variable m*v    in  deq1  

> 

mv:=unapply(solm,d,r,k);

mv := proc (d, r, k) options operator, arrow; (1/4)*((-d+r^2)/k) end proc

(4)
> 

deq2:=collect(4*k*subs(m=mv(d,r,k),deq1),diff(x(t),t$2));

deq2 := (-d+r^2)*(diff(x(t), `$`(t, 2)))+4*k*(r*(diff(x(t), t))+k*x(t)) = 0

(5)

With this differential equation, deq1   , the system is overdamped if  0 < d   , critically damped if d = 0   and underdamped if d < 0   . If we request Maple to solve deq1    for the undamped case deq2    and  for each of the damped cases,subject to initial condition  x(0) = x[0]   and = 0

> 

init:=x(0)=x[0],D(x)(0)=0: #initial conditions

rwe get:

 Undamped



With      r = 0      N/m    
With  r = 0  

> 

deq1a:=subs(r=0,deq1);

deq1a := m*(diff(x(t), `$`(t, 2)))+k*x(t) = 0

(6)
> 

dsolve({deq1a,init},x(t));

x(t) = x[0]*cos(k^(1/2)*t/m^(1/2))

(7)

A typical graph of  x(t)   with  m = k*g    and  k = 1   N/m    is shown in Figure 1. Figure 2 shows the animation of the undamped motion

> 

xk:=unapply(rhs(%),m,k,x[0],t): # defines x = x(m,k,x[0],t)

 OVERDAMPED

> 

interface(showassumed=0):

> 

assume(d>0):

> 

dsol:=dsolve({deq2,init},x(t));

dsol := x(t) = (1/2)*(x[0]*(r+d^(1/2))*exp(-2*(r-d^(1/2))*k*t/(-d+r^2))/d^(1/2))-(1/2)*(x[0]*(r-d^(1/2))*exp(-2*(r+d^(1/2))*k*t/(-d+r^2))/d^(1/2))

(8)
> 

subs(eq,dsol);

x(t) = (1/2)*(x[0]*(r+(r^2-4*m*k)^(1/2))*exp(-(r-(r^2-4*m*k)^(1/2))*t/(2*m))/(r^2-4*m*k)^(1/2))-(1/2)*(x[0]*(r-(r^2-4*m*k)^(1/2))*exp(-(r+(r^2-4*m*k)^(1/2))*t/(2*m))/(r^2-4*m*k)^(1/2))

(9)
> 

xo:=unapply(rhs(%),m,r,k,x[0],t): # defines x = x(m,r,k,x[0],t)

The solution consists of two exponential terms when   r > {4*m*k}^(1/2)   .  In all subsequent figures we take  m = 1 k*g    and  k   = 1 N/m  . Figure 3 shows some typical graphs of the position function for the overdamped case. Figure 4 shows the animation of the motion with the damping constant   r = 4 N*s/m  

 

CRITICALLY  DAMPED  :

> 

lhs(dsol)=limit(rhs(dsol),d=0);

x(t) = (2*x[0]*k*t+r*x[0])/((exp(k*t/r))^2*r)

(10)
> 

subs(r=sqrt(4*m*k),%);

x(t) = (1/2)*((2*x[0]*k*t+2*(m*k)^(1/2)*x[0])/((exp((1/2)*(k*t/(m*k)^(1/2))))^2*(m*k)^(1/2)))

(11)
> 

xc:=unapply(rhs(%),m,k,x[0],t): #defines x = x(m,k,x[0],t)

The solution consists of one exponential term when r = {4*m*k}^(1/2)    . The graph in Figure 5 resemble those of the overdamped case in Figure 3. The animation of a critically damped motion with  r =  2 N*s/m    is shown in Figure 6.

 

 

UNDER DAMPED :

 

> 

assume(d<0):

 dsolve({deq2,init},x(t));

x(t) = r*x[0]*exp(-2*k*r*t/(-d+r^2))*sin(2*(-d)^(1/2)*k*t/(-d+r^2))/(-d)^(1/2)+x[0]*exp(-2*k*r*t/(-d+r^2))*cos(2*(-d)^(1/2)*k*t/(-d+r^2))

(12)
> 

subs(eq,%);

x(t) = r*x[0]*exp(-r*t/(2*m))*sin((-r^2+4*m*k)^(1/2)*t/(2*m))/(-r^2+4*m*k)^(1/2)+x[0]*exp(-r*t/(2*m))*cos((-r^2+4*m*k)^(1/2)*t/(2*m))

(13)
> 

xu:=unapply(rhs(%),m,r,k,x[0],t): #defines x = x(m,r,k,x[0],t)

The solution represents exponentially damped oscillations of the mass-spring system about its equilibrium position as shown in both Figure 7 and Figure 8.

Undamped Motion

 

> 

plot(xk(1,1,1,t),t=0..19,labels=[`t`,`x(t)`],labelfont=[TIMES,BOLD,12]);

Figure 1 : Undamped motion.  m = 1 kg , k =  1  N/m  ,  x(0) = 2m .

> 

mass_spring(1,0,1,0,5,2,0,14,0,60,scaling=constrained);

Warning, the name arrow has been redefined

Warning, the name arrow has been redefined

diff(x(t), `$`(t, 2))+x(t) = 0

x(t) = 2*cos(t)

Figure 2 Animation of an undamped motion. m = 1kg , k = 1  N/m   ,  x(0) = 2m .

 

The small rectangle shows the position of the mass on the displacement curve.

Overdamped Motion

If we select m = 1kg , k = 1 N/m    and three different d -values: 12,  45,  96 {N*s/m}^2   , r  named  rk becomes:

> 

rk:=(d,m,k)->simplify(sqrt(d+4*m*k)):

 r=[seq(rk(d,1,1),d=[12,45,96])];

r = [4, 7, 10]

(14)
> 

plt1:=d->plot(xo(1,rk(d,1,1),1,2,t),t=0..20):
 txt1:=d->textplot([4,xo(1,rk(d,1,1),1,2,4),convert(r = rk(d,1,1),string)]):
pltod:=display(seq({plt1(d),txt1(d)},d=[12,45,96]),labels=[`t`,`x(t)`],labelfont=[TIMES,BOLD,12]):%;

Figure 3 Overdamped motion, m = 1kg  ,   , k = 1 N/m    , x(0) = 2m : damping constant

Figure 3 shows three typical graphs of the position function  x(t)  for the overdamped case. Figure 4 shows the animation of the motion with  r = 4 N*s/m   .

> 

mass_spring(1,4,1,0,5,2,0,14,0,60,scaling=constrained);

Warning, the name arrow has been redefined

Warning, the name arrow has been redefined

diff(x(t), `$`(t, 2))+4*(diff(x(t), t))+x(t) = 0

x(t) = (2/3)*(exp((-2+3^(1/2))*t)*3^(1/2))+exp((-2+3^(1/2))*t)-(2/3)*(exp((-2-3^(1/2))*t)*3^(1/2))+exp((-2-3^(1/2))*t)

Figure 4  Animation of an overdamped motion of a mass-spring system with dashpot .
                      r = 4 N*s/m   . m = 1kg  ,    k = 1 N/m    ,  x(0) = 2m

 Critically Damped Motion

> 

plt2:=plot(xc(1,1,2,t),t=0..20):
 txt2:=textplot([3,xc(1,1,2,3),convert(r = rk(0,1,1),string)]):
 pltcd:=display(plt2,txt2):%;

Figure 5 Critically damped motion, m = 1kg  ,    k = 1 N/m    ,  x(0) = 2m
 r : damping constant

> 

mass_spring(1,2,1,0,5,2,0,14,0,60,scaling=constrained);

Warning, the name arrow has been redefined

Warning, the name arrow has been redefined

diff(x(t), `$`(t, 2))+2*(diff(x(t), t))+x(t) = 0

x(t) = (2+2*t)*exp(-t)

Figure 6 Animation of an critically damped motion.  m = 1kg  ,    k = 1 N/m    ,  x(0) = 2m
   r = 2 N*s/m   

In this critically damped case, then resistance of the dashpot is just large enough to damp out any oscillations.

Underdamped Motion

We select m = 1kg  ,    k = 1 N/m   ,   and two r values  N*s/(21*5*m)   .

> 

plt3:=plot([xu(1,1/2,1,2,t),xu(1,1/5,1,2,t)],t=0..20):
 txt3:=textplot([ [6.5,xu(1,1/2,1,2,6.5),`r=1/2`],[6.5,xu(1,1/5,1,2,6.5),`r=1/5`]]):
pltod:=display(plt3,txt3,labels=[`t`,`x(t)`],labelfont=[TIMES,BOLD,12]):%;

Figure 7  Underdamped motion ,  m = 1kg  ,    k = 1 N/m    ,  x(0) = 2m
      r : damping constant

> 

mass_spring(1,1/5,1,0,5,2,0,20,0,60,scaling=constrained);

Warning, the name arrow has been redefined

Warning, the name arrow has been redefined

diff(x(t), `$`(t, 2))+(1/5)*(diff(x(t), t))+x(t) = 0

x(t) = (2/33)*(11^(1/2)*exp(-t/10)*sin(3*11^(1/2)*t/10))+2*exp(-t/10)*cos(3*11^(1/2)*t/10)

Figure 8  Animation of an underdamped motion of a mass-spring system with dashpot  m = 1kg  ,    k = 1 N/m  

r = N*s/(5*m)  

(15)

The action of the dashpot exponentially damps the oscillations in accord with  the time-varying amplitude. The dashpot also decreases the frequency of the motion from 1 in the undamped to case  3*{11^(1/2)}/10   in the underdamped motion with the same mass and force constant.


Comparison of the graphs in Figure 2, 5 and 7 (Figure 9) shows that when the motion is critically damped the mass reaches its equilibrium posistion in a shorter time than when the damping is larger. And for any damping constant less than  r  = 2 N*s/m    in our particular example the motion becomes oscillatory.

> 

display(plt1(45),plt2,plt3,txt1(45),txt2,txt3,labels=[`t`,`x(t)`],labelfont=[TIMES,BOLD,12]);

Figure 9  Overdamped, critically damped and damped oscillatary motion of a mass-spring system.
    m = 1kg  ,    k = 1 N/m    ,  x(0) = 2m
       r : damping constant

 

 

Forced Harmonic Oscillations

Mechanical systems are often acted upon by an external force of large magnitude that acts for only a short period of time. For example a spring-mass system could be given a sharp blow at some specific time t . In engineering practice it is convenient to use the Dirac delta function as a mathematical model for such a blow. But it is my experience that it is often difficult for engineering students to get a real understanding of what this impulse function stands for. With the use of Maple it is easy and instructive to simulate Dirac's delta function  delta(t-a)  

The response of a mass-spring system to an impulse function

As said above the Dirac delta function delta(t-a)   could serve as a mathematical model for an external force of large magntiude that acts for only a very short period of time.

> 

delta(t-a)=piecewise(t<a,0,t=a,infinity,t>a,0);

delta(t-a) = piecewise(t < a, 0, t = a, infinity, a < t, 0)

(16)

Where

> 

Int(delta(t-a),t=-infinity..infinity)=int(Dirac(t-a),t=-infinity..infinity);

Int(delta(t-a), t = -infinity .. infinity) = 1

(17)

Obviously no real function can satisfy beeing zero except at a single point and have an integral equal to one.

In Maple delta(t-a)    = Dirac (t - a)  is expressed as the derivative of the Heaviside standard unit step function  u( t-a )

> 

alias(u=Heaviside):
 Diff(u(t-a),t)=diff(u(t-a),t);

Diff(u(t-a), t) = Dirac(a-t)

(18)

It is instructive to use Maple to model such an instantaneous unit impulse by starting with the function

> 

delta[epsilon](t-a)=piecewise(t<a,0,t<a+epsilon,1/epsilon, 0);

delta[epsilon](t-a) = piecewise(t < a, 0, t < a+epsilon, 1/epsilon, 0)

(19)
> 

h:=unapply(convert(rhs(%),Heaviside),t,a,epsilon);

h := proc (t, a, epsilon) options operator, arrow; u(t-a)/epsilon-u(t-a-epsilon)*u(t-a)/epsilon end proc

(20)

A plot of the rectangular pulse for a   = Pi    and epsilon = 3   gives:

> 

plot(h(t,Pi,3),t=0..10,thickness=3);

Figure 10  Rectangular pulse delta[epsilon](t-a)      with   a   = Pi    and  epsilon = 3    (width epsilon   and height 1/epsilon   )

The area of the rectangular pulse is equal to 1 .

> 

assume(epsilon>0):
 A:=Int(h(t,a,epsilon),t=a..a+epsilon):
 A=value(A);

Int(u(t-a)/epsilon-u(t-a-epsilon)*u(t-a)/epsilon, t = a .. a+epsilon) = 1

(21)

where u(t)  is Heavisides unit step function. We can compare the response of a damped mass-spring system to a rectangular pulse delta[epsilon](t-a)    as with the response of the Dirac's delta function used by Maple.

The output of the Maple code impulse_func  in Figure 11 animate the response to the rectangular pulse delta[epsilon](t-a)    by a linear second order differential equation with constant coefficients and with damping.

 

> 

impulse_func(1,3,1,Pi,3,10);

Warning, the name arrow has been redefined

diff(y(t), `$`(t, 2))+3*(diff(y(t), t))+y(t) = (5/3)*u(t-Pi)-(5/3)*u(t-Pi-3), epsilon = 3

diff(y(t), `$`(t, 2))+3*(diff(y(t), t))+y(t) = 5*Dirac(t-Pi)

Figure 11 The behavior of the response to a rectangular pulse by a linear second order

differential equation with damping as compared to the response to Dirac's delta function.


The last animated frame in Figure 11 shows little difference between the two responses when  epsilon = 1  

Figure 12 shows an animated response of the damped mass-spring system initially at rest. At time t = Pi    the system is suddenly given a sharp "hammerblow" modelled by   f(t) = 5*Dirac(t-Pi)  
 

> 

f:=t->5*Dirac(t-Pi):
 mass_spring(1,3,1,0,5,0,0,14,f,60,scaling=constrained);

Warning, the name arrow has been redefined

diff(x(t), `$`(t, 2))+3*(diff(x(t), t))+x(t) = 5*Dirac(t-Pi)

x(t) = Int(5^(1/2)*exp(-(5^(1/2)+3)*t/2)*(Dirac(_z1-Pi)*exp(t*5^(1/2)+(3/2-5^(1/2)/2)*_z1)-Dirac(_z1-Pi)*exp((5^(1/2)+3)*_z1/2)), _z1 = 0 .. t)

 

Figure 12 The motion of a mass-spring system with dashpot

under the influence of a sharp blow at t = Pi    provided by  5*Dirac(t-Pi)  

 

 

The displacement resonance in forced oscillations

In every oscillating system there is dissipation of mechanical energy, which results that the motion of the mass-spring system described in the previous sections dies out. If the oscillation are to be maintained, energy must be supplied to the system. In this section we shall assume that the system is acted on by a periodic driving force. Suppose that the mass-spring system is subjected to a periodic force  F[0] = cos(omega*t)   ,  where F[0]    is the maximum value of the applied force  f = omega/(2*Pi)    and  is its frequency. The equation of motion is then:

                                                          m*(diff(x(t), `$`(t, 2)))+r*(diff(x(t), t))+k*x(t) = F[0]*sin(omega*t)  

where m  is the mass of the system, x(t)  is the distance of the body of mass from its equilibrium position, r is the damping constant and  k  is the force constant. The frequency response of this system can be readily obtained from the system transfer function  H(s) by replacing  s by  i*omega  

> 

f:='f':
 with(inttrans): L:=x->laplace(x,t,s): invL:=X->invlaplace(X,s,t):
 alias(X(s)=L(x(t)),F(s)=L(f(t))):
 dlign:=m*diff(x(t),t$2)+r*diff(x(t),t)+k*x(t)=f(t);

Warning, the name hilbert has been redefined

dlign := m*(diff(x(t), `$`(t, 2)))+r*(diff(x(t), t))+k*x(t) = f(t)

(22)
> 

simplify(L(dlign));

m*s^2*X(s)-m*(D(x))(0)-m*s*x(0)+r*s*X(s)-r*x(0)+k*X(s) = F(s)

(23)
> 

X(s)=solve(%,X(s));

X(s) = (m*(D(x))(0)+m*s*x(0)+r*x(0)+F(s))/(m*s^2+r*s+k)

(24)
> 

subs(x(0)=0,D(x)(0)=0,%);

X(s) = F(s)/(m*s^2+r*s+k)

(25)

The transfer function of the mass-spring system is given by:

                                                                          H(s) = 1/(m*(s^2)+r*s+k)  

> 

H:=s->1/(m*s^2+r*s+k):

It is easy to show that the steady-state frequency response to an input   F(t) = F[0]*cos(omega*t)    becomes:

> 

x(t)=F[0]*abs('H'(I*omega))*cos(omega*t+arg('H'(I*omega)));

x(t) = F[0]*abs(H(omega*I))*cos(omega*t+arg(H(omega*I)))

(26)
> 

value(%);

x(t) = F[0]*cos(omega*t+arg(1/(-m*omega^2+r*omega*I+k)))/abs(-m*omega^2+r*omega*I+k)

(27)

The steady-state system response is also a cosine having the same frequency omega   as the input. And the amplitude of this response is  F[0]*abs(H(omega*I))  . The variation in both the magnitude  abs(H(omega*I))   and argument  arg(H(omega*I)  ) as the frequency  omega   of the input cosine is varied constitute the frequency response of the system, as the following example shows.

Lets us first solve the differential equation using Maples dsolve .

> 

f:=(omega,t)->2*cos(omega*t);

f := proc (omega, t) options operator, arrow; 2*cos(t*omega) end proc

(28)
> 

deq:=(m,r,k,omega)->m*diff(x(t),t$2)+r*diff(x(t),t)+k*x(t)=f(omega,t);

deq := proc (m, r, k, omega) options operator, arrow; m*(diff(x(t), `$`(t, 2)))+r*(diff(x(t), t))+k*x(t) = f(omega, t) end proc

(29)
> 

sol:=dsolve({deq(1,0.25,4,2),x(0)=0,D(x)(0)=0},x(t),method=laplace);

sol := x(t) = -(64/255)*(255^(1/2)*exp(-t/8)*sin(255^(1/2)*t/8))+4*sin(2*t)

(30)
> 

plot(rhs(sol),t=0..40);

                                     Figure 13  The displacement  x(t) of a mass-spring system
                                               undergoing forced oscillations plotted against the time.

 

                                                     m = k*g   ,   r = Float(25, -2)    ,        k = 4*N/m      ,        omega = 2/s    .
                                                                                                                                                                                             The graph shows that the transient solution dies out as  increases.

The transfer function of the system is given by :

> 

H:=(m,r,k,s)->1/(m*s^2+r*s+k);

H := proc (m, r, k, s) options operator, arrow; 1/(m*s^2+r*s+k) end proc

(31)

With  m =1kg  ,    s = omega*I          , and   k = 4*N/m  

 

We get :

> 

'H'(1,r,4,I*omega)=H(1,r,4,I*omega);

H(1, r, 4, omega*I) = 1/(-omega^2+r*omega*I+4)

(32)
> 

Habs:=unapply(evalc(abs(H(1,r,4,I*omega))),r,omega);

Habs := proc (r, omega) options operator, arrow; 1/(((omega^4-8*omega^2+16+r^2*omega^2)^(1/2))) end proc

(33)

There is some frequency called the resonance frequency at which the amplitude  A=2Habs  becomes a maximum.


This resonance frequency can be recognized in many vibrating systems unless the damping force  r  is too large.

> 

solve(omega^4-8*omega^2+16+r^2*omega^2=0,{r});

{r = (omega+2)*(omega-2)*I/omega}, {r = -I*(omega+2)*(omega-2)/omega}

(34)

With zero damping force r = 0 ,  and the resonance frequency omega   = 2   in our example. See also Figure 15.
The phase angle phi    is given by:

> 

phi:=unapply(evalc(argument(H(1,r,4,I*omega))),r,omega);

phi := proc (r, omega) options operator, arrow; arctan(-omega*r, -omega^2+4) end proc

(35)

Substituting the values  r = Float(25, -2)*N*s/m    ,   omega = 2/s     and  F[0] = 2*N   ,  gives the steady-state response,  x*s = x(t)  

> 

xs:=unapply(2*Habs(0.25,2)*cos(2*t+phi(0.25,2)),t):
'xs'(t)=xs(t); #steady state respons

xs(t) = 4.000000000*cos(2*t-1.570796327)

(36)
> 

sol; # Maple's solution

x(t) = -(64/255)*(255^(1/2)*exp(-t/8)*sin(255^(1/2)*t/8))+4*sin(2*t)

(37)

As check ( I done this check):

> 

control:='4*sin(2*t)-4*cos(2*t-Pi/2)'=simplify(4*sin(2*t)-4*cos(2*t-Pi/2));

control := 4*sin(2*t)-4*cos(2*t-Pi/2) = 0

(38)

The system is lagging because the phase shift between input and output  arg(H(iomega  )  is phi   = -1.57 ... .

Figure 14 shows this phase shift.

> 

with(plots):
plt1:=plot([f(2,t),xs(t)],t=30..39,thickness=[1,2]):
plttext1:=textplot([[39,xs(39),`steady-state response`],[39,f(2,39),`input`]],align=RIGHT):
display({plt1,plttext1},labels=[`t`,`x(t)`],labelfont=[TIMES,BOLD,12]);

Warning, the name arrow has been redefined

Figure 14   Phase shift between steady-state response and input displacement  

> 

A:='A':
 Habsi:=i->plot(Habs(i,omega),omega=1..3,A=0..2.5) :
 plt2:=i->display(Habsi(i),text2(i)):
text2:=i->textplot([2,Habs(i,2),convert(r=i,string)],align=ABOVE,font=[TIMES,BOLD,12]):
display(seq(plt2(i),i=[0.1,0.25,0.5,0.75,1,1.41]),labels=[`omega`,`A`],labelfont=[TIMES,BOLD,12]);

Figure 15  Variation of the amplitude  :  A = 2/(omega^4-8*(omega^2)+16+r^2*omega^2)^(1/2) 
  against the angular velocity  omega   of the forced oscillations-displacement resonance for different values of the damping constant r
 .Figure 15  shows that the smaller value of r , the sharper the resonance of the vibrating system to the applied force.

And the resonance frequency varies with the value of the damping constant r . The amplitude A  becomes infinity when the applied frequency becomes equal to the resonance frequency and the system has zero damping.

Figure 16  animates the displacement resonance of the forced oscillations and Figure 17 shows the variation of the phase shift against the angular velocity, which is animated in Figure 18.

Figure 19  shows what happens with the mass-spring system when the damping constant   r = Float(25, -2)*N/s   and Figure 20 when the damping constant is approximately equal to zero, r = Float(1, -2)*N/s    .

> 

display(seq(plt2(0.05*i),i=1..28),insequence=true,labels=[`omega`,`A`],labelfont=[TIMES,BOLD,12]);

> 

Harg:=i->plot(phi(i,omega),omega=1..3) :
text3:=i->textplot([2.4,phi(i,2.4),convert(r=i,string)],align={RIGHT,ABOVE},font=[TIMES,BOLD,12]):
 plt3:=i->display(Harg(i),text3(i)):
display(seq(plt3(i),i=[0.1,0.25,0.5,0.75,1,1.41]),labels=[`omega`,`phi`],labelfont=[TIMES,BOLD,12]);

Figure 17   Variation of the phase shift   phi = arctan(omega^2-1/omega)    against the angular velocity  omega   of the forced oscillations-displacement resonance for different values of the damping constant r
 
As the applied angular velocity is increased from zero to the natural angular velocity of the system, the phase angle  phi   increases from zero to Pi/2   . When  r = 0  the  angle  phi   changes abruptly from  0  to Pi    

> 

display(seq(plt3(0.05*i),i=1..28),labelfont=[TIMES,BOLD,12],insequence=true);

                                                            Figure 18 Animation of Figure 17

> 

f:=t->2*sin(2*t):
 mass_spring(1,0.25,4,0,5,0,0,20,f,60,scaling=constrained);

Warning, the name arrow has been redefined

diff(x(t), `$`(t, 2))+.25*(diff(x(t), t))+4*x(t) = 2*sin(2*t)

x(t) = (4/255)*(255^(1/2)*exp(-t/8)*sin(255^(1/2)*t/8))+4*exp(-t/8)*cos(255^(1/2)*t/8)-4*cos(2*t)

                              Figure 19  The displacement  x(t)  of the mass-spring system
undergoing forced oscillations with damping constant    r = .25                   

> 

mass_spring(1,0.01,4,0,5,0,0,11.5,f,60,scaling=constrained);

Warning, the name arrow has been redefined

Warning, the name arrow has been redefined

diff(x(t), `$`(t, 2))+0.1e-1*(diff(x(t), t))+4*x(t) = 2*sin(2*t)

x(t) = (100/159999)*(exp(-t/200)*sin(159999^(1/2)*t/200)*159999^(1/2))+100*exp(-t/200)*cos(159999^(1/2)*t/200)-100*cos(2*t)

> 

 

> 

 

                                                                          Figure 20  The displacement  x(t)  of the mass-spring system
                                                                                                               undergoing forced oscillations with damping constant     r =0.01 Ns/m

Summary

This paper discusses the use of the Computer Algebra System (CAS) Maple in calculating and animating the responses of free and force vibrations systems. CAS have changed the fundamental way in which many math and engineering courses are taught. The user can enters a mathematical expression, which the CAS (in this paper Maple) subsequently execute and one are freed from tedious computations. This allow the student to explore how the solution of an engineering problem depends on various parameters of the problem under study. The graphics capabilities are extremely helpful for visualizing the behavior of the system under investigation. Maple's plots and plottools packages enable the user to utilize a lot of plotting objects for animating purposes. This is highlighted and demonstrated in this paper. It is hoped that this features can attract the attention of both lectures and students to visualize dynamic response, typical in engineering, more realistically.