Application Center - Maplesoft

# Analysis and Simulation of Simple Dynamic Systems

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

Mechanical Systems.mws

Analysis and Simulation of Simple Dynamic Systems

Columbus Alternative High School

Columbus, Ohio 43211

U.S.A

email : thesoccerstar2k3@hotmail.com

This worksheet demonstrates the use of Maple for studying dynamic systems.

Background and Program Code Usage

Dynamic Systems

Dynamic Systems (as they are used in this worksheet) are mechanical systems comprising of masses and constraints.  Constraints are a set of rules that have to be followed as the masses are allowed to move freely.  For example, a very simple dynamic system might be that of a pendulum.  The constraint of a pendulum is that one mass has to stay at a fixed distance from a fixed mass (rigid-rod constraint). Examples of more complex system are a double pendulum and a triple pendulum.

A Simple Example:
Study of dynamic systems is done by analyzing the differential equations that govern every mass in the system.  For a very simple example, consider only one point object of mass 2 [kg] located at the origin, and the only force acting on it is a 2 [N] force 60 degrees above the x-axis.

This is indeed a very simple system and Newton's second law can be used easily to find the differential equations of motion.  The solutions to these equations will be the functions that represent the motion of the mass object.

The solution to the preceding problem can be found easily using analytic methods, but in order to do so, initial conditions for the system must be defined.  These are the values of functions at initial times needed to successfully integrate the equations and find a particular solution.  The initial conditions needed for each mass in a system are its initial x and y positions (x0, y0) and its initial velocities (vx0 , vy0)   In this case, the initial position is the origin (0,0) and the initial velocity is also (vx0 , vy0) = (0,0) because the object is not given an initial velocity.

When analyzing more complex systems however, the differential equations are not so simple.  The more masses, constraints, vibrations, pin joints etc. are present in the system, the more complex the equations are.  Solving them by hand is nearly impossible, but computers have the ability to find numerical answers given the proper instructions and the equations.

Program Design

This entire project is based on a program code written in Maple 7.  Maple enables solving of differential equations (DE's) in either symbolic or numeric form.  For most simple DE's, Maple is able to find analytic solutions but for complex DE's, Maple uses its numerical algorithms to find a numeric answer.  My objective as a programmer is to set up the differential equations that govern a dynamic system that the user defines through Maple.  Defining can be thought of as describing to the computer a system of masses, constraints and initial conditions, and the job of the computer, then, is to simulate and analyze the motion of the system.

Program Diagram:

Program Usage

The Maple program is made up of modules or procedures that perform specific tasks when called upon.  The user defines a dynamic system using the following commands in the proper format indicated below.  Also shown here are different commands that are used to utilize this program.

1)   (restart/initialize):
Initialize ( ):

2)   (define gravity and viscosity):   GravityViscosity (g, ):
3)   (define a  mass  object):  Mass (n, x0, y0, vx0, vy0, m):
4)   (define a Hooke's spring ):
Spring (A1, A2, k, L):
5)   (add horizontal forceto a mass object):
to a mass object)
7)   (anchor a mass object so:   Fixed (n):
that it remains fixed)

8)   (add constraint that allows a : SlideX (n):
mass object to only slide
horizontally )
9)   (add constraint that allows a :
SlideY (n):
vertically )

10)  (cause a mass object to vibrate:  VibrateX (n, A, )
horizontally)

11)  (cause a mass object to vibrate:   VibrateY (n, A, )
vertically)

