Application Center - Maplesoft
  Submit your work
Application Center Applications Mobile Robot Modeling and Simulation Preview

App Preview:

Mobile Robot Modeling and Simulation

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

Learn about Maple
Download Application

 

Image 

Mobile Robot Modeling and Simulation 

© Maplesoft, a division of Waterloo Maple Inc., 2005 

Plot 

Derive the model of a two wheel differential drive mobile robot and simulate its trajectory response to various inputs. A 3-D trajectory animation of the mobile robot has been created (shown above) based on the open loop system response of the derived mobile robot model. Various open loop and closed loop responses of this system have been generated in a separate application example application entitled, Mobile Robot . 

 

The robot has two independent wheels, each connected to a motor. This way the robot can move forward, back, and turn. 

Section Layout 

System Definition 

 

Parameters 

 

Variables 

Model Derivation 

 

Robot Chassis 

 

Robot with Wheels 

 

Robot with Wheel and Motor Dynamics 

Simulation 

Image 

Initialization 

Loading libraries ... done 

Defining Custom Functions ... done 

System Definition 

 

The entire system is divided into three subsystems: Robot Chassis, Wheels, and Motors. 

The parameters and variables for each of the subsystems are defined as follows. 

 

Parameter Definition 

Name 

Value 

Units 

Robot Parameters 

Robot Chassis Radius (r[c]) 

Embedded component 

Unit(m) 

Robot Chassis Mass (m[c]) 

Embedded component 

Unit(kg) 

Robot Chassis Moment of Inertia (I[zz]) 

Embedded component 

Unit(kg*m^2) 

Wheel Parameters 

Wheel Radius (r[w]) 

Embedded component 

Unit(m) 

Wheel Mass (m[w]) 

Embedded component 

Unit(kg) 

Wheel Moment of Inertia (I[w]) 

Embedded component 

Unit(kg*m^2) 

Motor Parameters 

DC Motor Armature Resistance (Ra[m]) 

Embedded component 

Unit(Omega) 

DC Motor Armature Inductance (La[m]) 

Embedded component 

Unit(H) 

DC Motor Torque Constant (KMotor[m]) 

Embedded component 

Unit(N*m/A) 

DC Motor EMF Constant (Kemf[m]) 

Embedded component 

Unit(N*m/A) 

Combining all of the parameters together and assigning it to the variable system_parameters. 

 

{r_c = .25, m_c = 5.0, I_zz = .15625, r_w = .1, m_w = .5, I_w = 0.25e-2, Ra_m = 1.0, La_m = .5, KMotor_m = 0.5e-1, Kemf_m = 0.5e-1}
{r_c = .25, m_c = 5.0, I_zz = .15625, r_w = .1, m_w = .5, I_w = 0.25e-2, Ra_m = 1.0, La_m = .5, KMotor_m = 0.5e-1, Kemf_m = 0.5e-1} 

Variable Definition 

Name 

Description 

Robot State 

x(t), y(t) 

Robot Position (x and y coordinates) 

theta(t) 

Robot Heading (counter clockwise from the x-axis when viewed from top) 

v(t), w(t) 

Robot Linear and Angular Velocity (w.r.t. Body Fixed Frame, positive going forward and about the vertical axis). 

Wheel Variables 

omega[L](t), omega[R](t) 

Wheel Angular Speed (left and right) 

tau[L](t), tau[R](t) 

Wheel Applied Torque (left and right) 

F[L](t), F[R](t) 

Wheel Traction Force (left and right) 

Motor Variables 

i[L](t), i[R](t) 

Motor Current (left and right) 

u[L](t), u[R](t) 

Motor Input Voltage (left and right) 

Model Derivation 

 

The models for each of the three subsystems are derived here in sequence. 

 

Robot Chassis 

First, we derive the model for the robot chassis only. 

Kinematics 

Velocity 

Define a mapping from the robot linear velocity v(t) and angular velocity w(t) to the linear velocity and angular velocity in the world coordinates. 

Vector[column](%id = 152543152) = Vector[column](%id = 145520312) 

Acceleration 

Differentiate the velocities to obtain the accelerations. 

Vector[column](%id = 145114524) = Vector[column](%id = 144793688) 

Dynamics 

Kinetic Energy 

 

The kinetic energy of the robot chassis is the sum of the energy from the linear and angular motion. 

 

