Description 

A description of the physical system is given, along with the relevant system 

variables and parameters. 

Concepts 

A short list of the concepts presented in the example. 

Modeling 

The system is modeled using ModelBuilder to create the components and 

their interconnections, and BuildEQs to generate the governing equations. 

Simulation 

The kinematic and dynamic equations are solved, usually numerically, to 

provide a computer simulation of the time response.  Plots and animations 

are used to demonstrate the system motion. 

References 

Papers, journals, and web links are provided for further information on the 

given example. 

Description 

A prismatic joint is used to constrain the mass to translational  motion along the global 

X axis.  The displacement of the center of mass, x(t), is measured from  the position where 

the spring is undeformed.  Acting on the mass is an external force F(t), as well as the 

restoring forces from the spring and damper. 

Concepts 

• Prismatic (slider) joint

• BuildEQs

• Solving linear ordinary differential equations

• Vibrating systems

• Plotting results

Modeling 

Simulation 

Since the governing ODE is linear, we can obtain analytical solutions using Maple's dsolve command. 

Let's start by looking at the unforced response case, i.e. with F(t) = 0. 

Numerical values are set for the system parameters, and the system is released from rest with x=1 at 

time=0.  With m=k=1, critical damping would occur for d=2 (see Thomson and Dahleh).  With d=1/2, the 

dynamic response is an underdamped vibration: 

>	ParaSubs:= [m= 1, d= 1/2, k= 1, F(t)= 0]:

 

>	SysODE:= eval(GetDynEQs()[1],ParaSubs):

 

>	UnforcedSoln:= dsolve({SysODE, x(0)=1, D(x)(0)=0});

 

The dynamic response can be plotted, in order to better visualize the response: 

>	assign(UnforcedSoln);

 

>	plot(x(t),t=0..10,title="Displacement versus time");

 

>	plot(diff(x(t),t),t=0..10,title="Speed versus time");

 

Now look at the response of the mass to a sinusoidal forcing function.  Start by unassigning 

the previous solution for x(t), and and then find the new solution: 

>	unassign('x(t)');

 

>	ParaSubs:= [m= 1, d= 1/2, k= 1, F(t)= sin(omega*t)]:

 

>	SysODE:= subs(op(ParaSubs),GetDynEQs()[1]);

 

>	ForcedSoln:= dsolve({SysODE, x(0)=1, D(x)(0)=0}, x(t));

 

>	assign(ForcedSoln);

 

Notice how the forced response depends upon the frequency ω of the forcing excitation.  One can 

plot the displacement response for low and high frequency excitations: 

>	plot(subs(omega=1,x(t)),t=0..10,title="Displacement versus time (omega=1)");

 

>	plot(subs(omega=10,x(t)),t=0..10,title="Displacement versus time (omega=10)");

 

When the forcing frequency is much higher than the natural frequency, one can see that the response 

is essentially the unforced, underdamped response, with a superimposed high-frequency oscillation. 

>

 

References 

Description 

The particle mass is launched with some initial velocity in the XY plane of motion, 

and gravity acts in the -Y direction. 

Concepts 

• Planar joint

• Multiple degrees of freedom (DOF)

• BuildSimCode

• Plotting results

Modeling 

Simulation 

As in the previous example, the governing ODEs are linear and can be solved analytically 

using dsolve. 

However, to demonstrate the features associated with BuildSimCode, optimized simulation 

code will be generated for this example and numerically integrated within Maple. 

Numerical values are set for the system parameters.  Then, an optimized Maple procedure 

(pXdot) is created, which contains a computational sequence for computing the derivatives 

Xdot needed for subsequent numerical integration by Maple's dsolve/numeric. 

>	ParaSubs:= [m= 1, G=10, Inertia= 1]:

 

>	BuildSimCode(Model,{"Xdot"}, "MapleProc", "Optimize", ParaSubs);

 

The projectile is launched with initial speed v0, at a 30 degree angle above the X axis.  The motion 

is simulated by numerically integrating the optimized procedure (pXdot): 

>

v0:= 10:

 

>	odeICs := array([v0*cos(Pi/6),v0*sin(Pi/6),10,0,0,0]);

 

>	odeVars := convert(convert([Model:-vX],Vector),list);

 