12)  (create differential equations CreateEquation s( )
for the specified dynamic system.  Issue this command
after all masses, constraints
etc. have been defined)
13)   solve the DE's that have been:
SolveEquations ( )
created and store the solution in a list format)
14)   (animate the simulation):
RunSimulation (tf, dt )
15)   (visualize graphs of position:
Visualize (v1, v2, tf)

Initializations

For running this program, the following packages are required:

 > restart: with(plots): with(plottools): with(ListTools); with(DEtools):

```Warning, the name changecoords has been redefined

```

```Warning, the name arrow has been redefined

```

```Warning, the assigned name Group now has a global binding

```

```Warning, the name translate has been redefined

```

Program Code

Execute the following code in order for this program to work correctly

 > # 'Initialize all variables and set all flags to 0' module #; Initialize:=proc()  global i,F,ic,m,noo,nos,nof,nslx,nsly,nvibx,nviby,picture;   for i from 1 to 250 do:     F[i]:=0:     ic[i]:=0:     m[i]:=0:     noo:=0:     nos:=0:     nof:=0:     nslx:=0:     nsly:=0:     nvibx:=0:     nviby:=0:     picture:=0:    end do: end proc:

 > # module for defining gravity and air resistance # GravityViscosity:=proc(GG,muu)

 > global g,mu;

 > g:=GG:mu:=muu;

 > end proc:

 > # module for defining mass objects and their initial conditions # Mass:=proc(n,x0,y0,vx0,vy0,m0)    global F,ic,m,noo:

 > F[n]:=array([-mu*diff(x[n](t),t),m0*g-mu*diff(y[n](t),t)]):

 > ic[n]:={x[n](0)=x0,y[n](0)=y0,D(x[n])(0)=vx0,D(y[n])(0)=vy0}:    m[n]:=m0;    noo:=noo+1

 > end proc:

 > # module for defining springs and their properties # Spring:=proc(n1,n2,k,L)

 > global d,F,ic,temp1,temp2,temp3,temp4,springarray,nos: d:=L; nos:=nos+1; springarray[nos]:=array([n1,n2,k,L]);

 > temp1:=(F[n1][1]): temp1:=temp1+(x[n1](t)-x[n2](t))/sqrt((x[n1](t)-x[n2](t))^2+(y[n1](t)-y[n2](t))^2)*(-k)*(sqrt((x[n1](t)-x[n2](t))^2+(y[n1](t)-y[n2](t))^2)-d)-mu*diff(x[n1](t),t): (F[n1][1]):=temp1:

 > temp2:=(F[n1][2]): temp2:=temp2+(y[n1](t)-y[n2](t))/sqrt((x[n1](t)-x[n2](t))^2+(y[n1](t)-y[n2](t))^2)*(-k)*(sqrt((x[n1](t)-x[n2](t))^2+(y[n1](t)-y[n2](t))^2)-d)-mu*diff(y[n1](t),t): (F[n1][2]):=temp2:

 > temp3:=(F[n2][1]): temp3:=temp3+(x[n2](t)-x[n1](t))/sqrt((x[n1](t)-x[n2](t))^2+(y[n1](t)-y[n2](t))^2)*(-k)*(sqrt((x[n1](t)-x[n2](t))^2+(y[n1](t)-y[n2](t))^2)-d)-mu*diff(x[n2](t),t): (F[n2][1]):=temp3:

 > temp4:=(F[n2][2]): temp4:=temp4+(y[n2](t)-y[n1](t))/sqrt((x[n1](t)-x[n2](t))^2+(y[n1](t)-y[n2](t))^2)*(-k)*(sqrt((x[n1](t)-x[n2](t))^2+(y[n1](t)-y[n2](t))^2)-d)-mu*diff(y[n2](t),t): (F[n2][2]):=temp4

 > end proc:

 > # module for adding a horizontal force to a mass object # AddForceX:=proc(n,fx)    global F;    F[n][1]:=F[n][1]+fx; end proc:

 > # module for adding a vertical force to a mass object # AddForceY:=proc(n,fy)    global F:    F[n][2]:=F[n][2]+fy; end proc:

 > # module for anchoring mass objects (fixing them in place) # Fixed:=proc(n0)    global F,pin,nof;    nof:=nof+1;    pin[nof]:=n0:    F[n0][1]:=0:    F[n0][2]:=0: end proc:

 > # module for allowing mass objects to slide horizontally only # SlideX:=proc(n0)    global F,nslx,sliderx;    nslx:=nslx+1;    sliderx[nslx]:=n0;    F[n0][2]:=0; end proc:

 > # module for allowing mass objects to slide vertically only # SlideY:=proc(n0)    global F,nsly,slidery;    nsly:=nsly+1;    slidery[nsly]:=n0;    F[n0][1]:=0; end proc:

 > # module for causing a mass object to vibrate horizontally # VibrateX:=proc(n,amp,om)    global F,nvibx,arrx;    nvibx:=nvibx+1;    arrx[nvibx]:=n;    F[n][1]:=m[n]*amp*om^2*(cos(om*t))/2;    F[n][2]:=0; end proc: # module for causing a mass object to vibrate vertically # VibrateY:=proc(n,amp,om)   global F,nviby,arry;    nviby:=nviby+1;    arry[nviby]:=n;    F[n][2]:=m[n]*amp*om^2*cos(om*t)/2;    F[n][1]:=0; end proc:

 > # module for creating the system differential equations # CreateEquations:=proc()   global eq,i,h,k,H,u;   for i from 1 to noo do;      eq[i]:={m[i]*diff(x[i](t),t,t)=F[i][1],m[i]*diff(y[i](t),t,t)=F[i][2]}:   end do; H:=[seq([eq[i][1],eq[i][2],ic[i][1],ic[i][2],ic[i][3],ic[i][4]],i=1..noo)]:   u:= Flatten(H);    end proc:

 > # module for solving the differential equations for the system # SolveEquations:=proc()    global Sys,eq,sol,i,X,Y,u;             Flatten(H);    u:=Flatten(H):    Sys:=u:    sol:=dsolve(Sys,numeric,output=listprocedure):    for i from 1 to noo do;        X[i]:=eval(x[i](t),sol);        Y[i]:=eval(y[i](t),sol);    end do: end proc:

 > # module for animating the dynamic system # RunSimulation:=proc(tf,delt)    global massobject,springobject1,h,k,picture,sliderobjectx,sliderobjecty,vibobjectx,vibobjecty,pinobject:    local x1,x2,x3,x4,x5,x6,x7,x8,x9,y1,y2,y3,y4,y5,y6,y7,y8,y9,cc,s,i,j,middlenumber,size,a1,a2,dis,z,jj:    cc:=1:    for i from 0 by delt to tf do:

 > for j from 1 to noo do:

 > middlenumber:=trunc((4*noo+1)/2)+1:

 > h[j]:=X[j](i);

 > k[j]:=Y[j](i);

 > massobject[j]:=disk([h[j],k[j]],.07,color=blue):

 > end do:        for z from 1 to nos do:            a1:=springarray[z][1];            a2:=springarray[z][2];            size:=springarray[z][4];            x1:=h[a1];            y1:=k[a1];            x2:=h[a2];            y2:=k[a2];            dis:=sqrt((x2-x1)^2-(y2-y1)^2);            x3:= (x1 + x2)/2: y3 := (y1 + y2)/2:            x4:= (x1 + x3)/2: y4 := (y1 + y3)/2:            x5:= (x3 + x2)/2: y5 := (y3 + y2)/2:            x6:= (x4 + x1)/2: y6 := (y4 + y1)/2:            x7:= (x3 + x4)/2: y7 := (y4 + y3)/2:            x8:= (x3 + x5)/2: y8 := (y3 + y5)/2:            x9:= (x5 + x2)/2: y9 := (y5 + y2)/2:               springobject1[z]:={point([x3,y3],color=red),point([x4,y4],color=red),point([x5,y5],color=red),point([x6,y6],color=red),point([x7,y7],color=red),point([x8,y8],color=red),point([x9,y9],color=red)};      end do:     for j from 1 to nslx do;     jj:=sliderx[j];     sliderobjectx[j]:=line([h[jj]-1/2,k[jj]],[h[jj]+1/2,k[jj]],linestyle=7)     end do:     for j from 1 to nsly do:      jj:=slidery[j];  sliderobjecty[j]:=line([h[jj],k[jj]+1/2],[h[jj],k[jj]-1/2],linestyle=7)     end do:       for j from 1 to nvibx do:      jj:=arrx[j];   vibobjectx[j]:= {line([h[jj]+.15,k[jj]],[h[jj]-.15,k[jj]],color=red), line([h[jj]+.15,k[jj]],[h[jj]+.10,k[jj]+.05],color=red), line([h[jj]+.15,k[jj]],[h[jj]+.10,k[jj]-.05],color=red), line([h[jj]-.15,k[jj]],[h[jj]-.10,k[jj]+.05],color=red), line([h[jj]-.15,k[jj]],[h[jj]-.10,k[jj]-.05],color=red)}    end do:      for j from 1 to nviby do:      jj:=arry[j];   vibobjecty[j]:= {line([h[jj],k[jj]+.15],[h[jj],k[jj]-.15],color=red),  line([h[jj],k[jj]+.15],[h[jj]+.05,k[jj]+.10],color=red),  line([h[jj],k[jj]+.15],[h[jj]-.05,k[jj]+.10],color=red),  line([h[jj],k[jj]-.15],[h[jj]+.05,k[jj]-.10],color=red),  line([h[jj],k[jj]-.15],[h[jj]-.05,k[jj]-.10],color=red)}    end do: for j from 1 to nof do;      jj:=pin[j]; pinobject[j]:={line([h[jj]-.1,k[jj]],[h[jj]+.1,k[jj]],thickness=2,color=red),line([h[jj],k[jj]+.1],[h[jj],k[jj]-.1],thickness=2,color=red)};   end do;

 > picture[cc]:=display(seq(massobject[q],q=1..noo),seq(springobject1[w],w=1..nos),seq(pinobject[w],w=1..nof),seq(sliderobjectx[w],w=1..nslx),seq(sliderobjecty[w],w=1..nsly),seq(vibobjectx[w],w=1..nvibx),seq(vibobjecty[w],w=1..nviby),scaling=constrained):

 > cc:=cc+1;        print("Solving system of equations at time (t)"=i);

 > end do: display(seq(picture[q],q=1..cc-1),scaling=constrained,axes=boxed,insequence=true):

 > end proc:

 > # module for visualizing a mass object's position vs. time and more # Visualize:=proc(a,b,tf)   local temp1,temp2,temp3;   temp3:=Flatten([seq([x[i](t),y[i](t)],i=1..noo)]);     temp1:=Flatten([seq([eq[i][1],eq[i][2]],i=1..noo)]);   temp2:=[Flatten([seq([ic[i][1],ic[i][2],ic[i][3],ic[i][4]],i=1..noo)])];   DEplot(temp1,temp3,t=0..tf,temp2,scene=[a,b],stepsize=.005,thickness=1,linecolour=red);   display(%,axesfont=[HELVETICA,7],labelfont=[HELVETICA,7],scaling=constrained); end proc:

Case Study 1 (Damped Oscillation)

This program is used for studying simple mechanical systems, especially systems that involve vibrations.  Consider a simple system consisting of two mass objects of mass 1 [kg]; one fixed and one free (with mass-numbers 1 and 2 respectively) , connected by a spring of stiffness 50 [N/m] and length one.  The free mass is given an initial velocity downwards of 0.5 [m/s].  Viscosity (air-resistance) is present and the coefficient of dampening is  = 0.3.  Gravity is present and the gravitational constant is the usual -9.8 [m/ ].  How will this system behave?

In the above picture of the model, the masses are blue circles and the spring is the set of dots connecting the masses.  Here is the Maple input for defining this system (the masses are located at (0,0) and (0,-1); the red cross indicates that the mass is anchored or fixed):

 > Initialize():

 > GravityViscosity(-9.8,0.3):

 > Mass(1,0,0,0,0,1): Mass(2,0,-1,0,-.5,1): Spring(1,2,50,1):

 > Fixed(1):

 > CreateEquations():

 > SolveEquations();

After entering these commands, Maple has created the differential equations that govern the motion of the system.  Issuing the

RunSimulation (finaltime, increment) command animates the system from t = 0 to t = finaltime with the given increment. Let's run this simulation upto time = 6 seconds with an increment of .1

 > RunSimulation(6,.1);

This program also allows you to plot the position of a mass object versus time.  Suppose you want to graph the y-position of the bottom mass (mass # 2) versus time upto time = 6 seconds.  This following command is used:

 > Visualize(t,y[2](t),6);

As predicted, the air-resistance causes the vibration of mass-2 to decrease with time.

Case Study 2 (Resonance)

Consider a similar model to the previous one; instead of the top mass being fixed, it is forced to vibrate vertically with an amplitude of .2 meters and a frequency of 6.  Also, the length of the spring is increased to 3 meters.  There is no air resistance.  Here is a visual representation of the model.  The red arrows indicate that the mass is vibrating:

The top mass (vibrating) is mass 1 and the bottom is mass 2.  Their coordinates are (0,0) and (0,-3) respectively.  Here is the Maple input for defining this model:

 > Initialize():

 > GravityViscosity(-9.8,0):

 > Mass(1,0,0,0,0,1): Mass(2,0,-3,0,0,1): Spring(1,2,50,3):

 > VibrateY(1,.2,6):

 > CreateEquations():

 > SolveEquations():

Animate the model until time = 4 seconds with .1 second increments:

 > RunSimulation(4,.1);

Graph the y-position of the bottom mass (mass 2) versus time upto time = 20 seconds:

 > Visualize(t,y[2](t),20);

It can be seen that the amplitude of the vibration of the bottom mass increases, then drecreases and increases etc.  within the 20 second period.  However, if the vibration of the top mass matched the natural freqency of the system, the bottom mass will continue vibrate with increasing amplitude (resonance).  The natural frequency of a spring-mass system is  so in this case that frequency would be  or about 7.07.  Let's try this model again only changing the vibrational frequency of the top mass to match the natural frequency so resonance can occur.

Maple Input:

 > Initialize():

 > GravityViscosity(-9.8,0):

 > Mass(1,0,0,0,0,1): Mass(2,0,-3,0,0,1): Spring(1,2,50,3):

 > VibrateY(1,.2,7.07):

 > CreateEquations():

 > SolveEquations():

Graph the y-position of mass 2 versus time upto time = 8 seconds:

 > Visualize(t,y[2](t),8);

The bottom mass object's vibrational amplitude increses with time: it is resonating.

Consider a mass object with a mass of 0.1 [kg] located at the origin and it is only allowed to slide horizontally as if it were a bead on a frictionless wire.  It is connected to a fixed mass (located at (-1,0)) and a free mass (located at (-1, -1)) by springs of stiffness 1 and 200 [N/m] respectively.  One of the two connections (with the stiffer spring) forms a pendulum in which the free mass is allowed to swing freely.   Our goal is to model this scenario.  Here is a picture of the model:

The masses are labeled for clarity.  The horizontal dotted line through mass-1 indicates that it is only allowed to slide horizontally.  Mass-1 is connected to mass-2 and mass-3.  The first connection (between 1 and 2) is made with a spring of stiffness 5 [N/m] and the second connection (between 1 and 3) is made with a spring of stiffnes 200 [N/m] (stiff springs are used as approximations of rigid rods for modeling a pendulm).  Ther gravitational constant is -9.8 and there is no air resistance (  = 0).

Maple Input:

 > Initialize():

 > GravityViscosity(-9.8,0):

 > Mass(1,0,0,0,0,.1): Mass(2,-1,0,0,0,.1): Mass(3,-1,-1,0,0,.1): Spring(1,2,5,1): Spring(1,3,200,1.414):

 > Fixed(2): SlideX(1):

 > CreateEquations():

 > SolveEquations():

Let's run this simulation until time = 3 seconds with .1 second increments:

 > RunSimulation(4,.1);

Here is a graph of the horizontal vibration of mass-1 versus time upto time = 8 seconds:

 > Visualize(t,x[1](t),8);

Here is a graph of the horizontal vibration of mass-3 (swinging mass) versus its vertical vibration upto time = 3 seconds:

 > Visualize(x[3](t),y[3](t),3);

Case Study 4 (Vibrating Tower)

Consider a system of 8  masses connected by springs. They are aligned vertically (resembling a tower) and the lowest of the masses are forced to vibrate horizontally with a given amplitude and frequency:

All mass objects (in blue) have a mass-number (from 1 to 8) and a mass of 0.1 [kg] and the springs have stiffness 400.  The bottom two masses are 1 & 2; the top 2 masses are 7 & 8.  The bottom two masses (1 & 2) are forced to vibrate at an amplitude of  .2 meters and a frequency of 10.   The goal is to see how much the tower vibrates and how the vibration can be minimized.  This program allows you to model this dynamic system and analyze its motion.  Here is the Maple input:

 > Initialize():

 > GravityViscosity(-10,0):

 > Mass(1,0,0,0,0,.1): Mass(2,-1,0,0,0,.1): Mass(3,0,2,0,0,.1): Mass(4,-1,2,0,0,.1): Mass(5,0,4,0,0,.1): Mass(6,-1,4,0,0,.1): Mass(7,0,6,0,0,.1): Mass(8,-1,6,0,0,.1): Spring(1,3,400,2): Spring(2,4,400,2): Spring(1,4,400,sqrt(5)): Spring(2,3,400,sqrt(5)): Spring(4,3,400,1): Spring(4,6,400,2): Spring(3,5,400,2): Spring(3,6,400,sqrt(5)): Spring(4,5,400,sqrt(5)): Spring(5,6,400,1): Spring(5,7,400,2): Spring(7,8,400,1): Spring(8,6,400,2): Spring(5,8,400,sqrt(5)): Spring(6,7,400,sqrt(5)):

 > VibrateX(1,.2,10): VibrateX(2,.2,10):

 > CreateEquations():

 > SolveEquations():

Let's run this simulation upto time = 3 seconds with increments of .1 and see the motion (this may take a LONG time):

 > RunSimulation(3,.1);

This program also allows you to plot the position of a mass object versus time.  Suppose you want to graph the x-position of the top right mass (mass # 8) versus time upto time = 2 seconds.  This following command is used:

 > Visualize(t,x[8](t),2);

The preceding is a graph of the x-position of the top right mass object versus time upto 2 seconds.  It can be seen that witin that time period, the maximum displacement of the mass object was about .4 meters.  Now you can try different values for stiffness of the springs and try to see which ones cause the tower to vibrate less.  There is a value of stiffness of the springs that, if given the right frequency of vibration of the two bottom masses, will cause the tower to resonate and become unstable!  Its up to you to experiment.

Practice Problem

The preceding examples are only a small part of what this program can do.  This worksheet may be considered as a sum of many worksheets that are present in this website because it allows for the study of not just one case but of infinately many cases.  You can define a system as simple or as complex as you like and this program will model it for you!  The beauty lies within its versatility.

Practice Problem :  Use this worksheet to solve the following problem:

Consider the above model of a bridge.  Number the masses (in blue) in any way you like and connect them as they are connected in the picture by springs of different stiffnesses.  All mass objects have a mass of 0.1 [kg].  The two masses that have gray arrowheads pointing to them are fixed masses.  Define all the masses and all the springs.  Set the gravitational constant to -9.8 and the coefficient of air-resistance to 0.1.  Try different values for the stiffness of the springs and notice which ones lead to a stable bridge and which ones cause the bridge to collapse. (For example, weak springs in the middle of the bridge will probably cause the structure to collapse).  Have patience using this program because animating the preceding model will take time.

I am a high school student and would be honored to hear criticisms of this worksheet from professors, mathematicians and scientists. Also, if you have any questions or comments about this program, please feel free to email me.  I will be glad to know and help you.

Also, I know Maple has brought about Maplets as a tool for integration of the Maple platform with Java applets.  I do not own a copy of Maplets but I know enough about it to say that it would be a great tool to combine with the program in this worksheet.  The user can define masses, springs and other constraints through the Maplets Graphical Interface as input and the program can display the animation of graphs of motion as output.  If anyone is interested in making such a Maplet and wants to use this program as the basis of their project, feel free to contact me.  I think it would be a great application for studying mechanics!

Disclaimer:  While every effort has been made to validate the solutions in this worksheet, Waterloo Maple Inc. and the contributors are not responsible for any errors contained and are not liable for any damages resulting from the use of this material.