1/2*m_c*v(t)^2+1/2*I_zz*w(t)^2 

 

States (Generalized Coordinates) 

 

Define the state rate vector of the robot chassis as its forward velocity and turning rate. 

 

Vector[column](%id = 145698720) 

 

Generalized Force 

 

For the robot chassis, the generalized forces are the forces generated by the wheel acting on the chassis. 

 

Vector[column](%id = 152426740) = Vector[column](%id = 151968704) 

 

Euler-Lagrange Equations 

 

Applying the Euler-Lagrange formulation:
Setting the difference between the time derivative of the partial derivative of the Lagrangian with respect to the generalized rate and the generated force to zero.
In this case, the Lagrangian is just the kinetic energy term (since there is no change in the potential energy in the system). 

 

Vector[column](%id = 153029620) = Vector[column](%id = 153350656) 

 

Solving for diff(v(t), t) and diff(w(t), t) results in the following set of equations. 

 

(Typesetting:-mprintslash)([diff(v(t), t) = (F_R(t)+F_L(t))/m_c], [diff(v(t), t) = (F_R(t)+F_L(t))/m_c]) 

 

(Typesetting:-mprintslash)([diff(w(t), t) = (F_R(t)-F_L(t))*r_c/I_zz], [diff(w(t), t) = (F_R(t)-F_L(t))*r_c/I_zz]) 

Robot Motion Equations 

 

To derive the equations of motion in terms of the Cartesian coordinates, we take the acceleration equation and substitute the Euler-Lagrange equations in  and to it. 

 