>	tr := time(): Soln := dsolve(numeric, implicit=false, procedure=pXdot, abserr = 1e-8, relerr=1e-7,                start=0, initial=odeICs, output=listprocedure, procvars=odeVars, range=0..1): Time_to_integrate := time()-tr;

 

Once the numerical  results are obtained by dsolve/numeric, they can be plotted using plots[odeplot]: 

>	plots[odeplot](Soln, [[odeVars[4],odeVars[5]]],t=0..1, refine=1, title= "Trajectory of Projectile, Side View", labels = ["X", "Y"], scaling= constrained, view=[0..10,0..2]);

 

>

plots[odeplot](Soln, [GetFrameMotion("r", "CoM", "Ground", "Ground")[1],GetFrameMotion("r", "CoM", "Ground", "Ground")[2],GetFrameMotion("r", "CoM", "Ground", "Ground")[3]],t=0..1, refine=1, labels=["x","y","z"], refine=2, axes=framed, orientation=[-80,20],title= "Trajectory of Projectile, 3D View");

 

Note that the angular speed of the mass remains constant, since there are no torques 

acting on the mass in this example: 

>	plots[odeplot](Soln, [t,odeVars[3]],0..1, refine=1, title= "Angular speed versus time", labels = ["Time", "theta_dot"], view=[0..1,0..20]);

 

>

 

References 

Description 

The robot consists of two rigid bodies connected by revolute joints with parallel axes. 

Thus, this "RR" robot is constrained to the XY plane,  with gravity G acting in the -Y 

direction.  The robot is driven by two motor torques,T1 and T2, acting at the two joints. 

Concepts 

• Serial robot manipulators

• GetFrameMotion

• Forward kinematics

• Inverse dynamics

• Forward dynamics

Modeling 

Simulation 

In this example, we assume that the joint coordinates are known functions of time.  From these functions, 

the trajectory of point P can be calculated --- this is known as forward kinematics.  Alternatively, in an inverse 

kinematic analysis, the trajectory of P is given and the joint coordinates are found by solving the previous 

expressions from GetFrameMotion. 

Once the joint coordinates are known, the motor torques required to produce the motion can be found by 

solving the dynamic equations.  This is known as inverse dynamics.  In a forward dynamics problem, the 

torques are given and the dynamic equations are numerically integrated to get the dynamic response (see 

forward dynamic simulations in the projectile and spinning top examples). 

Start by prescribing the joint motions and numerical values for the system parameters: 

>	beta_1:= t->Pi*t^2; beta_2:= t->Pi*t^2;  # both angles traverse 180 degrees in 1 second. L3:= 1: L4:= 1: L5:= 1: L6:= 1:

 

>	plot([beta_1(t),diff(beta_1(t),t)],t=0..1,legend=["beta_1","beta_1_dot"]);

 

We now have enough information to solve the forward kinematics problem for the trajectory of P, 

which is seen to spiral in towards the first joint axis: 

>	plot([map(eval,xP(t)),map(eval,yP(t)),t=0..1],scaling=constrained);

 

For dynamics problems, we also need values for masses, moments of inertia, and gravity.  Once they 

are known, the torques can be obtained from the dynamic equations.  Since the torques appear linearly 

in the dynamic equations (see vF, above), a linear solver (solve) is all that is needed. 

>	m1:= 1: m2:= 1: I1:= 1: I2:= 1: G:= 10:

 

>	Torque1:= solve(GetDynEQs()[1],T1);

 

>	Torque2:= solve(GetDynEQs()[2],T2);

 

>	plot([Torque1,Torque2],t=0..1,legend=["Torque 1","Torque 2"]);

 

In the first part of the response, centripetal effects are small because of small angular speeds.  Thus, 

the static weight and constant angular acceleration terms dominate.  In the latter part of the response, 

the increasing angular speeds begin to affect the required driving torques through the centripetal 

acceleration terms. 

>

 

References 

• J.J. Craig, Introduction to Robotics: Mechanics and Control, 3rd ed., Pearson Prentice Hall, 2005.

• W. Stadler, Analytical Robotics and Mechatronics, McGraw-Hill, 1995.

• R. Schilling, Fundamentals of Robotics: Analysis and Control, Prentice Hall, 1990.

Description 

The base of the top is assumed to remain stationary at the origin O; this assumption is 

