Name: Ceferino Kevin A.
Tan
EE 179.1 Section: M45
Course and Year: BSECE - 4
Laboratory Schedule: Monday 12-3
Activity 2: Mathematical Modeling of Physical
Systems
Abstract
Obviously, the operation of physical systems often involve a lot of
nonlinearity and unpredictability. We normally refer to these uncertainties as
nonidealities. Understanding the effects of these nonidealities can contribute
to significant improvement in physical systems in terms of accuracy and
long-term reliability.
To create effective machines or physical systems, understanding the physical
nonidealities and how we can mitigate these will yield much more accurate
and much more efficient systems. We then create mathematical models of
these nonidealities to take account of their effects in our physical systems.
We normally take account of the nonidealities of physical systems into our
calculations by the use of the ordinary differential equations, which take
modeled nonidealities, say, the coriolis effect-on a long-range bullet, or the
spring constant of a nonideal spring.
In this activity, we will make use of the capabilities of SciLab to perform
mathematical operations involving ordinary differential equations, which will
be crucial in building effective control systems of physical systems and/or
machinery.
1
Objectives
To grasp the important role of mathematical models of physical
systems in the design and analysis of control systems.
To learn how Scilab helps in solving such models.
List of equipment/software
Personal Computer
Scilab
Deliverables
Scilab scripts and their results for all the assignments and exercises
properly discussed and explained.
Analytical conclusion for this lab activity.
Name: Ceferino Kevin A. Tan
EE 179.1 Section: M45
Course and Year: BSECE - 4
Laboratory Schedule: Monday 12-3
Name: Ceferino Kevin A. Tan
EE 179.1 Section: M45
Course and Year: BSECE - 4
Laboratory Schedule: Monday 12-3
Modeling Physical Systems
4.1
Mass-Spring System Model
Consider the following Mass -Spring system shown in Figure 1 where K is the
spring force, fv is the friction coefficient, x(t) is the displacement and f(t) is
the applied force:
Figure 1. Mass-Spring System and its Free-Body
Diagram
Referring to the free-body diagram, we get
M
d2 x ( t )
dx ( t )
+f v
+ Kx(t)=f (t )
2
dt
dt
(1)
The second order linear differential equation (1) describes the relationship
between the displacement and the applied force. The differential equation
can then be used to study the time behavior of x(t) under various changes
of the applied force.
The objectives behind modeling the mass-damper system can be many and
may include:
Understanding the dynamics of such system.
Name: Ceferino Kevin A. Tan
EE 179.1 Section: M45
Course and Year: BSECE - 4
Laboratory Schedule: Monday 12-3
Studying the effect of each parameter on the system such as mass M,
the friction coefficient B, and the elastic characteristic Kx(t).
Designing a new component such as damper or spring.
Reproducing a problem in order to suggest a solution.
4.2
Using Scilab in Solving Ordinary Differential Equations
Scilab can help solve linear or nonlinear ordinary differential equations (ODE)
using ode tool. To show how you can solve ODE using Scilab, we will proceed
in two ways. We first see how we can solve a first order ODE and then see
how we can solve a second order ODE.
4. 3 First Order ODE: Speed Cruise Control Example
Assume a zero spring force which means that K = 0. Equation (1) becomes
M
d2 x ( t )
dx ( t )
+f v
=f (t)
2
dt
dt
(2)
or
M
dv ( t )
+f v v (t )=f (t )
dt
(3)
since
2
a(t )=
and
dv ( t ) d x (t )
=
dt
d t2
v (t )=
,
dx ( t )
dt
Equation (3) is a first order linear ODE and we can use Scilab to solve for this
differential equation.
We can model and solve for equation (3) by writing the following script in
Scilab. You can also input the codes directly into the Scilab Console.
Code:
Name: Ceferino Kevin A. Tan
EE 179.1 Section: M45
Course and Year: BSECE - 4
Laboratory Schedule: Monday 12-3
//declare constant values
M = 750
B = 30
Fa = 300
//the force f(t)
//differential equation definition
function dvdt = f(t,v)
dvdt = (Fa-B*v)/M;
endfunction
//initial conditions @ t=0, v=0
t0 = 0
v0 = 0
t = 0:0.1:5;
v = ode(v0, t0, t, f);
clf;
plot2d(t, v)
//read references for ode tool
Name: Ceferino Kevin A. Tan
EE 179.1 Section: M45
Course and Year: BSECE - 4
Laboratory Schedule: Monday 12-3
4. 4 Second Order ODE: Mass-Spring System Example
In reality, the spring force and/or friction force can have a more complicated
expression or could be represented by a graph or data table. For instance, a
nonlinear spring can be designed such that the elastic characteristic is
K xr (t)
where r > 1. Figure 2 is an example of a nonlinear spring.
Figure 2. Nonlinear Spring[1]
In
such
case,
equation
(1)
becomes
d x (t)
dx ( t )
+f v
+ K x r ( t)=f (t)
2
dt
dt
(4)
Equation (3) represents another possible model that describes the dynamic
behavior of the mass-damper system under external force. Unlike equation
(1), this is said to be a nonlinear differential equation.
We can also solve equation (3) using Scilab ode tool but this time, the second
order differential equation has to be decomposed in a set of first order
differential equations as follows:
Name: Ceferino Kevin A. Tan
EE 179.1 Section: M45
Let
x ( t )=X 1
d X1
= X2
dt
so
d X 2 f v
f (t )
K
=
X 2 X 1r +
dt
M
M
M
so
In vector form, let
Course and Year: BSECE - 4
Laboratory Schedule: Monday 12-3
[ ][]
X=
X1
x
= '
X2
x
dx(t)
=X 2 ,
dt
and
[ ][ ]
d X1
dX
dt
x'
=
=
dt
d X2
x' '
dt
so we get
X2
dX
= f v
f (t )
K
dt
X 2 X 1 r +
M
M
M
The Scilab ode tool can now be used:
Code:
//declare constant values
M = 750
B = 30
Fa = 300
//the force f(t)
K = 15
r = 1
//differential equation definition
function xprime = f(t,x)
xprime(1) = x(2);
xprime(2) = -(B/M)*x(2)-(K/M)*((x(1))^r)+(Fa/M);
endfunction
t0 = 0
t = 0:0.1:5;
x0 = 0; //set initial position
xprime0 = 0; //set initial velocity
x = ode([x0; xprime0], t0, t, f);
Name: Ceferino Kevin A. Tan
EE 179.1 Section: M45
Course and Year: BSECE - 4
Laboratory Schedule: Monday 12-3
clf;
plot2d(t, x(1,:))
// x(1,:) plots the position vs time
// x(n, :) to plot a solution in matrix
form
// with dimension n x T
Name: Ceferino Kevin A. Tan
EE 179.1 Section: M45
Course and Year: BSECE - 4
Laboratory Schedule: Monday 12-3
Assessment
5.1
Assignment
1. Show and discuss the graphs obtained from the sample simulations.
The first graph shows
the
velocity-time
relation of an object,
having a mass of 750
kg, being influenced
by a force of Fa =
300N moving at one
direction,
all
while
under the coefficient
of
friction
of
magnitude 30 N/m.
The objects given
initial condition is vo = 0 m/s at to = 0 s.
At first glance, the function appears to
be linear, when in fact, it is not.
We used the ode function (ordinary
differential
equation).
Ode
differential equations defined by:
dy
=f ( t , y )
dt
solves
Name: Ceferino Kevin A. Tan
EE 179.1 Section: M45
Course and Year: BSECE - 4
Laboratory Schedule: Monday 12-3
y (t 0)= y 0
The second sample code yields the graph
shown to the right.
This graph plots the position of any point
in the spring with respect to time.
Notice that the rate by which the position
changes over time increases, hence the
upward parabolic characteristic of the
graph.
Screenshots
Sample code 1
Name: Ceferino Kevin A. Tan
EE 179.1 Section: M45
Course and Year: BSECE - 4
Laboratory Schedule: Monday 12-3
Sample code 2
Name: Ceferino Kevin A. Tan
EE 179.1 Section: M45
5.2
Course and Year: BSECE - 4
Laboratory Schedule: Monday 12-3
Exercise 1
Referring to the mass-spring system example:
1. Plot the position and the speed of the system both with respect to time
in separate graphs.
Graph of position vs time
Graph of velocity vs time
2. Change the value of r to 2 and 3. Plot the position and speed both with
respect to time. Discuss the results.
When r=2
Position vs. Time
Name: Ceferino Kevin A. Tan
EE 179.1 Section: M45
Course and Year: BSECE - 4
Laboratory Schedule: Monday 12-3
When r=3
Velocity vs. Time
Position
It can be said
that vs.
Velocity vs.
Time
Time
Name: Ceferino Kevin A. Tan
EE 179.1 Section: M45
Course and Year: BSECE - 4
Laboratory Schedule: Monday 12-3
3. With r = 1, vary the value of K (multiply by 5 and 10) and discuss the
results.
With K = 75
With K = 150
A higher spring constant will result in a bouncier spring. That is, a
higher spring constant results in an increased frequency of
oscillation of the spring over a period of time. It is also noticeable
that the range of motion has decreased for increasing values of the
spring constant, K.
Name: Ceferino Kevin A. Tan
EE 179.1 Section: M45
5.3
Course and Year: BSECE - 4
Laboratory Schedule: Monday 12-3
Exercise 2
Consider the mechanical system depicted in the
The input is given by f(t), and the output
figure.
is
given by y(t). Determine the differential
equation governing the system and using Scilab,
write a script and plot the system response such
that forcing function f(t) = 1. Let m = 10, k
= 0.5.
a.
b.
c.
d.
The peak displacement of the block is about ____.
At time t = ____, the output reaches peak displacement.
If k = 10, the peak value of y is ____.
What happens to the system when b is changed?
= 1, and b
Name: Ceferino Kevin A. Tan
EE 179.1 Section: M45
Course and Year: BSECE - 4
Laboratory Schedule: Monday 12-3
References
CISE 302 Linear Control Systems Laboratory Manual. (2011, October).
Systems Engineering Department, King Fahd University of Petroleum &
Minerals.
Scilab Enterprises . (2013). Scilab for very beginners. 78000 Versailles,
France.
Sengupta, A. (15, April 2010). Introduction to ODEs in Scilab. Indian Institute
of Technology Bombay.