(Typesetting:-mprintslash)([diff(x(t), `$`(t, 2)) = -sin(theta(t))*(diff(theta(t), t))*v(t)+cos(theta(t))*(F_R(t)+F_L(t))/m_c], [diff(diff(x(t), t), t) = -sin(theta(t))*(diff(theta(t), t))*v(t)+cos(th...
(Typesetting:-mprintslash)([diff(x(t), `$`(t, 2)) = -sin(theta(t))*(diff(theta(t), t))*v(t)+cos(theta(t))*(F_R(t)+F_L(t))/m_c], [diff(diff(x(t), t), t) = -sin(theta(t))*(diff(theta(t), t))*v(t)+cos(th... 

 

(Typesetting:-mprintslash)([diff(y(t), `$`(t, 2)) = cos(theta(t))*(diff(theta(t), t))*v(t)+sin(theta(t))*(F_R(t)+F_L(t))/m_c], [diff(diff(y(t), t), t) = cos(theta(t))*(diff(theta(t), t))*v(t)+sin(thet...
(Typesetting:-mprintslash)([diff(y(t), `$`(t, 2)) = cos(theta(t))*(diff(theta(t), t))*v(t)+sin(theta(t))*(F_R(t)+F_L(t))/m_c], [diff(diff(y(t), t), t) = cos(theta(t))*(diff(theta(t), t))*v(t)+sin(thet... 

 

(Typesetting:-mprintslash)([diff(theta(t), `$`(t, 2)) = (F_R(t)-F_L(t))*r_c/I_zz], [diff(diff(theta(t), t), t) = (F_R(t)-F_L(t))*r_c/I_zz]) 

System Model For Robot Chassis Only 

The system equations for the robot chassis only is formed by combining the Cartesian equation of motion (kinematics) in , and and the Lagrange equations (dynamics) in and . 

 

(Typesetting:-mprintslash)([{diff(v(t), t) = (F_R(t)+F_L(t))/m_c, diff(w(t), t) = (F_R(t)-F_L(t))*r_c/I_zz, diff(x(t), `$`(t, 2)) = -sin(theta(t))*(diff(theta(t), t))*v(t)+cos(theta(t))*(F_R(t)+F_L(t)...
(Typesetting:-mprintslash)([{diff(v(t), t) = (F_R(t)+F_L(t))/m_c, diff(w(t), t) = (F_R(t)-F_L(t))*r_c/I_zz, diff(x(t), `$`(t, 2)) = -sin(theta(t))*(diff(theta(t), t))*v(t)+cos(theta(t))*(F_R(t)+F_L(t)...
(Typesetting:-mprintslash)([{diff(v(t), t) = (F_R(t)+F_L(t))/m_c, diff(w(t), t) = (F_R(t)-F_L(t))*r_c/I_zz, diff(x(t), `$`(t, 2)) = -sin(theta(t))*(diff(theta(t), t))*v(t)+cos(theta(t))*(F_R(t)+F_L(t)... 

 

For our system, we assume the system has zero initial conditions. 

 

{x(0) = 0, y(0) = 0, theta(0) = 0, v(0) = 0, w(0) = 0, (D(x))(0) = 0, (D(y))(0) = 0, (D(theta))(0) = 0} 

Robot With Wheels 

Next, we would like to add the wheel dynamics to the robot chassis to obtain the combined system dynamics. 

Kinematics 

Velocity 

 

Define the kinematics relationship between the linear and angular velocities and accelerations of the wheels. 

 

Vector[column](%id = 152337168) = Vector[column](%id = 151424224) 

 

Expressing as a set of differential equations. 

 

{v_L(t) = r_w*omega_L(t), v_R(t) = r_w*omega_R(t)} 

 

Acceleration 

 

Differentiate the velocities to obtain the accelerations. 

 

Vector[column](%id = 147104360) = Vector[column](%id = 145643852) 

 

Expressing as a set of differential equations. 

 

(Typesetting:-mprintslash)([{diff(v_L(t), t) = r_w*(diff(omega_L(t), t)), diff(v_R(t), t) = r_w*(diff(omega_R(t), t))}], [{diff(v_L(t), t) = r_w*(diff(omega_L(t), t)), diff(v_R(t), t) = r_w*(diff(omeg...
(Typesetting:-mprintslash)([{diff(v_L(t), t) = r_w*(diff(omega_L(t), t)), diff(v_R(t), t) = r_w*(diff(omega_R(t), t))}], [{diff(v_L(t), t) = r_w*(diff(omega_L(t), t)), diff(v_R(t), t) = r_w*(diff(omeg... 

 

Robot-Wheel Connection 

Velocity Relation 

 

Vector[column](%id = 145550288) = Vector[column](%id = 145525220) 

 

Acceleration Relation 

 

Vector[column](%id = 153228908) = Vector[column](%id = 152942636) 

 

Dynamics 

Given the kinematic relationships between the wheel and the chassis, we can now derive the combined system dynamics. 

Kinetic Energy 

The total system kinetic energy is the sum of kinetic energy of the chassis and the wheels. 

 

(1/2*I_zz*r_w^2/r_c^2+1/8*m_c*r_w^2+1/2*m_w*r_w^2+1/2*I_w)*omega_L(t)^2+(1/4*m_c*r_w^2-I_zz*r_w^2/r_c^2)*omega_R(t)*omega_L(t)+(1/2*I_zz*r_w^2/r_c^2+1/8*m_c*r_w^2+1/2*m_w*r_w^2+1/2*I_w)*omega_R(t)^2
(1/2*I_zz*r_w^2/r_c^2+1/8*m_c*r_w^2+1/2*m_w*r_w^2+1/2*I_w)*omega_L(t)^2+(1/4*m_c*r_w^2-I_zz*r_w^2/r_c^2)*omega_R(t)*omega_L(t)+(1/2*I_zz*r_w^2/r_c^2+1/8*m_c*r_w^2+1/2*m_w*r_w^2+1/2*I_w)*omega_R(t)^2 

States (Generalized Coordinates) 

Define the state rate vector of the combined system as the angular rates of the two wheels. 

 

Vector[column](%id = 152735072) 

Generalized Forces 

For the combined system, the generalized forces are the applied torques to the wheels. 

 

Vector[column](%id = 145540804) 

 

Euler-Lagrange Equations 

Apply the Euler-Lagrange Formulation to the combined system kinetic energy equation, we get: 

 

Vector[column](%id = 152146748) = Vector[column](%id = 147078532)
Vector[column](%id = 152146748) = Vector[column](%id = 147078532)
Vector[column](%id = 152146748) = Vector[column](%id = 147078532)
Vector[column](%id = 152146748) = Vector[column](%id = 147078532)
Vector[column](%id = 152146748) = Vector[column](%id = 147078532) 

 

Writing them as separate differential equations: 

 

(Typesetting:-mprintslash)([diff(omega_L(t), t) = 1/2*(r_c^2*m_c*r_w^2*tau_L(t)-r_c^2*m_c*r_w^2*tau_R(t)+4*r_c^2*m_w*r_w^2*tau_L(t)+4*r_c^2*I_w*tau_L(t)+4*I_zz*r_w^2*tau_L(t)+4*I_zz*r_w^2*tau_R(t))/(r...
(Typesetting:-mprintslash)([diff(omega_L(t), t) = 1/2*(r_c^2*m_c*r_w^2*tau_L(t)-r_c^2*m_c*r_w^2*tau_R(t)+4*r_c^2*m_w*r_w^2*tau_L(t)+4*r_c^2*I_w*tau_L(t)+4*I_zz*r_w^2*tau_L(t)+4*I_zz*r_w^2*tau_R(t))/(r...
(Typesetting:-mprintslash)([diff(omega_L(t), t) = 1/2*(r_c^2*m_c*r_w^2*tau_L(t)-r_c^2*m_c*r_w^2*tau_R(t)+4*r_c^2*m_w*r_w^2*tau_L(t)+4*r_c^2*I_w*tau_L(t)+4*I_zz*r_w^2*tau_L(t)+4*I_zz*r_w^2*tau_R(t))/(r...
(Typesetting:-mprintslash)([diff(omega_L(t), t) = 1/2*(r_c^2*m_c*r_w^2*tau_L(t)-r_c^2*m_c*r_w^2*tau_R(t)+4*r_c^2*m_w*r_w^2*tau_L(t)+4*r_c^2*I_w*tau_L(t)+4*I_zz*r_w^2*tau_L(t)+4*I_zz*r_w^2*tau_R(t))/(r... 

 

(Typesetting:-mprintslash)([diff(omega_R(t), t) = 1/2*(4*I_zz*r_w^2*tau_L(t)+4*I_zz*r_w^2*tau_R(t)-r_c^2*m_c*r_w^2*tau_L(t)+4*I_w*r_c^2*tau_R(t)+4*m_w*r_w^2*r_c^2*tau_R(t)+r_c^2*m_c*r_w^2*tau_R(t))/(r...
(Typesetting:-mprintslash)([diff(omega_R(t), t) = 1/2*(4*I_zz*r_w^2*tau_L(t)+4*I_zz*r_w^2*tau_R(t)-r_c^2*m_c*r_w^2*tau_L(t)+4*I_w*r_c^2*tau_R(t)+4*m_w*r_w^2*r_c^2*tau_R(t)+r_c^2*m_c*r_w^2*tau_R(t))/(r...
(Typesetting:-mprintslash)([diff(omega_R(t), t) = 1/2*(4*I_zz*r_w^2*tau_L(t)+4*I_zz*r_w^2*tau_R(t)-r_c^2*m_c*r_w^2*tau_L(t)+4*I_w*r_c^2*tau_R(t)+4*m_w*r_w^2*r_c^2*tau_R(t)+r_c^2*m_c*r_w^2*tau_R(t))/(r...
(Typesetting:-mprintslash)([diff(omega_R(t), t) = 1/2*(4*I_zz*r_w^2*tau_L(t)+4*I_zz*r_w^2*tau_R(t)-r_c^2*m_c*r_w^2*tau_L(t)+4*I_w*r_c^2*tau_R(t)+4*m_w*r_w^2*r_c^2*tau_R(t)+r_c^2*m_c*r_w^2*tau_R(t))/(r... 

 

System Model For Robot with Wheels 

The combined system model is formed by combining the Cartesian equation of motion (kinematics) of the robot chassis in , the acceleration coupling equations between the robot and the wheels in and the Lagrange equations (dynamics) of the combined system and : 

(Typesetting:-mprintslash)([{diff(v(t), t) = 1/2*r_w*(diff(omega_L(t), t))+1/2*r_w*(diff(omega_R(t), t)), diff(w(t), t) = (r_w*(diff(omega_R(t), t))-r_w*(diff(omega_L(t), t)))/r_c, diff(omega_L(t), t)...
(Typesetting:-mprintslash)([{diff(v(t), t) = 1/2*r_w*(diff(omega_L(t), t))+1/2*r_w*(diff(omega_R(t), t)), diff(w(t), t) = (r_w*(diff(omega_R(t), t))-r_w*(diff(omega_L(t), t)))/r_c, diff(omega_L(t), t)...
(Typesetting:-mprintslash)([{diff(v(t), t) = 1/2*r_w*(diff(omega_L(t), t))+1/2*r_w*(diff(omega_R(t), t)), diff(w(t), t) = (r_w*(diff(omega_R(t), t))-r_w*(diff(omega_L(t), t)))/r_c, diff(omega_L(t), t)...
(Typesetting:-mprintslash)([{diff(v(t), t) = 1/2*r_w*(diff(omega_L(t), t))+1/2*r_w*(diff(omega_R(t), t)), diff(w(t), t) = (r_w*(diff(omega_R(t), t))-r_w*(diff(omega_L(t), t)))/r_c, diff(omega_L(t), t)...
(Typesetting:-mprintslash)([{diff(v(t), t) = 1/2*r_w*(diff(omega_L(t), t))+1/2*r_w*(diff(omega_R(t), t)), diff(w(t), t) = (r_w*(diff(omega_R(t), t))-r_w*(diff(omega_L(t), t)))/r_c, diff(omega_L(t), t)...
(Typesetting:-mprintslash)([{diff(v(t), t) = 1/2*r_w*(diff(omega_L(t), t))+1/2*r_w*(diff(omega_R(t), t)), diff(w(t), t) = (r_w*(diff(omega_R(t), t))-r_w*(diff(omega_L(t), t)))/r_c, diff(omega_L(t), t)...
(Typesetting:-mprintslash)([{diff(v(t), t) = 1/2*r_w*(diff(omega_L(t), t))+1/2*r_w*(diff(omega_R(t), t)), diff(w(t), t) = (r_w*(diff(omega_R(t), t))-r_w*(diff(omega_L(t), t)))/r_c, diff(omega_L(t), t)...
(Typesetting:-mprintslash)([{diff(v(t), t) = 1/2*r_w*(diff(omega_L(t), t))+1/2*r_w*(diff(omega_R(t), t)), diff(w(t), t) = (r_w*(diff(omega_R(t), t))-r_w*(diff(omega_L(t), t)))/r_c, diff(omega_L(t), t)...
(Typesetting:-mprintslash)([{diff(v(t), t) = 1/2*r_w*(diff(omega_L(t), t))+1/2*r_w*(diff(omega_R(t), t)), diff(w(t), t) = (r_w*(diff(omega_R(t), t))-r_w*(diff(omega_L(t), t)))/r_c, diff(omega_L(t), t)... 

 

As with the robot chassis, the initial conditions for the wheel are also zero: 

{omega_L(0) = 0., omega_R(0) = 0.} 

 

Robot with Wheel and Motor Dynamics 

Finally, we include the motor dynamics into the whole system model. 

Motor Characteristics 

Applied Torque 

The motors connect to the robot through the applied torques to the wheels. 

 

Vector[column](%id = 146033196) = Vector[column](%id = 145637048) 

 

Expressing as a set of equations. 

 

{tau_L(t) = KMotor_m*i_L(t), tau_R(t) = KMotor_m*i_R(t)}
{tau_L(t) = KMotor_m*i_L(t), tau_R(t) = KMotor_m*i_R(t)} 

 

Dynamics 

The dynamics of the each motor is taken as a linear second order system in the motor current. 

 

(Typesetting:-mprintslash)([diff(i_L(t), t) = -(Ra_m*i_L(t)-u_L(t)+Kemf_m*omega_L(t))/La_m], [diff(i_L(t), t) = -(Ra_m*i_L(t)-u_L(t)+Kemf_m*omega_L(t))/La_m]) 

 

(Typesetting:-mprintslash)([diff(i_R(t), t) = -(Ra_m*i_R(t)-u_R(t)+Kemf_m*omega_R(t))/La_m], [diff(i_R(t), t) = -(Ra_m*i_R(t)-u_R(t)+Kemf_m*omega_R(t))/La_m]) 

System Model For Robot with Wheel and Motor Dynamics 

The combined system model is formed by combining the robot-wheel system equations in , the torque coupling equations in between and the motor dynamics equations in and : 

(Typesetting:-mprintslash)([{diff(x(t), `$`(t, 2)) = -sin(theta(t))*(diff(theta(t), t))*v(t)+cos(theta(t))*(diff(v(t), t)), diff(y(t), `$`(t, 2)) = cos(theta(t))*(diff(theta(t), t))*v(t)+sin(theta(t))...
(Typesetting:-mprintslash)([{diff(x(t), `$`(t, 2)) = -sin(theta(t))*(diff(theta(t), t))*v(t)+cos(theta(t))*(diff(v(t), t)), diff(y(t), `$`(t, 2)) = cos(theta(t))*(diff(theta(t), t))*v(t)+sin(theta(t))...
(Typesetting:-mprintslash)([{diff(x(t), `$`(t, 2)) = -sin(theta(t))*(diff(theta(t), t))*v(t)+cos(theta(t))*(diff(v(t), t)), diff(y(t), `$`(t, 2)) = cos(theta(t))*(diff(theta(t), t))*v(t)+sin(theta(t))...
(Typesetting:-mprintslash)([{diff(x(t), `$`(t, 2)) = -sin(theta(t))*(diff(theta(t), t))*v(t)+cos(theta(t))*(diff(v(t), t)), diff(y(t), `$`(t, 2)) = cos(theta(t))*(diff(theta(t), t))*v(t)+sin(theta(t))...
(Typesetting:-mprintslash)([{diff(x(t), `$`(t, 2)) = -sin(theta(t))*(diff(theta(t), t))*v(t)+cos(theta(t))*(diff(v(t), t)), diff(y(t), `$`(t, 2)) = cos(theta(t))*(diff(theta(t), t))*v(t)+sin(theta(t))...
(Typesetting:-mprintslash)([{diff(x(t), `$`(t, 2)) = -sin(theta(t))*(diff(theta(t), t))*v(t)+cos(theta(t))*(diff(v(t), t)), diff(y(t), `$`(t, 2)) = cos(theta(t))*(diff(theta(t), t))*v(t)+sin(theta(t))...
(Typesetting:-mprintslash)([{diff(x(t), `$`(t, 2)) = -sin(theta(t))*(diff(theta(t), t))*v(t)+cos(theta(t))*(diff(v(t), t)), diff(y(t), `$`(t, 2)) = cos(theta(t))*(diff(theta(t), t))*v(t)+sin(theta(t))...
(Typesetting:-mprintslash)([{diff(x(t), `$`(t, 2)) = -sin(theta(t))*(diff(theta(t), t))*v(t)+cos(theta(t))*(diff(v(t), t)), diff(y(t), `$`(t, 2)) = cos(theta(t))*(diff(theta(t), t))*v(t)+sin(theta(t))...
(Typesetting:-mprintslash)([{diff(x(t), `$`(t, 2)) = -sin(theta(t))*(diff(theta(t), t))*v(t)+cos(theta(t))*(diff(v(t), t)), diff(y(t), `$`(t, 2)) = cos(theta(t))*(diff(theta(t), t))*v(t)+sin(theta(t))...
(Typesetting:-mprintslash)([{diff(x(t), `$`(t, 2)) = -sin(theta(t))*(diff(theta(t), t))*v(t)+cos(theta(t))*(diff(v(t), t)), diff(y(t), `$`(t, 2)) = cos(theta(t))*(diff(theta(t), t))*v(t)+sin(theta(t))...
(Typesetting:-mprintslash)([{diff(x(t), `$`(t, 2)) = -sin(theta(t))*(diff(theta(t), t))*v(t)+cos(theta(t))*(diff(v(t), t)), diff(y(t), `$`(t, 2)) = cos(theta(t))*(diff(theta(t), t))*v(t)+sin(theta(t))... 

 

The initial conditions for the motors are also set to zero. 

{i_L(0) = 0., i_R(0) = 0.} 

 

Simulation 

System Input 

The inputs to the combined system are the applied voltage to the motor. 

For simulation purpose, they are defined as piecewise continuous functions: 

 

 

 

 

Motor Input Voltage 

 

 

Plot 

 

System Responses 

Robot X and Y Position 

 

 

Plot 

 

 

Robot Heading 

 

 

Plot 

 

Robot Chassis Linear and Angular Velocities 

 

Plot 

 

Trajectory Animation 

 

Plot 

 

 

Legal Notice: The copyright for this application is owned by Maplesoft. The application is intended to demonstrate the use of Maple to solve a particular problem. It has been made available for product evaluation purposes only and may not be used in any other context without the express permission of Maplesoft.   

Image