enforced by a spherical joint between the ground and the rigid spinning body.  The 3-D 

rotational motion of the top is represented using 3-1-3 Euler angles:  ζ for precession 

about Z, η for nutation about the rotated x frame, and ξ for spin about the body's axis of 

symmetry.  Gravity acts through the center of mass C, located a distance L from O. 

Concepts 

• 3D rotational motion

• Spherical (ball-and-socket) joint

• Angular velocity components versus Euler angle derivatives

• GetKinTrans

• Euler angle singularity

Modeling 

Simulation 

The BuildSimCode command is used to generate an optimized Maple procedure (pXdot) that 

computes the time derivatives (Xdot) of the state variables.  These nonlinear differential equations 

are numerically integrated using dsolve/numeric, given that the top is released from a horizontal 

position (spin axis parallel to -Y) with a spin rate of 10 rad/s.  Numerical values for the system 

parameters are passed as arguments to BuildSimCode: 

>	BuildSimCode(Model,{"Xdot"}, "MapleProc", "Optimize",[L=1, m=1, A=2, C=1, G=9.8]):

 

>	odeVars := convert(convert(Model:-vX,Vector),list); eulerICs := array([0,0,10,0,Pi/2,0]);

 

>	tr := time(): Soln := dsolve(numeric, implicit=false, procedure=pXdot, abserr = 1e-8, relerr=1e-7,                start=0, initial=eulerICs, output=listprocedure, procvars= odeVars, range=0..10): Time_to_integrate := time()-tr;

 

Note that the simulation is much faster than real time, i.e. the CPU time is significantly less 

than the 10s of simulated time.  As expected, the motion of the top is cuspidal: 

>	plots[odeplot](Soln, [t,odeVars[5]*180/Pi],0..10, refine=1, title= "Nutation (deg) versus Time", labels = ["time (s)", "eta (deg)"]);

 

>

L:=1:

 

>	plots[odeplot](Soln, [Model:-GetFrameMotion("r", "CoM", "Ground", "Ground")[1],Model:-GetFrameMotion("r", "CoM", "Ground", "Ground")[2],Model:-GetFrameMotion("r", "CoM", "Ground", "Ground")[3]],t=0..10, refine=1, labels=["x","y","z"], axes=framed, title="3D Motion of Center of Mass");

 

The numerical simulations are repeated, this time using the governing equations in terms 

of angular velocity components: 

>	BuildSimCode(ModelOmega,{"Xdot"}, "MapleProc", "Optimize",[L=1, m=1, A=2, C=1, G=9.8]):

 

>	odeVars := convert(convert([ModelOmega:-vX],Vector),list); omegaICs := array([0,0,10,0,Pi/2,0]);

 

>	tr := time():

 

>	Soln := dsolve(numeric, implicit=false, procedure=pXdot, abserr = 1e-8, relerr=1e-7,                start=0, initial=omegaICs, output=listprocedure, procvars= odeVars, range=0..10):

 

>	Time_to_integrate := time()-tr;

 

The simulation is usually faster using angular velocity components, since the governing equations are 

more simple.  Of course,  the predicted motion of the top is the same in both cases: 

>

plots[odeplot](Soln, [ModelOmega:-GetFrameMotion("r", "CoM", "Ground", "Ground")[1],ModelOmega:-GetFrameMotion("r", "CoM", "Ground", "Ground")[2],ModelOmega:-GetFrameMotion("r", "CoM", "Ground", "Ground")[3]],t=0..10, refine=1, labels=["x","y","z"], axes=framed, title="3D Motion of Center of Mass");

 

>

 

References 

• H. Goldstein, Classical Mechanics, 2nd ed., Addison-Wesley, 1980.

• L. Meirovitch, Methods of Analytical Dynamics, McGraw-Hill, 1970.

Description 

As shown in the Figure above, the crank has a mass m1 (and length L1), while 

the connecting rod has a mass m2 (and length L2).  Gravity in the -Y direction 

is included, as is the force F acting on the piston.  The equations will be generated 

in the 3 joint coordinates [θ,β,s].  Since the system has only 1 degree of freedom, 

2 algebraic constraint equations will be generated for these 3 coordinates, in 

addition to the 3 dynamic equations. 

Concepts 

• Dependent coordinates

• Kinematic constraint equations

• Differential-algebraic equations

Modeling